"Fossies" - the Fresh Open Source Software Archive

Member "PDL-2.080/Libtmp/Transform/Cartography/Cartography.pm" (7 May 2022, 89944 Bytes) of package /linux/misc/PDL-2.080.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Perl source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "Cartography.pm" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.079_vs_2.080.

    1 =head1 NAME
    2 
    3 PDL::Transform::Cartography - Useful cartographic projections
    4 
    5 =head1 SYNOPSIS
    6 
    7  # make a Mercator map of Earth
    8  use PDL::Transform::Cartography;
    9  $x = earth_coast();
   10  $x = graticule(10,2)->glue(1,$x);
   11  $t = t_mercator;
   12  $w = pgwin(xs);
   13  $w->lines($t->apply($x)->clean_lines());
   14 
   15 =head1 DESCRIPTION
   16 
   17 PDL::Transform::Cartography includes a variety of useful cartographic
   18 and observing projections (mappings of the surface of a sphere),
   19 including reprojected observer coordinates.  See L<PDL::Transform>
   20 for more information about image transforms in general.
   21 
   22 Cartographic transformations are used for projecting not just
   23 terrestrial maps, but also any nearly spherical surface including the
   24 Sun, the Celestial sphere, various moons and planets, distant stars,
   25 etc.  They also are useful for interpreting scientific images, which
   26 are themselves generally projections of a sphere onto a flat focal
   27 plane (e.g. the L</t_gnomonic> projection).
   28 
   29 Unless otherwise noted, all the transformations in this file convert
   30 from (theta,phi) coordinates on the unit sphere (e.g. (lon,lat) on a
   31 planet or (RA,dec) on the celestial sphere) into some sort of
   32 projected coordinates, and have inverse transformations that convert
   33 back to (theta,phi).  This is equivalent to working from the
   34 equidistant cylindrical (or L<"plate carree"|/t_carree>) projection, if
   35 you are a cartography wonk.
   36 
   37 The projected coordinates are generally in units of body radii
   38 (radians), so that multiplying the output by the scale of the map
   39 yields physical units that are correct wherever the scale is correct
   40 for that projection.  For example, areas should be correct everywhere
   41 in the authalic projections; and linear scales are correct along
   42 meridians in the equidistant projections and along the standard
   43 parallels in all the projections.
   44 
   45 The transformations that are authalic (equal-area), conformal
   46 (equal-angle), azimuthal (circularly symmetric), or perspective (true
   47 perspective on a focal plane from some viewpoint) are marked.  The
   48 first two categories are mutually exclusive for all but the 
   49 L<"unit sphere"|/t_unit_sphere> 3-D projection.
   50 
   51 Extra dimensions tacked on to each point to be transformed are, in
   52 general, ignored.  That is so that you can add on an extra index
   53 to keep track of pen color.  For example, L</earth_coast>
   54 returns a 3x<n> ndarray containing (lon, lat, pen) at each list location.
   55 Transforming the vector list retains the pen value as the first index
   56 after the dimensional directions.
   57 
   58 =head1 GENERAL NOTES ON CARTOGRAPHY
   59 
   60 Unless otherwise noted, the transformations and miscellaneous
   61 information in this section are taken from Snyder & Voxland 1989: "An
   62 Album of Map Projections", US Geological Survey Professional Paper
   63 1453, US Printing Office (Denver); and from Snyder 1987: "Map
   64 Projections - A Working Manual", US Geological Survey Professional
   65 Paper 1395, US Printing Office (Denver, USA).  You can obtain your own
   66 copy of both by contacting the U.S. Geological Survey, Federal Center,
   67 Box 25425, Denver, CO 80225 USA.
   68 
   69 The mathematics of cartography have a long history, and the details
   70 are far trickier than the broad overview.  For terrestrial (and, in
   71 general, planetary) cartography, the best reference datum is not a
   72 sphere but an oblate ellipsoid due to centrifugal force from the
   73 planet's rotation.  Furthermore, because all rocky planets, including
   74 Earth, have randomly placed mass concentrations that affect the
   75 gravitational field, the reference gravitational isosurface (sea level
   76 on Earth) is even more complex than an ellipsoid and, in general,
   77 different ellipsoids have been used for different locations at the
   78 same time and for the same location at different times.
   79 
   80 The transformations in this package use a spherical datum and hence
   81 include global distortion at about the 0.5% level for terrestrial maps
   82 (Earth's oblateness is ~1/300).  This is roughly equal to the
   83 dimensional precision of physical maps printed on paper (due to
   84 stretching and warping of the paper) but is significant at larger
   85 scales (e.g. for regional maps).  If you need more precision than
   86 that, you will want to implement and use the ellipsoidal
   87 transformations from Snyder 1987 or another reference work on geodesy.
   88 A good name for that package would be C<...::Cartography::Geodetic>.
   89 
   90 =head1 GENERAL NOTES ON PERSPECTIVE AND SCIENTIFIC IMAGES
   91 
   92 Cartographic transformations are useful for interpretation of
   93 scientific images, as all cameras produce projections of the celestial
   94 sphere onto the focal plane of the camera.  A simple (single-element)
   95 optical system with a planar focal plane generates
   96 L<gnomonic|/t_gnomonic> images -- that is to say, gnomonic projections
   97 of a portion of the celestial sphere near the paraxial direction.
   98 This is the projection that most consumer grade cameras produce.
   99 
  100 Magnification in an optical system changes the angle of incidence
  101 of the rays on the focal plane for a given angle of incidence at the
  102 aperture.  For example, a 10x telescope with a 2 degree field of view
  103 exhibits the same gnomonic distortion as a simple optical system with 
  104 a 20 degree field of view.  Wide-angle optics typically have magnification
  105 less than 1 ('fisheye lenses'), reducing the gnomonic distortion 
  106 considerably but introducing L<"equidistant azimuthal"|/t_az_eqd> distortion --
  107 there's no such thing as a free lunch!
  108 
  109 Because many solar-system objects are spherical,
  110 PDL::Transform::Cartography includes perspective projections for
  111 producing maps of spherical bodies from perspective views.  Those
  112 projections are L<"t_vertical"|/t_vertical> and
  113 L<"t_perspective"|/t_perspective>.  They map between (lat,lon) on the
  114 spherical body and planar projected coordinates at the viewpoint.  
  115 L<"t_vertical"|/t_vertical> is the vertical perspective projection 
  116 given by Snyder, but L<"t_perspective"|/t_perspective> is a fully
  117 general perspective projection that also handles magnification correction.
  118 
  119 =head1 TRANSVERSE & OBLIQUE PROJECTIONS; STANDARD OPTIONS
  120 
  121 Oblique projections rotate the sphere (and graticule) to an arbitrary
  122 angle before generating the projection; transverse projections rotate
  123 the sphere exactly 90 degrees before generating the projection.  
  124 
  125 Most of the projections accept the following standard options,
  126 useful for making transverse and oblique projection maps.  
  127 
  128 =over 3
  129 
  130 =item o, origin, Origin [default (0,0,0)]
  131 
  132 The origin of the oblique map coordinate system, in (old-theta, old-phi) 
  133 coordinates.
  134 
  135 =item r, roll, Roll [default 0.0]
  136 
  137 The roll angle of the sphere about the origin, measured CW from (N = up)
  138 for reasonable values of phi and CW from (S = up) for unreasonable
  139 values of phi.  This is equivalent to observer roll angle CCW from the
  140 same direction.
  141 
  142 =item u, unit, Unit [default 'degree']
  143 
  144 This is the name of the angular unit to use in the lon/lat coordinate system.
  145 
  146 =item b, B
  147 
  148 The "B" angle of the body -- used for extraterrestrial maps.  Setting
  149 this parameter is exactly equivalent to setting the phi component
  150 of the origin, and in fact overrides it.
  151 
  152 =item l,L
  153 
  154 The longitude of the central meridian as observed -- used for extraterrestrial
  155 maps.  Setting this parameter is exactly equivalent to setting the theta
  156 component of the origin, and in fact overrides it.
  157 
  158 =item p,P
  159 
  160 The "P" (or position) angle of the body -- used for extraterrestrial maps.
  161 This parameter is a synonym for the roll angle, above.
  162 
  163 =item bad, Bad, missing, Missing [default nan]
  164 
  165 This is the value that missing points get.  Mainly useful for the
  166 inverse transforms.  (This should work fine if set to BAD, if you have
  167 bad-value support compiled in).  The default nan is asin(1.2), calculated
  168 at load time.
  169 
  170 =back
  171 
  172 =head1 EXAMPLES
  173 
  174 Draw a Mercator map of the world on-screen:
  175 
  176    $w = pgwin(xs);
  177    $w->lines(earth_coast->apply(t_mercator)->clean_lines);
  178 
  179 Here, C<earth_coast()> returns a 3xn ndarray containing (lon, lat, pen) 
  180 values for the included world coastal outline; C<t_mercator> converts
  181 the values to projected Mercator coordinates, and C<clean_lines> breaks
  182 lines that cross the 180th meridian.
  183 
  184 Draw a Mercator map of the world, with lon/lat at 10 degree intervals:
  185 
  186    $w = pgwin(xs)
  187    $x = earth_coast()->glue(1,graticule(10,1));
  188    $w->lines($x->apply(t_mercator)->clean_lines);
  189 
  190 This works just the same as the first example, except that a map graticule
  191 has been applied with interline spacing of 10 degrees lon/lat and 
  192 inter-vertex spacing of 1 degree (so that each meridian contains 181 points,
  193 and each parallel contains 361 points).
  194 
  195 =head1 NOTES
  196 
  197 Currently angular conversions are rather simpleminded.  A list of
  198 common conversions is present in the main constructor, which inserts a
  199 conversion constant to radians into the {params} field of the new
  200 transform.  Something like Math::Convert::Units should be used instead
  201 to generate the conversion constant. 
  202 
  203 A cleaner higher-level interface is probably needed (see the examples);
  204 for example, earth_coast could return a graticule if asked, instead of 
  205 needing one to be glued on.
  206 
  207 The class structure is somewhat messy because of the varying needs of
  208 the different transformations.  PDL::Transform::Cartography is a base
  209 class that interprets the origin options and sets up the basic
  210 machinery of the Transform.  The conic projections have their
  211 own subclass, PDL::Transform::Conic, that interprets the standard
  212 parallels.  Since the cylindrical and azimuthal projections are pretty
  213 simple, they are not subclassed.
  214 
  215 The perl 5.6.1 compiler is quite slow at adding new classes to the
  216 structure, so it does not makes sense to subclass new transformations
  217 merely for the sake of pedantry.
  218 
  219 =head1 AUTHOR
  220 
  221 Copyright 2002, Craig DeForest (deforest@boulder.swri.edu).  This
  222 module may be modified and distributed under the same terms as PDL
  223 itself.  The module comes with NO WARRANTY.
  224 
  225 The included digital world map is derived from the 1987 CIA World Map,
  226 translated to ASCII in 1988 by Joe Dellinger (geojoe@freeusp.org) and
  227 simplified in 1995 by Kirk Johnson (tuna@indra.com) for the program
  228 XEarth.  The map comes with NO WARRANTY.  An ASCII version of the map,
  229 and a sample PDL function to read it, may be found in the Demos
  230 subdirectory of the PDL source distribution.
  231 
  232 =head1 FUNCTIONS
  233 
  234 The module exports both transform constructors ('t_<foo>') and some
  235 auxiliary functions (no leading 't_').
  236 
  237 =cut
  238 
  239 # Import PDL::Transform into the calling package -- the cartography
  240 # stuff isn't much use without it.
  241 use PDL::Transform;
  242 
  243 package PDL::Transform::Cartography;
  244 use strict;
  245 use warnings;
  246 use PDL::Core ':Internal'; # Load 'topdl' (internal routine)
  247 
  248 use Exporter ();
  249 use PDL::LiteF;
  250 use PDL::Math;
  251 use PDL::Transform;
  252 use PDL::MatrixOps;
  253 use Carp;
  254 
  255 our @ISA = ( 'Exporter','PDL::Transform' );
  256 our $VERSION = "0.6";
  257 $VERSION = eval $VERSION;
  258 our @EXPORT_OK = qw(
  259   graticule earth_image earth_coast earth_shape clean_lines t_unit_sphere
  260   t_orthographic t_rot_sphere t_caree t_carree t_mercator t_utm t_sin_lat
  261   t_sinusoidal t_conic t_albers t_lambert t_stereographic t_gnomonic
  262   t_az_eqd t_az_eqa t_vertical t_perspective t_hammer t_aitoff
  263   t_raster2fits t_raster2float
  264 );
  265 our @EXPORT = @EXPORT_OK;
  266 our %EXPORT_TAGS = (Func=>\@EXPORT_OK);
  267 
  268 ##############################
  269 # Steal _opt from PDL::Transform.
  270 *PDL::Transform::Cartography::_opt = \&PDL::Transform::_opt;
  271 use overload '""' => \&_strval;
  272 
  273 our $PI = $PDL::Transform::PI;
  274 our $DEG2RAD = $PDL::Transform::DEG2RAD;
  275 our $RAD2DEG = $PDL::Transform::RAD2DEG;
  276 
  277 sub _strval {
  278   my($me) = shift;
  279   $me->stringify();
  280 }
  281 
  282 ######################################################################
  283 
  284 =head2 graticule
  285 
  286 =for usage
  287 
  288    $lonlatp     = graticule(<grid-spacing>,<line-segment-size>);   
  289    $lonlatp     = graticule(<grid-spacing>,<line-segment-size>,1);
  290 
  291 =for ref
  292 
  293 (Cartography) PDL constructor - generate a lat/lon grid.
  294 
  295 Returns a grid of meridians and parallels as a list of vectors suitable
  296 for sending to
  297 L<PDL::Graphics::PGPLOT::Window::lines|PDL::Graphics::PGPLOT::Window/lines>
  298 for plotting.
  299 The grid is in degrees in (theta, phi) coordinates -- this is (E lon, N lat) 
  300 for terrestrial grids or (RA, dec) for celestial ones.  You must then 
  301 transform the graticule in the same way that you transform the map.
  302 
  303 You can attach the graticule to a vector map using the syntax:
  304 
  305     $out = graticule(10,2)->glue(1,$map);
  306 
  307 In array context you get back a 2-element list containing an ndarray of
  308 the (theta,phi) pairs and an ndarray of the pen values (1 or 0) suitable for
  309 calling
  310 L<PDL::Graphics::PGPLOT::Window::lines|PDL::Graphics::PGPLOT::Window/lines>.
  311 In scalar context the two elements are combined into a single ndarray.
  312 
  313 The pen values associated with the graticule are negative, which will cause
  314 L<PDL::Graphics::PGPLOT::Window::lines|PDL::Graphics::PGPLOT::Window/lines>
  315 to plot them as hairlines.
  316 
  317 If a third argument is given, it is a hash of options, which can be:
  318 
  319 =over 3
  320 
  321 =item nan - if true, use two columns instead of three, and separate lines with a 'nan' break
  322 
  323 =item lonpos - if true, all reported longitudes are positive (0 to 360) instead of (-180 to 180).
  324 
  325 =item dup - if true, the meridian at the far boundary is duplicated.
  326 
  327 =back 
  328 
  329 =cut
  330 
  331 sub graticule {
  332     my $grid = shift;
  333     my $step = shift;
  334     my $hash = shift; $hash = {} unless defined($hash); # avoid // for ancient compatibility
  335     my $two_cols = $hash->{nan} || 0;
  336     my $lonpos = $hash->{lonpos} || 0;
  337     my $dup = $hash->{dup} || 0;
  338 
  339     $grid = 10 unless defined($grid);
  340     $grid = $grid->at(0) if(ref($grid) eq 'PDL');
  341     
  342     $step = $grid/2 unless defined($step);
  343     $step = $step->at(0) if(ref($step) eq 'PDL');
  344 
  345     # Figure number of parallels and meridians
  346     my $np = 2 * int(90/$grid);
  347     my $nm = 2 * int(180/$grid);
  348     
  349     # First do parallels.
  350     my $xp = xvals(360/$step + 1 + !!$two_cols, $np + 1) * $step       - 180 * (!$lonpos);
  351     my $yp = yvals(360/$step + 1 + !!$two_cols, $np + 1) * 180/$np     -  90;
  352     $xp->slice("-1,:") .= $yp->slice("-1,:") .= asin(pdl(1.1)) if($two_cols);
  353 
  354     # Next do meridians.
  355     my $xm = yvals( 180/$step + 1 + !!$two_cols,  $nm + !!$dup  ) * 360/$nm   - 180 * (!$lonpos);
  356     my $ym = xvals( 180/$step + 1 + !!$two_cols,  $nm + !!$dup  ) * $step     - 90;
  357     $xm->slice("-1,:") .= $ym->slice("-1,:") .= asin(pdl(1.1)) if($two_cols);
  358 
  359     
  360     if($two_cols) {
  361     return pdl( $xp->flat->append($xm->flat),
  362             $yp->flat->append($ym->flat)
  363         )->mv(1,0);
  364     } else {
  365      our $pp = (zeroes($xp)-1); $pp->slice("(-1)") .= 0;
  366      our $pm = (zeroes($xm)-1); $pm->slice("(-1)") .= 0;
  367 
  368     if(wantarray) {
  369         return (  pdl( $xp->flat->append($xm->flat),
  370                $yp->flat->append($ym->flat)
  371               )->mv(1,0),
  372               
  373               $pp->flat->append($pm->flat)
  374         );
  375     } else {
  376         return pdl( $xp->flat->append($xm->flat),
  377             $yp->flat->append($ym->flat),
  378             $pp->flat->append($pm->flat)
  379              )->mv(1,0);
  380     }  
  381     barf "This can't happen";
  382     }
  383 }
  384 
  385 =head2 earth_coast
  386 
  387 =for usage
  388 
  389   $x = earth_coast()
  390 
  391 =for ref
  392 
  393 (Cartography) PDL constructor - coastline map of Earth
  394 
  395 Returns a vector coastline map based on the 1987 CIA World Coastline
  396 database (see author information).  The vector coastline data are in
  397 plate carree format so they can be converted to other projections via
  398 the L<apply|PDL::Transform/apply> method and cartographic transforms,
  399 and are suitable for plotting with the
  400 L<lines|PDL::Graphics::PGPLOT::Window/lines>
  401 method in the PGPLOT
  402 output library:  the first dimension is (X,Y,pen) with breaks having 
  403 a pen value of 0 and hairlines having negative pen values.  The second 
  404 dimension broadcasts over all the points in the data set.
  405 
  406 The vector map includes lines that pass through the antipodean
  407 meridian, so if you want to plot it without reprojecting, you should
  408 run it through L</clean_lines> first:
  409 
  410     $w = pgwin();
  411     $w->lines(earth_coast->clean_lines);     # plot plate carree map of world
  412     $w->lines(earth_coast->apply(t_gnomonic))# plot gnomonic map of world
  413 
  414 C<earth_coast> is just a quick-and-dirty way of loading the file
  415 "earth_coast.vec.fits" that is part of the normal installation tree.
  416 
  417 =cut
  418 
  419 sub earth_coast {
  420     my $fn = "PDL/Transform/Cartography/earth_coast.vec.fits";
  421     local $_;
  422     require PDL::IO::FITS;
  423     foreach(@INC) {
  424     my $file = "$_/$fn";
  425     return PDL::IO::FITS::rfits($file) if(-e $file);
  426     }
  427     barf("earth_coast: $fn not found in \@INC.\n");
  428 }
  429 
  430 =head2 earth_image
  431 
  432 =for usage
  433 
  434  $rgb = earth_image()
  435 
  436 =for ref
  437 
  438 (Cartography) PDL constructor - RGB pixel map of Earth 
  439 
  440 Returns an RGB image of Earth based on data from the MODIS instrument
  441 on the NASA EOS/Terra satellite.  (You can get a full-resolution
  442 image from L<http://earthobservatory.nasa.gov/Newsroom/BlueMarble/>).
  443 The image is a plate carree map, so you can convert it to other
  444 projections via the L<map|PDL::Transform/map> method and cartographic
  445 transforms.
  446 
  447 This is just a quick-and-dirty way of loading the earth-image files that
  448 are distributed along with PDL.
  449 
  450 =cut
  451 
  452 sub earth_image {
  453   my($nd) = shift;
  454   my $f;
  455   require PDL::IO::Pic;
  456   my $dir = "PDL/Transform/Cartography/earth_";
  457   $f = (($nd//'') =~ m/^n/i) ? "${dir}night.jpg" : "${dir}day.jpg";
  458   local $_;
  459   my $im;
  460   my $found = 0;
  461   foreach(@INC) {
  462     my $file = "$_/$f";
  463     if(-e $file) {
  464       $found = 1;
  465       $im = PDL::IO::Pic::rpic($file);
  466     }
  467     last if defined($im);
  468   }
  469   barf("earth_image: $f not found in \@INC\n")
  470     unless defined($found);
  471   barf("earth_image: couldn't load $f; you may need to install netpbm.\n")
  472     unless defined($im);
  473   t_raster2fits()->apply($im);
  474 }
  475 
  476 =head2 earth_shape
  477 
  478 =for usage
  479 
  480  $fits_shape = earth_shape()
  481 
  482 =for ref
  483 
  484 (Cartography) PDL constructor - height map of Earth
  485 
  486 Returns a height map of Earth based on data from the General
  487 Bathymetric Chart of the Oceans (GEBCO) produced by the British
  488 Oceanographic Data Centre. (You can get a full-resolution image from
  489 L<http://visibleearth.nasa.gov/view.php?id=73934>).  The image is a
  490 plate carree map, so you can convert it to other projections via the
  491 L<map|PDL::Transform/map> method and cartographic transforms.
  492 The data is from 8-bit grayscale (so only 256 levels), but is returned
  493 in a similar format to L</earth_image>. The range represents a span of
  494 6400m, so Everest and the Marianas Trench are not accurately represented.
  495 
  496 To turn this into a C<float>, (C<lonlatradius,x,y>) with C<x>
  497 and C<y> in radians, and the radius as a C<float> as a proportion of the
  498 Earth's mean radius, use L</t_raster2float>.
  499 The Earth is treated here as a perfect sphere with sea
  500 level at radius 6,371km.
  501 
  502   Value       Hex value   Float    From centre in km   Float as radius
  503   Base        00          0.0      6370.69873km        0.99995
  504   Sea level   0C          0.04705  6371km              1.0
  505   Highest     FF          1.0      6377.09863km        1.00096
  506 
  507 Code:
  508 
  509   $shape = earth_shape();
  510   $floats = t_raster2float()->apply($shape->mv(2,0));
  511   $lonlatradius = $floats->slice('0:2'); # r g b all same
  512   $lonlatradius->slice('(2)') *= float((6377.09863 - 6370.69873) / 6371);
  513   $lonlatradius->slice('(2)') += float(6370.69873 / 6371);
  514 
  515 =cut
  516 
  517 sub earth_shape {
  518   my($nd) = shift;
  519   require PDL::IO::Pic;
  520   my $f = "PDL/Transform/Cartography/earth_height-2048x1024.jpg";
  521   local $_;
  522   my $im;
  523   my $found = 0;
  524   foreach (@INC) {
  525     my $file = "$_/$f";
  526     if(-e $file) {
  527       $found = 1;
  528       $im = PDL::IO::Pic::rpic($file);
  529     }
  530     last if defined($im);
  531   }
  532   barf("earth_shape: $f not found in \@INC\n")
  533     unless defined($found);
  534   barf("earth_shape: couldn't load $f; you may need to install netpbm.\n")
  535     unless defined($im);
  536   $im = $im->dummy(0,3); # fake RGB
  537   t_raster2fits()->apply($im);
  538 }
  539 
  540 =head2 clean_lines
  541 
  542 =for usage
  543 
  544  $x = clean_lines(t_mercator->apply(scalar(earth_coast())));
  545  $x = clean_lines($line_pen, [threshold]);
  546  $x = $lines->clean_lines;
  547 
  548 =for ref
  549 
  550 (Cartography) PDL method - remove projection irregularities
  551 
  552 C<clean_lines> massages vector data to remove jumps due to singularities
  553 in the transform.
  554 
  555 In the first (scalar) form, C<$line_pen> contains both (X,Y) points and pen 
  556 values suitable to be fed to
  557 L<lines|PDL::Graphics::PGPLOT::Window/lines>:
  558 in the second (list) form, C<$lines> contains the (X,Y) points and C<$pen>
  559 contains the pen values.  
  560 
  561 C<clean_lines> assumes that all the outline polylines are local --
  562 that is to say, there are no large jumps.  Any jumps larger than a
  563 threshold size are broken by setting the appropriate pen values to 0.
  564 
  565 The C<threshold> parameter sets the relative size of the largest jump, relative
  566 to the map range (as determined by a min/max operation).  The default size is
  567 0.1.
  568 
  569 NOTES
  570 
  571 This almost never catches stuff near the apex of cylindrical maps,
  572 because the anomalous vectors get arbitrarily small.  This could be 
  573 improved somewhat by looking at individual runs of the pen and using
  574 a relative length scale that is calibrated to the rest of each run.
  575 it is probably not worth the computational overhead.
  576 
  577 =cut
  578 
  579 *PDL::clean_lines = *PDL::clean_lines = \&clean_lines;
  580 sub clean_lines {
  581     my($lines) = shift;
  582     my($x) = shift;
  583     my($y) = shift;
  584     my($l,$p,$th);
  585 
  586     $th = 0.1;
  587 
  588     if(defined($y)) {
  589     # separate case with thresh
  590     $l = $lines;
  591     $p = $x->is_inplace?$x:$x->copy;
  592     $th = $y;
  593     } else {
  594     if(!defined($x)) {
  595         # duplex case no thresh
  596         $l = $lines->slice("0:1");
  597         $p = $lines->is_inplace ? $lines->slice("(2)") : $lines->slice("(2)")->sever;
  598     } elsif(UNIVERSAL::isa($x,'PDL') && 
  599         $lines->slice("(0)")->nelem == $x->nelem) {
  600         # Separate case no thresh
  601         $l = $lines;
  602         $p = $x->is_inplace ? $x : $x->copy;;
  603     } else {
  604         # duplex case with thresh
  605         $l = $lines->slice("0:1");
  606         $p = $lines->is_inplace ? $lines->slice("(2)") : $lines->slice("(2)")->sever;
  607         $th = $x;
  608     }
  609     }
  610 
  611     my $pok = (($p != 0) & isfinite($p));
  612     # Kludge to work around minmax bug (nans confuse it!)
  613     my($l0) = $l->slice("(0)");
  614     my($x0,$x1) = $l0->where(isfinite($l0) & $pok)->minmax;
  615     my($xth) = abs($x1-$x0) * $th;
  616 
  617     my($l1) = $l->slice("(1)");
  618     ($x0,$x1) = $l1->where(isfinite($l1) & $pok)->minmax;
  619     my($yth) = abs($x1-$x0) * $th;
  620 
  621     my $diff = abs($l->slice(":,1:-1") - $l->slice(":,0:-2"));
  622 
  623     $diff->where(!isfinite($diff)) .= 2*($xth + $yth); 
  624     $p->where(($diff->slice("(0)") > $xth) | ($diff->slice("(1)") > $yth)) .= 0;
  625     if(wantarray){
  626     return($l,$p);
  627     } else {
  628     return $l->append($p->dummy(0,1));
  629     }
  630 }    
  631 
  632 
  633 
  634 ######################################################################
  635 
  636 ###
  637 # Units parser
  638 # Get unit, return conversion factor to radii, or undef if no match found.
  639 #
  640 sub _uconv{
  641   ###
  642   # Replace this with a more general units resolver call!
  643   ###
  644   local($_) = shift;
  645   my($silent) =shift;
  646 
  647   my($x) = 
  648     ( m/^deg/i    ? $DEG2RAD :
  649       m/^arcmin/i ? $DEG2RAD / 60 :
  650       m/^arcsec/i ? $DEG2RAD / 3600 :
  651       m/^hour/i   ? $DEG2RAD * 15 :    # Right ascension
  652       m/^min/i    ? $DEG2RAD * 15/60 : # Right ascension
  653       m/^microrad/i ? 1e-6 :
  654       m/^millirad/i ? 1e-3 :
  655       m/^rad(ian)?s?$/i ? 1.0 :
  656       m/^meter/ ? 1.0/6371000 :                # Assuming Earth cartography!
  657       m/^kilometer/ ? 1.0/6371 :
  658       m/^km/ ? 1.0/6371 :
  659       m/^Mm/ ? 1.0/6.371 :
  660       m/^mile/  ? 1.0/(637100000/2.54/12/5280) :
  661       undef
  662       );
  663   print STDERR "Cartography: unrecognized unit '$_'\n"    
  664     if( (!defined $x) && !$silent && ($PDL::debug || $PDL::verbose));
  665   $x;
  666 }
  667 
  668 ###
  669 #
  670 # Cartography general constructor -- called by the individual map
  671 # constructors.  Not underscored because it's certainly OK to call from
  672 # outside -- but the *last* argument is the name of the transform.
  673 #
  674 # The options list is put into the {options} field of the newly constructed
  675 # Transform -- fastidious subclass constructors will want to delete it before 
  676 # returning.
  677 #
  678 
  679 
  680 sub _new { __PACKAGE__->new(@_); } # not exported
  681 sub new {
  682     my($class) = shift;
  683     my($name) = pop;
  684 
  685     my($o) = $_[0];
  686     $o = {@_}
  687       unless(ref $o eq 'HASH');
  688 
  689     my($me) = __PACKAGE__->SUPER::new;
  690     $me->{idim} = $me->{odim} = 2;
  691     $me->{name} = $name;
  692  
  693     ####
  694     # Parse origin and units arguments
  695     # 
  696     my $or = _opt($o,['o','origin','Origin'],zeroes(2));
  697     if($or->nelem != 2) {
  698     croak("PDL::Transform::Cartography: origin must have 2 elements\n");
  699     }
  700 
  701     my($l) = _opt($o,['l','L']);
  702     my($b_angle) = _opt($o,['b','B']);
  703 
  704     $or->slice("0") .= $l if defined($l);
  705     $or->slice("1") .= $b_angle if defined($b_angle);
  706 
  707     my $roll = topdl(_opt($o,['r','roll','Roll','P'],0));
  708     my $unit = _opt($o,['u','unit','Unit'],'degrees');
  709 
  710     $me->{params}->{conv} = my $conv = _uconv($unit);
  711     $me->{params}->{u} = $unit;
  712 
  713     $me->{itype} = ['longitude','latitude'];
  714     $me->{iunit} = [$me->{params}->{u},$me->{params}->{u}];
  715 
  716     my($ou) = _opt($o,['ou','ounit','OutputUnit'],undef);
  717     $me->{params}->{ou} = $ou;
  718     if(defined $ou) {
  719       if(!(ref $ou)) {
  720     $me->{params}->{oconv} = _uconv($ou);
  721       } else {
  722     my @oconv;
  723     map {push(@oconv,_uconv($_))} @$ou;
  724     $me->{params}->{oconv} = topdl([@oconv]);
  725       }
  726     } else {
  727       $me->{params}->{oconv} = undef;
  728     }
  729     $me->{ounit} = $me->{params}->{ou};
  730 
  731     $me->{params}->{o} = $or * $conv;
  732     $me->{params}->{roll} = $roll * $conv;
  733 
  734     $me->{params}->{bad} = _opt($o,['b','bad','Bad','missing','Missing'],
  735                   asin(pdl(1.1)));
  736 
  737     # Get the standard parallel (in general there's only one; the conics
  738     # have two but that's handled by _c_new)
  739     $me->{params}->{std} = topdl(_opt($me->{options},
  740              ['s','std','standard','Standard'],
  741              0))->at(0) * $me->{params}->{conv};
  742 
  743     $me->{options} = $o;
  744     $me;
  745 }
  746 
  747 # Compose self with t_rot_sphere if necessary -- useful for 
  748 # finishing off the transformations that accept the origin and roll 
  749 # options.
  750 sub PDL::Transform::Cartography::_finish {
  751   my($me) = shift;
  752   if( ($me->{params}->{o}->slice("0") != 0) ||
  753       ($me->{params}->{o}->slice("1") != 0) ||
  754       ($me->{params}->{roll} != 0) 
  755       ) {
  756       my $out = t_compose($me,t_rot_sphere($me->{options}));
  757       $out->{itype} = $me->{itype};
  758       $out->{iunit} = $me->{iunit};
  759       $out->{otype} = $me->{otype};
  760       $out->{ounit} = $me->{ounit};
  761       $out->{odim} = 2;
  762       $out->{idim} = 2;
  763       return $out;
  764     } 
  765   return $me;
  766 }
  767 
  768 =head2 t_raster2float
  769 
  770 =for usage
  771 
  772   $t = t_raster2float();
  773 
  774 =for ref
  775 
  776 (Cartography) Convert a raster (3,x,y) to C<float> (lonlatrgb,x,y)
  777 
  778 Assumes C<bytes> input, and radians and C<float> output, with the first
  779 2 coordinates suitable for use as plate carree.
  780 
  781 =cut
  782 
  783 sub t_raster2float {
  784   my ($me) = _new(@_, 'Raster FITS plate carree to OpenGL-ready float conversion');
  785   $me->{odim} = 3;
  786   $me->{params}->{itype} = ['RGB','X','Y'];
  787   $me->{params}->{iunit} = ['index','pixels','pixels'];
  788   $me->{odim} = 2;
  789   $me->{params}->{otype} = ['LonLatRGB','X','Y'];
  790   $me->{params}->{ounit} = ['Float','index','index'];
  791   $me->{func} = sub {
  792     my($d,$o) = @_;
  793     my (undef, $x, $y, @otherdims) = $d->dims;
  794     my $type = float;
  795     my $out_xy = zeroes(byte, $x, $y, @otherdims);
  796     my $out = zeroes($type, 5, $x, $y, @otherdims);
  797     $out->slice($_->[0]) .= $_->[1]
  798       for ['(0)', $out_xy->xlinvals(-$PI, $PI)],
  799         ['(1)', $out_xy->ylinvals(-$PI/2, $PI/2)],
  800         ['2:4', $d->convert($type) / 255];
  801     $out;
  802   };
  803   $me->{inv} = sub {
  804     my($d,$o) = @_;
  805     my $type = byte;
  806     my (undef, $x, $y, @otherdims) = $d->dims;
  807     my $out = zeroes($type, 3, $x, $y, @otherdims);
  808     $out .= $d->slice('2:4') * 255;
  809     $out;
  810   };
  811   $me;
  812 }
  813 
  814 =head2 t_raster2fits
  815 
  816 =for usage
  817 
  818   $t = t_raster2fits();
  819 
  820 =for ref
  821 
  822 (Cartography) Convert a raster (3,x,y) to FITS plate carree (x,y,3)
  823 
  824 Adds suitable C<hdr>. Assumes degrees. Used by L</earth_image>.
  825 
  826 =cut
  827 
  828 sub t_raster2fits {
  829   my ($me) = _new(@_, 'Raster to FITS plate carree conversion');
  830   $me->{odim} = 3;
  831   $me->{params}->{itype} = ['RGB','X','Y'];
  832   $me->{params}->{iunit} = ['RGB','pixels','pixels'];
  833   $me->{params}->{otype} = ['X','Y','RGB'];
  834   $me->{params}->{ounit} = ['pixels','pixels','RGB'];
  835   $me->{func} = sub {
  836     my($d,$o) = @_;
  837     my $out = $d->mv(0,2);
  838     my $h = $out->fhdr;
  839     $h->{SIMPLE} = 'T';
  840     $h->{NAXIS} = $d->ndims;
  841     local $_;
  842     $h->{"NAXIS".($_+1)} = $out->dim($_) for 0..$out->ndims-1;
  843     $h->{"CRVAL".($_+1)} = 0 for 0..$out->ndims-1;
  844     $h->{"CRPIX".($_+1)} = $_<2 ? ($out->dim($_)+1)/2 : 1 for 0..$out->ndims-1;
  845     my ($lon, $lat) = $out->dims;
  846     $h->{CTYPE1}='Longitude';   $h->{CUNIT1}='degrees'; $h->{CDELT1}=360/$lon;
  847     $h->{CTYPE2}='Latitude';    $h->{CUNIT2}='degrees'; $h->{CDELT2}=180/$lat;
  848     $h->{CTYPE3}='RGB';         $h->{CUNIT3}='index';   $h->{CDELT3}=1.0;
  849     $h->{COMMENT}='Plate Carree Projection';
  850     $h->{HISTORY}='PDL conversion from raster image',
  851     $out->hdrcpy(1);
  852     $out;
  853   };
  854   $me->{inv} = sub {
  855     my($d,$o) = @_;
  856     $d->mv(2,0);
  857   };
  858   $me;
  859 }
  860 
  861 ######################################################################
  862 
  863 =head2 t_unit_sphere
  864 
  865 =for usage
  866 
  867   $t = t_unit_sphere(<options>);
  868 
  869 =for ref
  870 
  871 (Cartography) 3-D globe projection (conformal; authalic)
  872 
  873 This is similar to the inverse of
  874 L<t_spherical|PDL::Transform/t_spherical>,
  875 but the
  876 inverse transform projects 3-D coordinates onto the unit sphere,
  877 yielding only a 2-D (lon/lat) output.  Similarly, the forward
  878 transform deprojects 2-D (lon/lat) coordinates onto the surface of a
  879 unit sphere.
  880 
  881 The cartesian system has its Z axis pointing through the pole of the 
  882 (lon,lat) system, and its X axis pointing through the equator at the 
  883 prime meridian.
  884 
  885 Unit sphere mapping is unusual in that it is both conformal and authalic.
  886 That is possible because it properly embeds the sphere in 3-space, as a 
  887 notional globe.
  888 
  889 This is handy as an intermediate step in lots of transforms, as 
  890 Cartesian 3-space is cleaner to work with than spherical 2-space.
  891 
  892 Higher dimensional indices are preserved, so that "rider" indices (such as 
  893 pen value) are propagated.
  894 
  895 There is no oblique transform for t_unit_sphere, largely because 
  896 it's so easy to rotate the output using t_linear once it's out into 
  897 Cartesian space.  In fact, the other projections implement oblique
  898 transforms by
  899 L<wrapping|PDL::Transform/t_wrap>
  900 L<t_linear|PDL::Transform/t_linear> with
  901 L</t_unit_sphere>.
  902 
  903 OPTIONS:
  904 
  905 =over 3
  906 
  907 =item radius, Radius (default 1.0) 
  908 
  909 The radius of the sphere, for the inverse transform.  (Radius is ignored
  910 in the forward transform).  Defaults to 1.0 so that the resulting Cartesian
  911 coordinates are in units of "body radii".
  912 
  913 =back
  914 
  915 =cut
  916 
  917 sub t_unit_sphere {
  918   my($me) = _new(@_,'Unit Sphere Projection'); 
  919   $me->{odim} = 3;
  920 
  921   $me->{params}->{otype} = ['X','Y','Z'];
  922   $me->{params}->{ounit} = ['body radii','body radii','body radii'];
  923 
  924   $me->{params}->{r} = topdl(_opt($me->{options},
  925                 ['r','radius','Radius'],
  926                 1.0)
  927                )->at(0);
  928                     
  929   
  930   $me->{func} = sub {
  931     my($d,$o) = @_;
  932     my(@dims) = $d->dims;
  933     $dims[0] ++;
  934     my $out = zeroes(@dims);
  935     
  936     my($thetaphi) = ((defined $o->{conv} && $o->{conv} != 1.0) ? 
  937              $d->slice("0:1") * $o->{conv} : $d->slice("0:1")
  938              );
  939 
  940     my $th = $thetaphi->slice("(0)");
  941     my $ph = $thetaphi->slice("(1)");
  942 
  943     # use x as a holding tank for the cos-phi multiplier
  944     $out->slice("(0)") .= $o->{r} * cos($ph) ;
  945     $out->slice("(1)") .= $out->slice("(0)") * sin($th);
  946     $out->slice("(0)") *= cos($th);
  947 
  948     $out->slice("(2)") .= $o->{r} * sin($ph);
  949 
  950     if($d->dim(0) > 2) {
  951     $out->slice("3:-1") .= $d->slice("2:-1");
  952     }
  953 
  954     $out;
  955 
  956   };
  957 
  958   $me->{inv} = sub {
  959     my($d,$o) = @_;
  960     
  961     my($d0,$d1,$d2) = ($d->slice("(0)"),$d->slice("(1)"),$d->slice("(2)"));
  962     my($r) = sqrt(($d->slice("0:2")*$d->slice("0:2"))->sumover);
  963     my(@dims) = $d->dims;
  964     $dims[0]--;
  965     my($out) = zeroes(@dims);
  966     
  967     $out->slice("(0)") .= atan2($d1,$d0);
  968     $out->slice("(1)") .= asin($d2/$r);
  969 
  970     if($d->dim(0) > 3) {
  971     $out->slice("2:-1") .= $d->slice("3:-1");
  972     }
  973 
  974     $out->slice("0:1") /= $o->{conv}
  975       if(defined $o->{conv} && $o->{conv} != 1.0);
  976 
  977     $out;
  978   };
  979 
  980     
  981   $me;
  982 }
  983 
  984 ######################################################################
  985 
  986 =head2 t_rot_sphere
  987 
  988 =for usage
  989 
  990     $t = t_rot_sphere({origin=>[<theta>,<phi>],roll=>[<roll>]});
  991 
  992 =for ref
  993 
  994 (Cartography) Generate oblique projections
  995 
  996 You feed in the origin in (theta,phi) and a roll angle, and you get back 
  997 out (theta', phi') coordinates.  This is useful for making oblique or 
  998 transverse projections:  just compose t_rot_sphere with your favorite
  999 projection and you get an oblique one.
 1000 
 1001 Most of the projections automagically compose themselves with t_rot_sphere
 1002 if you feed in an origin or roll angle.
 1003 
 1004 t_rot_sphere converts the base plate carree projection (straight lon, straight
 1005 lat) to a Cassini projection.
 1006 
 1007 OPTIONS
 1008 
 1009 =over 3
 1010 
 1011 =item STANDARD POSITIONAL OPTIONS
 1012 
 1013 =back
 1014 
 1015 =cut
 1016 
 1017 # helper routine for making the rotation matrix
 1018 sub _rotmat {
 1019   my($th,$ph,$r) = @_;
 1020   
 1021   pdl( [ cos($th) ,  -sin($th),    0  ],   # apply theta
 1022        [ sin($th) ,   cos($th),    0  ],
 1023        [  0,          0,           1  ] )
 1024     x
 1025   pdl( [ cos($ph),    0,  -sin($ph)  ], # apply phi
 1026      [ 0,           1,    0       ],
 1027      [ sin($ph),   0,  cos($ph)  ] )
 1028     x
 1029   pdl( [ 1,         0 ,       0      ], # apply roll last
 1030      [ 0,    cos($r),   -sin($r)   ], 
 1031      [ 0,    sin($r),    cos($r)   ])
 1032     ;
 1033 }
 1034 
 1035 sub t_rot_sphere {
 1036     my($me) = _new(@_,'Spherical rotation');
 1037 
 1038 
 1039     my($th,$ph) = $me->{params}->{o}->list;
 1040     my($r) = $me->{params}->{roll}->at(0);
 1041 
 1042     my($rotmat) = _rotmat($th,$ph,$r);
 1043 
 1044     my $out =  t_wrap( t_linear(m=>$rotmat, d=>3), t_unit_sphere());
 1045     $out->{itype} = $me->{itype};
 1046     $out->{iunit} = $me->{iunit};
 1047     $out->{otype} = ['rotated longitude','rotated latitude'];
 1048     $out->{ounit} = $me->{iunit};
 1049 
 1050     $out;
 1051 }
 1052 
 1053 
 1054 ######################################################################
 1055 
 1056 =head2 t_orthographic
 1057 
 1058 =for usage
 1059 
 1060     $t = t_orthographic(<options>);
 1061 
 1062 =for ref
 1063 
 1064 (Cartography) Ortho. projection (azimuthal; perspective)
 1065 
 1066 This is a perspective view as seen from infinite distance.  You
 1067 can specify the sub-viewer point in (lon,lat) coordinates, and a rotation
 1068 angle of the map CW from (north=up).  This is equivalent to specify
 1069 viewer roll angle CCW from (north=up).
 1070 
 1071 t_orthographic is a convenience interface to t_unit_sphere -- it is implemented
 1072 as a composition of a t_unit_sphere call, a rotation, and a slice.
 1073 
 1074 [*] In the default case where the near hemisphere is mapped, the
 1075 inverse exists.  There is no single inverse for the whole-sphere case,
 1076 so the inverse transform superimposes everything on a single
 1077 hemisphere.  If you want an invertible 3-D transform, you want
 1078 L</t_unit_sphere>.
 1079 
 1080 OPTIONS
 1081 
 1082 =over 3
 1083 
 1084 =item STANDARD POSITIONAL OPTIONS
 1085 
 1086 =item m, mask, Mask, h, hemisphere, Hemisphere [default 'near']
 1087 
 1088 The hemisphere to keep in the projection (see L<PDL::Transform::Cartography>).
 1089 
 1090 =back 
 1091 
 1092 NOTES
 1093 
 1094 Alone of the various projections, this one does not use
 1095 L</t_rot_sphere> to handle the standard options, because
 1096 the cartesian coordinates of the rotated sphere are already correctly
 1097 projected -- t_rot_sphere would put them back into (theta', phi')
 1098 coordinates.
 1099 
 1100 =cut
 1101 
 1102 sub t_orthographic {
 1103     my($me) = _new(@_,'Orthographic Projection');
 1104     
 1105     $me->{otype} = ['projected X','projected Y'];
 1106     $me->{ounit} = ['body radii','body radii'];
 1107 
 1108     my $m= _opt($me->{options},
 1109         ['m','mask','Mask','h','hemi','hemisphere','Hemisphere'],
 1110         1);
 1111     if($m=~m/^b/i) {
 1112     $me->{params}->{m} = 0;
 1113     } elsif($m=~m/^n/i) {
 1114     $me->{params}->{m} = 1;
 1115     } elsif($m=~m/^f/i) {
 1116     $me->{params}->{m} = 2;
 1117     } else {
 1118     $me->{params}->{m} = $m;
 1119     }
 1120 
 1121     my $origin= $me->{params}->{o} * $RAD2DEG;
 1122     my $roll = $me->{params}->{roll} * $RAD2DEG;
 1123 
 1124     $me->{params}->{t_int} = 
 1125     t_compose(
 1126           t_linear(rot=>[90 - $origin->at(1),
 1127                  0,
 1128                  90 + $origin->at(0)],
 1129                d=>3),
 1130           t_unit_sphere(u=>$me->{params}->{u})
 1131           );
 1132     
 1133     $me->{params}->{t_int} = 
 1134     t_compose(
 1135           t_linear(rot=>[0,0,$roll->at(0)],d=>3),
 1136           $me->{params}->{t_int}
 1137           )
 1138         if($roll->at(0));
 1139 
 1140     $me->{name} = "orthographic";
 1141 
 1142     $me->{idim} = 2;
 1143     $me->{odim} = 2;
 1144 
 1145     $me->{func} = sub {
 1146     my ($d,$o) = @_ ;
 1147     my ($out) = $o->{t_int}->apply($d);
 1148     if($o->{m}) {
 1149         my $idx;
 1150         $idx = whichND($out->slice("(2)") < 0)
 1151         if($o->{m} == 1);
 1152         $idx = whichND($out->slice("(2)") > 0)
 1153         if($o->{m} == 2);
 1154         if(defined $idx && ref $idx eq 'PDL' && $idx->nelem){
 1155           $out->slice("(0)")->range($idx) .= $o->{bad};
 1156           $out->slice("(1)")->range($idx) .= $o->{bad};
 1157         }
 1158     }
 1159 
 1160     my($d0) = $out->dim(0);
 1161 
 1162     # Remove the Z direction
 1163     ($d0 > 3) ? $out->slice(pdl(0,1,3..$d0-1)) : $out->slice("0:1");
 1164 
 1165     };
 1166 
 1167     # This is slow to run, quick to code -- could be made better by
 1168     # having its own 2-d inverse instead of calling the internal one.
 1169     $me->{inv} = sub {
 1170     my($d,$o) = @_;
 1171     my($d1) = $d->slice("0:1");
 1172     my(@dims) = $d->dims;
 1173     $dims[0]++;
 1174     my($out) = zeroes(@dims);
 1175     $out->slice("0:1") .= $d1;
 1176     $out->slice("3:-1") .= $d->slice("2:-1")
 1177         if($dims[0] > 3);
 1178 
 1179     $out->slice("(2)") .= sqrt(1 - ($d1*$d1)->sumover);
 1180     $out->slice("(2)") *= -1 if($o->{m} == 2);
 1181 
 1182     $o->{t_int}->invert($out);
 1183     };
 1184 
 1185     $me;
 1186 }
 1187 
 1188 ######################################################################
 1189 
 1190 =head2 t_carree
 1191 
 1192 =for usage
 1193 
 1194     $t = t_carree(<options>);
 1195 
 1196 =for ref
 1197 
 1198 (Cartography) Plate Carree projection (cylindrical; equidistant)
 1199 
 1200 This is the simple Plate Carree projection -- also called a "lat/lon plot".
 1201 The horizontal axis is theta; the vertical axis is phi.  This is a no-op
 1202 if the angular unit is radians; it is a simple scale otherwise. 
 1203 
 1204 OPTIONS
 1205 
 1206 =over 3
 1207 
 1208 =item STANDARD POSITIONAL OPTIONS
 1209 
 1210 =item s, std, standard, Standard (default 0)
 1211 
 1212 The standard parallel where the transformation is conformal.  Conformality
 1213 is achieved by shrinking of the horizontal scale to match the 
 1214 vertical scale (which is correct everywhere).
 1215 
 1216 =back
 1217 
 1218 =cut
 1219 
 1220 @PDL::Transform::Cartography::Carree::ISA = ('PDL::Transform::Cartography');
 1221 
 1222 *t_caree = *t_carree; # compat alias for poor spellers
 1223 sub t_carree {
 1224     my($me) = _new(@_,'Plate Carree Projection');
 1225     my $p = $me->{params};
 1226 
 1227     $me->{otype} = ['projected longitude','latitude'];
 1228     $me->{ounit} = ['proj. body radii','body radii'];
 1229 
 1230     $p->{stretch} = cos($p->{std});
 1231 
 1232     $me->{func} = sub {
 1233     my($d,$o) = @_;
 1234     my($out) = $d->is_inplace ? $d : $d->copy;
 1235     $out->slice("0:1") *= $o->{conv};
 1236     $out->slice("0") *= $p->{stretch};
 1237     $out;
 1238     };
 1239     
 1240     $me->{inv} = sub {
 1241     my($d,$o)= @_;
 1242     my($out) = $d->is_inplace ? $d : $d->copy;
 1243     $out->slice("0:1") /= $o->{conv};
 1244     $out->slice("0") /= $p->{stretch};
 1245     $out;
 1246     };
 1247     
 1248     $me->_finish;
 1249 }
 1250 
 1251 ######################################################################
 1252 
 1253 =head2 t_mercator
 1254 
 1255 =for usage
 1256 
 1257     $t = t_mercator(<options>);
 1258 
 1259 =for ref
 1260 
 1261 (Cartography) Mercator projection (cylindrical; conformal)
 1262 
 1263 This is perhaps the most famous of all map projections: meridians are mapped
 1264 to parallel vertical lines and parallels are unevenly spaced horizontal lines.
 1265 The poles are shifted to +/- infinity.  The output values are in units of 
 1266 globe-radii for easy conversion to kilometers; hence the horizontal extent
 1267 is -pi to pi.
 1268 
 1269 You can get oblique Mercator projections by specifying the C<origin> or
 1270 C<roll> options; this is implemented via L</t_rot_sphere>.
 1271 
 1272 OPTIONS
 1273 
 1274 =over 3
 1275 
 1276 =item STANDARD POSITIONAL OPTIONS
 1277 
 1278 =item c, clip, Clip (default 75 [degrees])
 1279 
 1280 The north/south clipping boundary of the transformation.  Because the poles are
 1281 displaced to infinity, many applications require a clipping boundary.  The
 1282 value is in whatever angular unit you set with the standard 'units' option.
 1283 The default roughly matches interesting landforms on Earth.
 1284 For no clipping at all, set b=>0.  For asymmetric clipping, use a 2-element
 1285 list ref or ndarray.
 1286 
 1287 =item s, std, Standard (default 0)
 1288 
 1289 This is the parallel at which the map has correct scale.  The scale
 1290 is also correct at the parallel of opposite sign.  
 1291 
 1292 =back
 1293 
 1294 =cut
 1295 
 1296 
 1297 @PDL::Transform::Cartography::Mercator::ISA = ('PDL::Transform::Cartography');
 1298 
 1299 sub t_mercator {
 1300     my($me) = _new(@_,'Mercator Projection');
 1301     my $p = $me->{params};
 1302 
 1303 # This is a lot of shenanigans just to get the clip parallels, but what the
 1304 # heck -- it's not a hot spot and it saves copying the input data (for 
 1305 # nondestructive clipping).
 1306     $p->{c} = _opt($me->{options},
 1307            ['c','clip','Clip'],
 1308            undef);
 1309     if(defined($p->{c})) {
 1310     $p->{c} = topdl($p->{c});
 1311     $p->{c} *= $p->{conv};
 1312     } else {
 1313     $p->{c} = pdl($DEG2RAD * 75);
 1314     }
 1315     $p->{c} = abs($p->{c}) * pdl(-1,1) if($p->{c}->nelem == 1);
 1316 
 1317     $p->{c} = log(tan(($p->{c}/2) + $PI/4));       
 1318     $p->{c} = [$p->{c}->list];
 1319 
 1320     $p->{std} = topdl(_opt($me->{options},
 1321              ['s','std','standard','Standard'],
 1322              0))->at(0) * $p->{conv};
 1323 
 1324     if($p->{std} == 0) {
 1325       $me->{otype} = ['longitude','tan latitude'];
 1326       $me->{ounit} = ['radians',' '] unless(defined $me->{ounit});
 1327     } else {
 1328       $me->{otype} = ['proj. longitude','proj. tan latitude'];
 1329       $me->{ounit} = ['radians',' '] unless(defined $me->{ounit});
 1330     }
 1331 
 1332     $p->{stretch} = cos($p->{std});
 1333 
 1334     $me->{func} = sub {
 1335     my($d,$o) = @_;
 1336 
 1337     my($out) = $d->is_inplace ? $d : $d->copy;
 1338 
 1339     $out->slice("0:1") *= $o->{conv};
 1340 
 1341     $out->slice("(1)") .= log(tan($out->slice("(1)")/2 + $PI/4));
 1342     $out->slice("(1)") .= $out->slice("(1)")->clip(@{$o->{c}})
 1343         unless($o->{c}->[0] == $o->{c}->[1]);
 1344 
 1345     $out->slice("0:1") *= $o->{stretch};
 1346     $out->slice("0:1") /= $o->{oconv} if(defined $o->{oconv});
 1347                 
 1348     $out;
 1349     };
 1350 
 1351     $me->{inv} = sub {
 1352     my($d,$o) = @_;
 1353     my($out) = $d->is_inplace? $d : $d->copy;
 1354 
 1355     $out->slice("0:1") *= $o->{oconv} if defined($o->{oconv});
 1356     $out->slice("0:1") /= $o->{stretch};
 1357     $out->slice("(1)") .= (atan(exp($out->slice("(1)"))) - $PI/4)*2;
 1358     $out->slice("0:1") /= $o->{conv};
 1359 
 1360     $out;
 1361     };
 1362 
 1363     $me->_finish;
 1364 }    
 1365 
 1366 ######################################################################
 1367 
 1368 =head2 t_utm
 1369 
 1370 =for usage
 1371 
 1372   $t = t_utm(<zone>,<options>);
 1373 
 1374 =for ref
 1375 
 1376 (Cartography) Universal Transverse Mercator projection (cylindrical)
 1377 
 1378 This is the internationally used UTM projection, with 2 subzones 
 1379 (North/South).  The UTM zones are parametrized individually, so if you
 1380 want a Zone 30 map you should use C<t_utm(30)>.  By default you get
 1381 the northern subzone, so that locations in the southern hemisphere get 
 1382 negative Y coordinates.  If you select the southern subzone (with the 
 1383 "subzone=>-1" option), you get offset southern UTM coordinates.  
 1384 
 1385 The 20-subzone military system is not yet supported.  If/when it is
 1386 implemented, you will be able to enter "subzone=>[a-t]" to select a N/S
 1387 subzone.
 1388 
 1389 Note that UTM is really a family of transverse Mercator projections
 1390 with different central meridia.  Each zone properly extends for six
 1391 degrees of longitude on either side of its appropriate central meridian,
 1392 with Zone 1 being centered at -177 degrees longitude (177 west).
 1393 Properly speaking, the zones only extend from 80 degrees south to 84 degrees
 1394 north; but this implementation lets you go all the way to 90 degrees.
 1395 The default UTM coordinates are meters.  The origin for each zone is
 1396 on the equator, at an easting of -500,000 meters.
 1397 
 1398 The default output units are meters, assuming that you are wanting a
 1399 map of the Earth.  This will break for bodies other than Earth (which
 1400 have different radii and hence different conversions between lat/lon
 1401 angle and meters).
 1402 
 1403 The standard UTM projection has a slight reduction in scale at the
 1404 prime meridian of each zone: the transverse Mercator projection's
 1405 standard "parallels" are 180km e/w of the central meridian.  However,
 1406 many Europeans prefer the "Gauss-Kruger" system, which is virtually
 1407 identical to UTM but with a normal tangent Mercator (standard parallel
 1408 on the prime meridian).  To get this behavior, set "gk=>1".
 1409 
 1410 Like the rest of the PDL::Transform::Cartography package, t_utm uses a
 1411 spherical datum rather than the "official" ellipsoidal datums for the
 1412 UTM system.
 1413 
 1414 This implementation was derived from the rather nice description by 
 1415 Denis J. Dean, located on the web at:
 1416 http://www.cnr.colostate.edu/class_info/nr502/lg3/datums_coordinates/utm.html
 1417 
 1418 OPTIONS
 1419 
 1420 =over 3
 1421 
 1422 =item STANDARD OPTIONS 
 1423 
 1424 (No positional options -- Origin and Roll are ignored) 
 1425 
 1426 =item ou, ounit, OutputUnit (default 'meters')
 1427 
 1428 (This is likely to become a standard option in a future release) The
 1429 unit of the output map.  By default, this is 'meters' for UTM, but you
 1430 may specify 'deg' or 'km' or even (heaven help us) 'miles' if you
 1431 prefer.
 1432 
 1433 =item sz, subzone, SubZone (default 1)
 1434 
 1435 Set this to -1 for the southern hemisphere subzone.  Ultimately you
 1436 should be able to set it to a letter to get the corresponding military
 1437 subzone, but that's too much effort for now.
 1438 
 1439 =item gk, gausskruger (default 0)
 1440 
 1441 Set this to 1 to get the (European-style) tangent-plane Mercator with
 1442 standard parallel on the prime meridian.  The default of 0 places the
 1443 standard parallels 180km east/west of the prime meridian, yielding better 
 1444 average scale across the zone.  Setting gk=>1 makes the scale exactly 1.0
 1445 at the central meridian, and >1.0 everywhere else on the projection. 
 1446 The difference in scale is about 0.3%.
 1447 
 1448 =back
 1449 
 1450 =cut
 1451 
 1452 sub t_utm {
 1453   my $zone = (int(shift)-1) % 60 + 1;
 1454   my($x) = _new(@_,"UTM-$zone");
 1455   my $opt = $x->{options};
 1456 
 1457   ## Make sure that there is a conversion (default is 'meters')
 1458   $x->{ounit} = ['meter','meter'] unless defined($x->{ounit});
 1459   $x->{ounit} = [$x->{ounit},$x->{ounit}] unless ref($x->{ounit});
 1460   $x->{params}->{oconv} = _uconv($x->{ounit}->[0]);
 1461 
 1462   ## Define our zone and NS offset 
 1463   my $subzone = _opt($opt,['sz', 'subzone', 'SubZone'],1);
 1464   my $offset = zeroes(2);
 1465   $offset->slice("0") .= 5e5*(2*$PI/40e6)/$x->{params}->{oconv};
 1466   $offset->slice("1") .= ($subzone < 0) ? $PI/2/$x->{params}->{oconv} : 0;
 1467 
 1468   my $merid = ($zone * 6) - 183;
 1469 
 1470   my $gk = _opt($opt,['gk','gausskruger'],0);
 1471   
 1472   my($me) = t_compose(t_linear(post=>$offset,
 1473                    rot=>-90
 1474                    ),
 1475 
 1476               t_mercator(o=>[$merid,0], 
 1477                  r=>90, 
 1478                  ou=>$x->{ounit}, 
 1479                  s=>$gk ? 0 : ($RAD2DEG * (180/6371))
 1480                 )
 1481               );
 1482 
 1483 
 1484   my $s = ($zone < 0) ? "S Hemisphere " : "";
 1485   $me->{otype} = ["UTM-$zone Easting","${s}Northing"];
 1486   $me->{ounit} = $x->{ounit};
 1487 
 1488   return $me;
 1489 }
 1490 
 1491 ######################################################################
 1492 
 1493 =head2 t_sin_lat
 1494 
 1495 =for usage
 1496 
 1497     $t = t_sin_lat(<options>);
 1498 
 1499 =for ref
 1500 
 1501 (Cartography) Cyl. equal-area projection (cyl.; authalic)
 1502 
 1503 This projection is commonly used in solar Carrington plots; but not much
 1504 for terrestrial mapping.
 1505 
 1506 OPTIONS
 1507 
 1508 =over 3
 1509 
 1510 =item STANDARD POSITIONAL OPTIONS
 1511 
 1512 =item s,std, Standard (default 0)
 1513 
 1514 This is the parallel at which the map is conformal.  It is also conformal
 1515 at the parallel of opposite sign.  The conformality is achieved by matched
 1516 vertical stretching and horizontal squishing (to achieve constant area).
 1517 
 1518 =back
 1519 
 1520 =cut
 1521 
 1522 @PDL::Transform::Cartography::SinLat::ISA = ('PDL::Transform::Cartography');
 1523 sub t_sin_lat {
 1524     my($me) = _new(@_,"Sine-Latitude Projection");
 1525 
 1526     $me->{params}->{std} = topdl(_opt($me->{options},
 1527                 ['s','std','standard','Standard'],
 1528                 0))->at(0) * $me->{params}->{conv};
 1529 
 1530     if($me->{params}->{std} == 0) {
 1531       $me->{otype} = ['longitude','sin latitude'];
 1532       $me->{ounit} = ['radians',' ']; # nonzero but blank!
 1533     } else {
 1534       $me->{otype} = ['proj. longitude','proj. sin latitude'];
 1535       $me->{ounit} = ['radians',' '];
 1536     }
 1537 
 1538     $me->{params}->{stretch} = sqrt(cos($me->{params}->{std}));
 1539 
 1540     $me->{func} = sub {
 1541     my($d,$o) = @_;
 1542     my($out) = $d->is_inplace ? $d : $d->copy;
 1543 
 1544     $out->slice("0:1") *= $me->{params}->{conv};
 1545     $out->slice("(1)") .= sin($out->slice("(1)")) / $o->{stretch};
 1546     $out->slice("(0)") *= $o->{stretch};
 1547     $out;
 1548     };
 1549 
 1550     $me->{inv} = sub {
 1551     my($d,$o) = @_;
 1552     my($out) = $d->is_inplace ? $d : $d->copy;
 1553     $out->slice("(1)") .= asin($out->slice("(1)") * $o->{stretch});
 1554     $out->slice("(0)") /= $o->{stretch};
 1555     $out->slice("0:1") /= $me->{params}->{conv};
 1556     $out;
 1557     };
 1558 
 1559     $me->_finish;
 1560 }
 1561 
 1562 ######################################################################
 1563 
 1564 =head2 t_sinusoidal
 1565 
 1566 =for usage
 1567 
 1568     $t = t_sinusoidal(<options>);
 1569 
 1570 =for ref
 1571 
 1572 (Cartography) Sinusoidal projection (authalic)
 1573 
 1574 Sinusoidal projection preserves the latitude scale but scales
 1575 longitude according to sin(lat); in this respect it is the companion to
 1576 L</t_sin_lat>, which is also authalic but preserves the
 1577 longitude scale instead.  
 1578 
 1579 OPTIONS
 1580 
 1581 =over 3
 1582 
 1583 =item STANDARD POSITIONAL OPTIONS
 1584 
 1585 =back
 1586 
 1587 =cut
 1588 
 1589 sub t_sinusoidal {
 1590   my($me) = _new(@_,"Sinusoidal Projection");
 1591   $me->{otype} = ['longitude','latitude'];
 1592   $me->{ounit} = [' ','radians'];
 1593   
 1594   $me->{func} = sub {
 1595     my($d,$o) = @_;
 1596     my($out) = $d->is_inplace ? $d : $d->copy;
 1597     $out->slice("0:1") *= $o->{conv};
 1598 
 1599     $out->slice("(0)") *= cos($out->slice("(1)"));
 1600     $out;
 1601   };
 1602 
 1603   $me->{inv} = sub {
 1604     my($d,$o) = @_;
 1605     my($out) = $d->is_inplace ? $d : $d->copy;
 1606     my($x) = $out->slice("(0)");
 1607     my($y) = $out->slice("(1)");
 1608     
 1609     $x /= cos($out->slice("(1)"));
 1610 
 1611     my($rej) = ( (abs($x)>$PI) | (abs($y)>($PI/2)) )->flat;
 1612     $x->flat->slice($rej) .= $o->{bad};
 1613     $y->flat->slice($rej) .= $o->{bad};
 1614     
 1615     $out->slice("0:1") /= $o->{conv};
 1616     $out;
 1617   };
 1618 
 1619   $me->_finish;
 1620 }
 1621    
 1622 
 1623 
 1624 
 1625 
 1626 ######################################################################
 1627 #
 1628 # Conic projections are subclassed for easier stringification and
 1629 # parsing of the standard parallels.  The constructor gets copied
 1630 # into the current package for ease of hackage.
 1631 #
 1632 # This is a little kludgy -- it's intended for direct calling
 1633 # rather than method calling, and it puts its own class name on the
 1634 # front of the argument list.  But, hey, it works...
 1635 #
 1636 @PDL::Transform::Cartography::Conic::ISA = ('PDL::Transform::Cartography');
 1637 sub _c_new {
 1638     my($def_std) = pop;
 1639     my($me) = new('PDL::Transform::Cartography::Conic',@_); 
 1640 
 1641     my($p) = $me->{params};
 1642     $p->{std} = _opt($me->{options},['s','std','standard','Standard'],
 1643              $def_std);
 1644     $p->{std} = topdl($p->{std}) * $me->{params}->{conv};
 1645     $p->{std} = topdl([$PI/2 * ($p->{std}<0 ? -1 : 1), $p->{std}->at(0)])
 1646     if($p->{std}->nelem == 1);
 1647     
 1648     $me->{params}->{cylindrical} = 1
 1649     if(approx($p->{std}->slice("0"),-$p->{std}->slice("1")));
 1650 
 1651     $me;
 1652 }
 1653 
 1654 sub PDL::Transform::Cartography::Conic::stringify {
 1655     my($me) = shift;
 1656     my($out) = $me->SUPER::stringify;
 1657 
 1658     $out .= sprintf("\tStd parallels: %6.2f,%6.2f %s\n",
 1659             $me->{params}->{std}->at(0) / $me->{params}->{conv}, 
 1660             $me->{params}->{std}->at(1) / $me->{params}->{conv}, 
 1661             $me->{params}->{u});
 1662     $out;
 1663 }
 1664 
 1665 ######################################################################
 1666 
 1667 =head2 t_conic
 1668 
 1669 =for usage
 1670 
 1671     $t = t_conic(<options>)
 1672 
 1673 =for ref
 1674 
 1675 (Cartography) Simple conic projection (conic; equidistant)
 1676 
 1677 This is the simplest conic projection, with parallels mapped to
 1678 equidistant concentric circles.  It is neither authalic nor conformal.
 1679 This transformation is also referred to as the "Modified Transverse
 1680 Mercator" projection in several maps of Alaska published by the USGS;
 1681 and the American State of New Mexico re-invented the projection in
 1682 1936 for an official map of that State.
 1683 
 1684 OPTIONS
 1685 
 1686 =over 3
 1687 
 1688 =item STANDARD POSITIONAL OPTIONS
 1689 
 1690 =item s, std, Standard (default 29.5, 45.5)
 1691 
 1692 The locations of the standard parallel(s) (where the cone intersects
 1693 the surface of the sphere).  If you specify only one then the other is
 1694 taken to be the nearest pole.  If you specify both of them to be one
 1695 pole then you get an equidistant azimuthal map.  If you specify both
 1696 of them to be opposite and equidistant from the equator you get a
 1697 Plate Carree projection.
 1698 
 1699 =back
 1700 
 1701 =cut
 1702 
 1703 sub t_conic {
 1704     my($me) = _c_new(@_,"Simple Conic Projection",[29.5,45.5]);
 1705 
 1706     my($p) = $me->{params};
 1707 
 1708     if($p->{cylindrical}) {
 1709     print STDERR "Simple conic: degenerate case; using Plate Carree\n"
 1710         if($PDL::verbose);
 1711     return t_carree($me->{options});
 1712     }
 1713 
 1714     $p->{n} = ((cos($p->{std}->slice("(0)")) - cos($p->{std}->slice("(1)")))
 1715            /
 1716            ($p->{std}->slice("(1)") - $p->{std}->slice("(0)")));
 1717 
 1718     $p->{G} = cos($p->{std}->slice("(0)"))/$p->{n} + $p->{std}->slice("(0)");
 1719 
 1720     $me->{otype} = ['Conic X','Conic Y'];
 1721     $me->{ounit} = ['Proj. radians','Proj. radians'];
 1722 
 1723     $me->{func} = sub {
 1724     my($d,$o) = @_;
 1725     my($out) = $d->is_inplace ? $d : $d->copy;
 1726 
 1727     my($rho) = $o->{G} - $d->slice("(1)") * $o->{conv};
 1728     my($theta) = $o->{n} * $d->slice("(0)") * $o->{conv};
 1729     
 1730     $out->slice("(0)") .= $rho * sin($theta);
 1731     $out->slice("(1)") .= $o->{G} - $rho * cos($theta);
 1732 
 1733     $out;
 1734     };
 1735 
 1736     $me->{inv} = sub {
 1737     my($d,$o) = @_;
 1738     my($out) = $d->is_inplace ? $d : $d->copy;
 1739 
 1740     my($x) = $d->slice("(0)");
 1741     my($y) = $o->{G} - $d->slice("(1)");
 1742     my($rho) = sqrt($x*$x + $y*$y);
 1743     $rho *= -1 if($o->{n}<0);
 1744 
 1745     my($theta) = ($o->{n} < 0) ? atan2(-$x,-$y) : atan2($x,$y);
 1746     
 1747     $out->slice("(1)") .= $o->{G} - $rho;
 1748     $out->slice("(1)")->where(($out->slice("(1)") < -$PI/2) | ($out->slice("(1)") > $PI/2))
 1749         .= $o->{bad};
 1750 
 1751     $out->slice("(0)") .= $theta / $o->{n};
 1752     $out->slice("(0)")->where(($out->slice("(0)") < -$PI) | ($out->slice("(0)") > $PI/2))
 1753         .= $o->{bad};
 1754 
 1755     $out->slice("0:1") /= $o->{conv};
 1756 
 1757     $out;
 1758     };
 1759 
 1760     $me->_finish;
 1761 }
 1762 
 1763 
 1764 
 1765 
 1766 ######################################################################
 1767 
 1768 =head2 t_albers
 1769 
 1770 =for usage
 1771 
 1772     $t = t_albers(<options>)
 1773 
 1774 =for ref
 1775 
 1776 (Cartography) Albers conic projection (conic; authalic)
 1777 
 1778 This is the standard projection used by the US Geological Survey for
 1779 sectionals of the 50 contiguous United States of America.  
 1780 
 1781 The projection reduces to the Lambert equal-area conic (infrequently
 1782 used and not to be confused with the Lambert conformal conic,
 1783 L</t_lambert>!)  if the pole is used as one of the two
 1784 standard parallels.
 1785 
 1786 Notionally, this is a conic projection onto a cone that intersects the
 1787 sphere at the two standard parallels; it works best when the two parallels
 1788 straddle the region of interest.
 1789 
 1790 OPTIONS
 1791 
 1792 =over 3
 1793 
 1794 =item STANDARD POSITIONAL OPTIONS
 1795 
 1796 =item s, std, standard, Standard (default (29.5,45.5))
 1797 
 1798 The locations of the standard parallel(s).  If you specify only one then 
 1799 the other is taken to be the nearest pole and a Lambert Equal-Area Conic
 1800 map results.  If you specify both standard parallels to be the same pole,
 1801 then the projection reduces to the Lambert Azimuthal Equal-Area map as
 1802 aq special case.  (Note that L</t_lambert> is Lambert's
 1803 Conformal Conic, the most commonly used of Lambert's projections.)
 1804 
 1805 The default values for the standard parallels are those chosen by Adams
 1806 for maps of the lower 48 US states: (29.5,45.5).  The USGS recommends
 1807 (55,65) for maps of Alaska and (8,18) for maps of Hawaii -- these latter
 1808 are chosen to also include the Canal Zone and Philippine Islands farther
 1809 south, which is why both of those parallels are south of the Hawaiian islands.
 1810 
 1811 The transformation reduces to the cylindrical equal-area (sin-lat)
 1812 transformation in the case where the standard parallels are opposite and
 1813 equidistant from the equator, and in fact this is implemented by a call to
 1814 t_sin_lat.
 1815 
 1816 =back
 1817 
 1818 =cut
 1819 
 1820 sub t_albers  {
 1821     my($me) = _c_new(@_,"Albers Equal-Area Conic Projection",[29.5,45.5]);
 1822 
 1823     my($p) = $me->{params};
 1824 
 1825     if($p->{cylindrical}) {
 1826     print STDERR "Albers equal-area conic: degenerate case; using equal-area cylindrical\n"
 1827         if($PDL::verbose);
 1828     return t_sin_lat($me->{options});
 1829     }
 1830 
 1831     $p->{n} = sin($p->{std})->sumover / 2;
 1832     $p->{C} = (cos($p->{std}->slice("(1)"))*cos($p->{std}->slice("(1)")) +
 1833              2 * $p->{n} * sin($p->{std}->slice("(1)")) );
 1834     $p->{rho0} = sqrt($p->{C}) / $p->{n}; 
 1835 
 1836     $me->{otype} = ['Conic X','Conic Y'];
 1837     $me->{ounit} = ['Proj. radians','Proj. radians'];
 1838 
 1839     $me->{func} = sub {
 1840     my($d,$o) = @_;
 1841     my($out) = $d->is_inplace ? $d : $d->copy;
 1842 
 1843     my($rho) = sqrt( $o->{C} - 2 * $o->{n} * sin($d->slice("(1)") * $o->{conv}) ) / $o->{n};
 1844     my($theta) = $o->{n} * $d->slice("(0)") * $o->{conv};
 1845 
 1846     $out->slice("(0)") .= $rho * sin($theta);
 1847     $out->slice("(1)") .= $p->{rho0} - $rho * cos($theta);
 1848     $out;
 1849     };
 1850 
 1851     $me->{inv} = sub {
 1852     my($d,$o) = @_;
 1853 
 1854     my($out) = $d->is_inplace ? $d : $d->copy;
 1855 
 1856     my($x) = $d->slice("(0)");
 1857     my($y) = $o->{rho0} - $d->slice("(1)");
 1858 
 1859     my($theta) = ($o->{n} < 0) ? atan2 -$x,-$y : atan2 $x, $y;
 1860     my($rho) = sqrt( $x*$x + $y*$y ) * $o->{n};
 1861 
 1862     $out->slice("(1)") .= asin( ( $o->{C} - ( $rho * $rho ) ) / (2 * $o->{n}) );
 1863 
 1864     $out->slice("(0)") .= $theta / $o->{n};
 1865     $out->slice("(0)")->where(($out->slice("(0)")>$PI) | ($out->slice("(0)")<-$PI)) .= $o->{bad};
 1866 
 1867     $out->slice("0:1") /= $o->{conv};
 1868 
 1869     $out;
 1870     };
 1871 
 1872     $me->_finish;
 1873 }
 1874 
 1875 ######################################################################
 1876 
 1877 =head2 t_lambert
 1878 
 1879 =for usage
 1880 
 1881     $t = t_lambert(<options>);
 1882 
 1883 =for ref
 1884 
 1885 (Cartography) Lambert conic projection (conic; conformal)
 1886 
 1887 Lambert conformal conic projection is widely used in aeronautical
 1888 charts and state base maps published by the USA's FAA and USGS.  It's
 1889 especially useful for mid-latitude charts.  In particular, straight lines
 1890 approximate (but are not exactly) great circle routes of up to ~2 radians.
 1891 
 1892 The default standard parallels are 33 and 45 to match the USGS state
 1893 1:500,000 base maps of the United States.  At scales of 1:500,000 and
 1894 larger, discrepancies between the spherical and ellipsoidal projections
 1895 become important; use care with this projection on spheres.
 1896 
 1897 OPTIONS
 1898 
 1899 =over 3
 1900 
 1901 =item STANDARD POSITIONAL OPTIONS
 1902 
 1903 =item s, std, standard, Standard (default (33,45))
 1904 
 1905 The locations of the standard parallel(s) for the conic projection.
 1906 The transform reduces to the Mercator projection in the case where the
 1907 standard parallels are opposite and equidistant from the equator, and
 1908 in fact this is implemented by a call to t_mercator.
 1909 
 1910 =item c, clip, Clip (default [-75,75])
 1911 
 1912 Because the transform is conformal, the distant pole is displaced to
 1913 infinity.  Many applications require a clipping boundary.  The value
 1914 is in whatever angular unit you set with the standard 'unit' option.
 1915 For consistency with L</t_mercator>, clipping works the same
 1916 way even though in most cases only one pole needs it.  Set this to 0
 1917 for no clipping at all.
 1918 
 1919 =back
 1920 
 1921 =cut
 1922 
 1923 sub t_lambert {
 1924     my($me)= _c_new(@_,"Lambert Conformal Conic Projection",[33,45]);
 1925     my($p) = $me->{params};
 1926 
 1927     if($p->{cylindrical}){
 1928     print STDERR "Lambert conformal conic: std parallels are opposite & equal; using Mercator\n" 
 1929         if($PDL::verbose);
 1930     return t_mercator($me->{options});
 1931     }
 1932     
 1933     # Find clipping parallels
 1934     $p->{c} = _opt($me->{options},['c','clip','Clip'],undef);
 1935     if(defined($p->{c})) {
 1936     $p->{c} = topdl($p->{c});
 1937     } else {
 1938     $p->{c} = topdl([-75,75]);
 1939     }
 1940     $p->{c} = abs($p->{c}) * topdl([-1,1]) if($p->{c}->nelem == 1);
 1941     $p->{c} = [$p->{c}->list];
 1942 
 1943     # Prefrobnicate
 1944     if(approx($p->{std}->slice("(0)"),$p->{std}->slice("(1)"))) {
 1945     $p->{n} = sin($p->{std}->slice("(0)"));
 1946     } else {
 1947     $p->{n} = (log(cos($p->{std}->slice("(0)"))/cos($p->{std}->slice("(1)")))
 1948            / 
 1949            log( tan( $PI/4 + $p->{std}->slice("(1)")/2 )
 1950             / 
 1951             tan( $PI/4 + $p->{std}->slice("(0)")/2 )
 1952             )
 1953            );
 1954     }
 1955 
 1956     $p->{F} = ( cos($p->{std}->slice("(0)"))
 1957         *
 1958         ( tan( $PI/4 + $p->{std}->slice("(0)")/2 ) ** $p->{n} ) / $p->{n}
 1959         );
 1960 
 1961     $p->{rho0} = $p->{F};
 1962 
 1963     $me->{otype} = ['Conic X','Conic Y'];
 1964     $me->{ounit} = ['Proj. radians','Proj. radians'];
 1965 
 1966     $me->{func} = sub {
 1967     my($d,$o) = @_;
 1968     my($out) = $d->is_inplace ? $d : $d->copy;
 1969 
 1970     my($cl) = ( ($o->{c}->[0] == $o->{c}->[1]) ? 
 1971             $d->slice("(1)")*$o->{conv} :
 1972             ($d->slice("(1)")->clip(@{$o->{c}}) * $o->{conv})
 1973             );
 1974 
 1975     my($rho) = $o->{F} / ( tan($PI/4 + ($cl)/2 ) ** $o->{n} );
 1976     my($theta) = $o->{n} * $d->slice("(0)") * $o->{conv};
 1977 
 1978     $out->slice("(0)") .= $rho * sin($theta);
 1979     $out->slice("(1)") .= $o->{rho0} - $rho * cos($theta);
 1980     $out;
 1981     };
 1982 
 1983     $me->{inv} = sub {
 1984     my($d,$o) = @_;
 1985     my($out) = $d->is_inplace ? $d : $d->copy;
 1986     
 1987     my($x) = $d->slice("(0)");
 1988     my($y) = $o->{rho0} - $d->slice("(1)");
 1989 
 1990     my($rho) = sqrt($x * $x + $y * $y);
 1991     $rho *= -1 if($o->{n} < 0);
 1992     my($theta) = ($o->{n} < 0) ? atan2(-$x,-$y):(atan2 $x,$y);
 1993 
 1994 
 1995     $out->slice("(0)") .= $theta / $o->{n};
 1996     $out->slice("(0)")->where(($out->slice("(0)") > $PI) | ($out->slice("(0)") < -$PI)) .= $o->{bad};
 1997 
 1998 
 1999     $out->slice("(1)") .= 2 * atan(($o->{F}/$rho)**(1.0/$o->{n})) - $PI/2;
 2000     $out->slice("(1)")->where(($out->slice("(1)") > $PI/2) | ($out->slice("(1)") < -$PI/2)) .= $o->{bad};
 2001 
 2002     $out->slice("0:1") /= $o->{conv};
 2003 
 2004     $out;
 2005     };
 2006 
 2007 
 2008     $me->_finish;
 2009 }
 2010 
 2011 ######################################################################
 2012 
 2013 =head2 t_stereographic
 2014 
 2015 =for usage
 2016 
 2017     $t = t_stereographic(<options>);
 2018 
 2019 =for ref
 2020 
 2021 (Cartography) Stereographic projection (az.; conf.; persp.)
 2022 
 2023 The stereographic projection is a true perspective (planar) projection
 2024 from a point on the spherical surface opposite the origin of the map.  
 2025 
 2026 OPTIONS
 2027 
 2028 =over 3
 2029 
 2030 =item STANDARD POSITIONAL OPTIONS
 2031 
 2032 =item c, clip, Clip (default 120)
 2033 
 2034 This is the angular distance from the center to the edge of the 
 2035 projected map.  The default 120 degrees gives you most of the opposite
 2036 hemisphere but avoids the hugely distorted part near the antipodes.
 2037 
 2038 =back
 2039 
 2040 =cut
 2041 
 2042 sub t_stereographic {
 2043     my($me) = _new(@_,"Stereographic Projection");
 2044     
 2045     $me->{params}->{k0} = 1.0;
 2046     $me->{params}->{c} = _opt($me->{options},
 2047                   ['c','clip','Clip'],
 2048                   120) * $me->{params}->{conv};
 2049 
 2050     $me->{otype} = ['Stereo X','Stereo Y'];
 2051     $me->{ounit} = ['Proj. body radii','Proj. radians'];
 2052 
 2053     $me->{func} = sub {
 2054     my($d,$o) = @_;
 2055     my($out) = $d->is_inplace ? $d : $d->copy;
 2056 
 2057     my($th,$ph) = ($out->slice("(0)") * $o->{conv},
 2058                $out->slice("(1)") * $o->{conv});
 2059 
 2060     my($cph) = cos($ph); # gets re-used 
 2061     my($k) = 2 * $o->{k0} / (1 + cos($th) * $cph);
 2062     $out->slice("(0)") .= $k * $cph * sin($th);
 2063     $out->slice("(1)") .= $k * sin($ph);
 2064 
 2065     my($cl0) = 2*$o->{k0} / (1 + cos($o->{c}));
 2066     $out->slice("(0)")->where($k>$cl0) .= $o->{bad};
 2067     $out->slice("(1)")->where($k>$cl0) .= $o->{bad};
 2068     $out;
 2069     };
 2070 
 2071     $me->{inv} = sub {
 2072     my($d,$o) = @_;
 2073     my($out) = $d->is_inplace ? $d : $d->copy;
 2074     
 2075     my($x) = $d->slice("(0)");
 2076     my($y) = $d->slice("(1)");
 2077     my($rho) = sqrt($x*$x + $y*$y);
 2078     my($c) = 2 * atan2($rho,2*$o->{k0});
 2079     
 2080     $out->slice("(0)") .= atan2($x * sin($c), $rho * cos($c));
 2081     $out->slice("(1)") .= asin($y * sin($c) / $rho);
 2082     
 2083     $out ->slice("0:1") /= $o->{conv};
 2084     $out;
 2085     };
 2086 
 2087     $me->_finish;
 2088 }
 2089      
 2090 ######################################################################
 2091 
 2092 =head2 t_gnomonic
 2093 
 2094 =for usage
 2095 
 2096     $t = t_gnomonic(<options>);
 2097 
 2098 =for ref
 2099 
 2100 (Cartography) Gnomonic (focal-plane) projection (az.; persp.)
 2101 
 2102 The gnomonic projection projects a hemisphere onto a tangent plane.
 2103 It is useful in cartography for the property that straight lines are
 2104 great circles; and it is useful in scientific imaging because 
 2105 it is the projection generated by a simple optical system with a flat
 2106 focal plane.
 2107 
 2108 OPTIONS
 2109 
 2110 =over 3
 2111 
 2112 =item STANDARD POSITIONAL OPTIONS
 2113 
 2114 =item c, clip, Clip (default 75)
 2115 
 2116 This is the angular distance from the center to the edge of the 
 2117 projected map.  The default 75 degrees gives you most of the 
 2118 hemisphere but avoids the hugely distorted part near the horizon.
 2119 
 2120 =back
 2121 
 2122 =cut
 2123 
 2124 sub t_gnomonic {
 2125     my($me) = _new(@_,"Gnomonic Projection");
 2126     
 2127     $me->{params}->{k0} = 1.0;  # Useful for standard parallel (TBD: add one)
 2128 
 2129     $me->{params}->{c} = topdl(_opt($me->{options},
 2130                   ['c','clip','Clip'],
 2131                   75) * $me->{params}->{conv});
 2132 
 2133     $me->{params}->{c} .= $me->{params}->{c}->clip(undef,(90-1e-6)*$me->{params}->{conv});
 2134 
 2135     $me->{otype} = ['Tangent-plane X','Tangent-plane Y'];
 2136     $me->{ounit} = ['Proj. radians','Proj. radians'];
 2137 
 2138     $me->{func} = sub {
 2139     my($d,$o) = @_;
 2140     my($out) = $d->is_inplace ? $d : $d->copy;
 2141 
 2142     my($th,$ph) = ($out->slice("(0)") * $o->{conv},
 2143                $out->slice("(1)") * $o->{conv});
 2144 
 2145     my($cph) = cos($ph); # gets re-used 
 2146 
 2147     my($k) = $o->{k0} / (cos($th) * $cph);
 2148     my($cl0) = $o->{k0} / (cos($o->{c}));
 2149 
 2150     $out->slice("(0)") .= $k * $cph * sin($th);
 2151     $out->slice("(1)") .= $k * sin($ph);
 2152 
 2153     my $idx = whichND(($k > $cl0)  | ($k < 0) | (!isfinite($k)));
 2154     if($idx->nelem) {
 2155       $out->slice("(0)")->range($idx) .= $o->{bad};
 2156       $out->slice("(1)")->range($idx) .= $o->{bad};
 2157     }
 2158 
 2159     $out;
 2160     };
 2161 
 2162     $me->{inv} = sub {
 2163     my($d,$o) = @_;
 2164     my($out) = $d->is_inplace ? $d : $d->copy;
 2165 
 2166     my($x) = $d->slice("(0)");
 2167     my($y) = $d->slice("(1)");
 2168 
 2169     my($rho) = sqrt($x*$x + $y*$y);
 2170     my($c) = atan($rho/$o->{k0});
 2171 
 2172     $out->slice("(0)") .= atan2($x * sin($c), $rho * cos($c));
 2173     $out->slice("(1)") .= asin($y * sin($c) / $rho);
 2174 
 2175     my $idx = whichND($rho==0);
 2176 
 2177     if($idx->nelem) {
 2178       $out->slice("(0)")->range($idx) .= 0;
 2179       $out->slice("(1)")->range($idx) .= 0;
 2180     }
 2181     $out->slice("0:1") /= $o->{conv};
 2182     $out;
 2183     };
 2184 
 2185     $me->_finish;
 2186 }
 2187 
 2188 ######################################################################
 2189 
 2190 =head2 t_az_eqd
 2191 
 2192 =for usage
 2193 
 2194   $t = t_az_eqd(<options>);
 2195 
 2196 =for ref 
 2197 
 2198 (Cartography) Azimuthal equidistant projection (az.; equi.)
 2199 
 2200 Basic azimuthal projection preserving length along radial lines from
 2201 the origin (meridians, in the original polar aspect).  Hence, both
 2202 azimuth and distance are correct for journeys beginning at the origin.
 2203 
 2204 Applied to the celestial sphere, this is the projection made by
 2205 fisheye lenses; it is also the projection into which C<t_vertical>
 2206 puts perspective views.
 2207 
 2208 The projected plane scale is normally taken to be planetary radii;
 2209 this is useful for cartographers but not so useful for scientific
 2210 observers.  Setting the 't=>1' option causes the output scale to shift
 2211 to camera angular coordinates (the angular unit is determined by the
 2212 standard 'Units' option; default is degrees).
 2213 
 2214 OPTIONS
 2215 
 2216 =over 3
 2217 
 2218 =item STANDARD POSITIONAL OPTIONS
 2219 
 2220 =item c, clip, Clip (default 180 degrees)
 2221 
 2222 The largest angle relative to the origin.  Default is the whole sphere.
 2223 
 2224 =back
 2225 
 2226 =cut
 2227 
 2228 sub t_az_eqd {
 2229   my($me) = _new(@_,"Equidistant Azimuthal Projection");
 2230 
 2231   $me->{params}->{c} = topdl(_opt($me->{options},
 2232                 ['c','clip','Clip'],
 2233                 180) * $me->{params}->{conv});
 2234 
 2235   $me->{otype} = ['X distance','Y distance'];
 2236   $me->{ounit} = ['radians','radians'];
 2237 
 2238   $me->{func} = sub {
 2239     my($d,$o) = @_;
 2240     my($out) = $d->is_inplace ? $d : $d->copy;
 2241     
 2242     my($ph) = $d->slice("(1)") * $o->{conv};
 2243     my($th) = $d->slice("(0)") * $o->{conv};
 2244 
 2245     my $cos_c = cos($ph) * cos($th);
 2246     my $c = acos($cos_c);
 2247     my $k = $c / sin($c);
 2248     $k->where($c==0) .= 1;
 2249     
 2250     my($x,$y) = ($out->slice("(0)"), $out->slice("(1)"));
 2251 
 2252     $x .= $k * cos($ph) * sin($th);
 2253     $y .= $k * sin($ph);
 2254 
 2255     my $idx = whichND($c > $o->{c});
 2256     if($idx->nelem) {
 2257       $x->range($idx) .= $o->{bad};
 2258       $y->range($idx) .= $o->{bad};
 2259     }
 2260 
 2261     $out;
 2262   };
 2263 
 2264   $me->{inv} = sub {
 2265     my($d,$o) = @_;
 2266     my($out) = $d->is_inplace ? $d : $d->copy;
 2267     my($x) = $d->slice("(0)");
 2268     my($y) = $d->slice("(1)");
 2269 
 2270     my $rho = sqrt(($d->slice("0:1")*$d->slice("0:1"))->sumover);
 2271     # Order is important -- ((0)) overwrites $x if is_inplace!
 2272     $out->slice("(0)") .= atan2( $x * sin($rho), $rho * cos $rho );
 2273     $out->slice("(1)") .= asin( $y * sin($rho) / $rho );
 2274 
 2275     my $idx = whichND($rho == 0);
 2276     if($idx->nelem) {
 2277       $out->slice("(0)")->range($idx) .= 0;
 2278       $out->slice("(1)")->range($idx) .= 0;
 2279     }
 2280 
 2281     $out->slice("0:1") /= $o->{conv};
 2282 
 2283     $out;
 2284   };
 2285 
 2286   $me->_finish;
 2287 }
 2288 
 2289 
 2290 ######################################################################
 2291 
 2292 =head2 t_az_eqa
 2293 
 2294 =for usage
 2295 
 2296   $t = t_az_eqa(<options>);
 2297 
 2298 =for ref
 2299 
 2300 (Cartography) Azimuthal equal-area projection (az.; auth.)
 2301 
 2302 OPTIONS
 2303 
 2304 =over 3
 2305 
 2306 =item STANDARD POSITIONAL OPTIONS
 2307 
 2308 =item c, clip, Clip (default 180 degrees)
 2309 
 2310 The largest angle relative to the origin.  Default is the whole sphere.
 2311 
 2312 =back
 2313 
 2314 =cut
 2315 
 2316 sub t_az_eqa {
 2317   my($me) = _new(@_,"Equal-Area Azimuthal Projection");
 2318   
 2319   $me->{params}->{c} = topdl(_opt($me->{options},
 2320                 ['c','clip','Clip'],
 2321                 180) * $me->{params}->{conv});
 2322 
 2323   $me->{otype} = ['Azimuthal X','Azimuthal Y'];
 2324   $me->{ounit} = ['Proj. radians','Proj. radians'];
 2325 
 2326   $me->{func} = sub {
 2327     my($d,$o) = @_;
 2328     my($out) = $d->is_inplace ? $d : $d->copy;
 2329     
 2330     my($ph) = $d->slice("(1)") * $o->{conv};
 2331     my($th) = $d->slice("(0)") * $o->{conv};
 2332                 
 2333     my($c) = acos(cos($ph) * cos($th));
 2334     my($rho) = 2 * sin($c/2);
 2335     my($k) = 1.0/cos($c/2);
 2336 
 2337     my($x,$y) = ($out->slice("(0)"),$out->slice("(1)"));
 2338     $x .= $k * cos($ph) * sin($th);
 2339     $y .= $k * sin($ph);
 2340 
 2341     my $idx = whichND($c > $o->{c});
 2342     if($idx->nelem) {
 2343       $x->range($idx) .= $o->{bad};
 2344       $y->range($idx) .= $o->{bad};
 2345     }
 2346 
 2347     $out;
 2348   };
 2349 
 2350   $me->{inv} = sub {
 2351     my($d,$o) = @_;
 2352     my($out) = $d->is_inplace ? $d : $d->copy;
 2353 
 2354     my($x,$y) = ($d->slice("(0)"),$d->slice("(1)"));
 2355     my($ph,$th) = ($out->slice("(0)"),$out->slice("(1)"));
 2356     my($rho) = sqrt($x*$x + $y*$y);
 2357     my($c) = 2 * asin($rho/2);
 2358 
 2359     $ph .= asin($d->slice("(1)") * sin($c) / $rho);
 2360     $th .= atan2($x * sin($c),$rho * cos($c));
 2361 
 2362     $ph /= $o->{conv};
 2363     $th /= $o->{conv};
 2364     
 2365     $out;
 2366   };
 2367 
 2368   $me->_finish;
 2369 }
 2370 
 2371 
 2372 ######################################################################
 2373 
 2374 =head2 t_aitoff
 2375 
 2376 C<t_aitoff> in an alias for C<t_hammer>
 2377 
 2378 =head2 t_hammer
 2379 
 2380 =for ref
 2381 
 2382 (Cartography) Hammer/Aitoff elliptical projection (az.; auth.)
 2383 
 2384 The Hammer/Aitoff projection is often used to display the Celestial
 2385 sphere.  It is mathematically related to the Lambert Azimuthal Equal-Area
 2386 projection (L</t_az_eqa>), and maps the sphere to an ellipse of unit 
 2387 eccentricity, with vertical radius sqrt(2) and horizontal radius of 
 2388 2 sqrt(2).
 2389 
 2390 OPTIONS
 2391 
 2392 =over 3
 2393 
 2394 =item STANDARD POSITIONAL OPTIONS
 2395 
 2396 =back
 2397 
 2398 =cut
 2399 
 2400 *t_aitoff = *t_aitoff = \&t_hammer;
 2401 
 2402 sub t_hammer {
 2403   my($me) = _new(@_,"Hammer/Aitoff Projection");
 2404   
 2405   $me->{otype} = ['Longitude','Latitude'];
 2406   $me->{ounit} = [' ',' '];
 2407   $me->{odim} = 2;
 2408   $me->{idim} = 2;
 2409 
 2410   $me->{func} = sub {
 2411     my($d,$o) = @_;
 2412     my($out) = $d->is_inplace ? $d : $d->copy;
 2413     $out->slice("0:1") *= $o->{conv};
 2414     my($th) = $out->slice("(0)");
 2415     my($ph) = $out->slice("(1)");
 2416     my($t) = sqrt( 2 / (1 + cos($ph) * cos($th/2)));
 2417     $th .= 2 * $t * cos($ph) * sin($th/2);
 2418     $ph .= $t * sin($ph);
 2419     $out;
 2420   }
 2421   ;
 2422 
 2423   $me->{inv} = sub {
 2424     my($d,$o) = @_;
 2425     my($out) = $d->is_inplace ? $d : $d->copy;
 2426     my($x) = $out->slice("(0)");
 2427     my($y) = $out->slice("(1)");
 2428 
 2429     my($rej) = which(($x*$x/8 + $y*$y/2)->flat > 1);
 2430     
 2431     my($zz);
 2432     my($z) = sqrt( $zz = (1 - $x*$x/16 - $y*$y/4) );
 2433     $x .= 2 * atan( ($z * $x) / (4 * $zz - 2) );
 2434     $y .= asin($y * $z);
 2435     
 2436     $out->slice("0:1") /= $o->{conv};
 2437 
 2438     $x->flat->slice($rej) .= $o->{bad};
 2439     $y->flat->slice($rej) .= $o->{bad};
 2440 
 2441     $out;
 2442   };
 2443 
 2444   $me->_finish;
 2445 }
 2446 
 2447 
 2448 ######################################################################
 2449 
 2450 =head2 t_zenithal
 2451 
 2452 Vertical projections are also called "zenithal", and C<t_zenithal> is an
 2453 alias for C<t_vertical>.
 2454 
 2455 =head2 t_vertical
 2456 
 2457 =for usage
 2458 
 2459     $t = t_vertical(<options>);
 2460 
 2461 =for ref
 2462 
 2463 (Cartography) Vertical perspective projection (az.; persp.)
 2464 
 2465 Vertical perspective projection is a generalization of L<gnomonic|/t_gnomonic>
 2466 and L<stereographic|/t_stereographic> projection, and a special case of 
 2467 L<perspective|/t_perspective> projection.  It is a projection from the 
 2468 sphere onto a tangent plane from a point at the camera location.
 2469 
 2470 OPTIONS
 2471 
 2472 =over 3
 2473 
 2474 =item STANDARD POSITIONAL OPTIONS
 2475 
 2476 =item m, mask, Mask, h, hemisphere, Hemisphere [default 'near']
 2477 
 2478 The hemisphere to keep in the projection (see L<PDL::Transform::Cartography>).
 2479 
 2480 =item r0, R0, radius, d, dist, distance [default 2.0]
 2481 
 2482 The altitude of the focal plane above the center of the sphere.  The default
 2483 places the point of view one radius above the surface.
 2484 
 2485 =item t, telescope, Telescope, cam, Camera (default '')
 2486 
 2487 If this is set, then the central scale is in telescope or camera 
 2488 angular units rather than in planetary radii.  The angular units are 
 2489 parsed as with the normal 'u' option for the lon/lat specification.
 2490 If you specify a non-string value (such as 1) then you get telescope-frame
 2491 radians, suitable for working on with other transformations.
 2492 
 2493 =item f, fish, fisheye (default '')
 2494 
 2495 If this is set then the output is in azimuthal equidistant coordinates
 2496 instead of in tangent-plane coordinates.  This is a convenience function
 2497 for '(t_az_eqd) x !(t_gnomonic) x (t_vertical)'.
 2498 
 2499 =back
 2500 
 2501 =cut
 2502 
 2503 sub t_vertical {
 2504     my($me) = _new(@_,'Vertical Perspective');
 2505     my $p = $me->{params};
 2506     
 2507     my $m= _opt($me->{options},
 2508         ['m','mask','Mask','h','hemi','hemisphere','Hemisphere'],
 2509         1);
 2510 
 2511     $me->{otype} = ['Perspective X','Perspective Y'];
 2512     $me->{ounit} = ['Body radii','Body radii'];
 2513 
 2514     if($m=~m/^b/i) {
 2515     $p->{m} = 0;
 2516     } elsif($m=~m/^n/i) {
 2517     $p->{m} = 1;
 2518     } elsif($m=~m/^f/i) {
 2519     $p->{m} = 2;
 2520     } else {
 2521     $p->{m} = $m;
 2522     }
 2523 
 2524     $p->{r0} = _opt($me->{options},
 2525             ['r0','R0','radius','Radius',
 2526              'd','dist','distance','Distance'],
 2527             2.0
 2528             );
 2529     
 2530     if($p->{r0} == 0) {
 2531       print "t_vertical: r0 = 0; using t_gnomonic instead\n"
 2532     if($PDL::verbose);
 2533       return t_gnomonic($me->{options});
 2534     }
 2535     
 2536     if($p->{r0} == 1) {
 2537       print "t_vertical: r0 = 1; using t_stereographic instead\n"
 2538     if($PDL::verbose);
 2539       return t_stereographic($me->{options});
 2540     }
 2541     
 2542     
 2543     $p->{t} = _opt($me->{options},
 2544            ['t','tele','telescope','Telescope',
 2545             'cam','camera','Camera'],
 2546            undef);
 2547     
 2548     $p->{f} = _opt($me->{options},
 2549            ['f','fish','fisheye','Fisheye'],
 2550            undef);
 2551     
 2552     $p->{t} = 'rad'
 2553       if($p->{f} && !defined($p->{t}));
 2554       
 2555     $p->{tconv} = _uconv($p->{t},1) || _uconv('rad')
 2556       if(defined $p->{t});
 2557 
 2558     $me->{func} = sub {
 2559     my($d,$o) = @_;
 2560     my($out) = $d->is_inplace ? $d : $d->copy;
 2561     my($th) = $d->slice("(0)")*$o->{conv};
 2562     my($ph) = $d->slice("(1)")*$o->{conv};
 2563 
 2564     my($cph) = cos($ph);
 2565 
 2566     my($cos_c) = $cph * cos($th);
 2567 
 2568     my($k) = (($o->{r0} - 1) / 
 2569           ($o->{r0} - $cos_c));
 2570 
 2571     # If it's a telescope perspective, figure the apparent size
 2572     # of the globe and scale accordingly.
 2573     if($o->{t}) {
 2574       my($theta) = asin(1/$o->{r0});
 2575     }
 2576     
 2577     $out->slice("0:1") /= ($o->{r0} - 1.0) * ($o->{f} ? 1.0 : $o->{tconv})
 2578       if($o->{t});
 2579 
 2580 
 2581 
 2582     $out->slice("(0)") .= $cph * sin($th);
 2583     $out->slice("(1)") .= sin($ph);
 2584 
 2585     # Handle singularity at the origin
 2586     $k->where(($out->slice("(0)") == 0) & ($out->slice("(1)") == 0)) .= 0;
 2587     $out->slice("0:1") *= $k->dummy(0,2);
 2588 
 2589     if($o->{m}) {
 2590         my $idx;
 2591         $idx = whichND($cos_c < 1.0/$o->{r0})
 2592         if($o->{m} == 1);
 2593         $idx = whichND($cos_c > 1.0/$o->{r0})
 2594         if($o->{m} == 2);
 2595 
 2596         if(defined $idx && ref $idx eq 'PDL' && $idx->nelem){
 2597           $out->slice("(0)")->range($idx) .= $o->{bad};
 2598           $out->slice("(1)")->range($idx) .= $o->{bad};
 2599         }
 2600     }
 2601 
 2602 
 2603     $out;
 2604     };
 2605 
 2606     $me->{inv} = sub {
 2607     my($d,$o) = @_;
 2608     my($out) = $d->is_inplace ? $d : $d->copy;
 2609 
 2610     # Reverse the hemisphere if the mask is set to 'far'
 2611     my($P) = ($o->{m} == 2) ? -$o->{r0} : $o->{r0};
 2612 
 2613     $out->slice("0:1") *= ($P - 1.0) * ($o->{f} ? 1.0 : $o->{tconv})
 2614             if($o->{t});
 2615 
 2616     my($rho) = sqrt(sumover($d->slice("0:1") * $d->slice("0:1")));
 2617     my($sin_c) = ( (  $P - sqrt( 1 - ($rho*$rho * ($P+1)/($P-1)) ) ) /
 2618                ( ($P-1)/$rho + $rho/($P-1) )
 2619                );
 2620 
 2621     my($cos_c) = sqrt(1 - $sin_c*$sin_c);
 2622 
 2623     # Switch c's quadrant where necessary, by inverting cos(c).
 2624     if($P<0) {
 2625       my $idx = whichND($rho > ($P-1/$P));
 2626       $cos_c->range($idx) *= -1
 2627         if($idx->nelem > 0);
 2628     }
 2629 
 2630     
 2631     $out->slice("(0)") .= atan( $d->slice("(0)") * $sin_c / ($rho * $cos_c) );
 2632     $out->slice("(1)") .= asin( $d->slice("(1)") * $sin_c / $rho );
 2633 
 2634     $out->slice("0:1") /= $o->{conv};
 2635 
 2636     $out;
 2637     };
 2638           
 2639 
 2640     # Compose on both front and back as necessary.
 2641     return t_compose( t_scale(1.0/$p->{tconv}), 
 2642               t_az_eqd, 
 2643               t_gnomonic->inverse, 
 2644               $me->_finish )
 2645       if($p->{f}); 
 2646 
 2647     $me->_finish;
 2648   }
 2649 
 2650 *t_zenithal = *t_zenithal = \&t_vertical;
 2651 
 2652 ######################################################################
 2653 
 2654 =head2 t_perspective
 2655 
 2656 =for usage
 2657 
 2658     $t = t_perspective(<options>);
 2659 
 2660 =for ref
 2661 
 2662 (Cartography) Arbitrary perspective projection 
 2663 
 2664 Perspective projection onto a focal plane from an arbitrary location
 2665 within or without the sphere, with an arbitrary central look direction,
 2666 and with correction for magnification within the optical system.
 2667 
 2668 In the forward direction, t_perspective generates perspective views of
 2669 a sphere given (lon/lat) mapping or vector information.  In the reverse
 2670 direction, t_perspective produces (lon/lat) maps from aerial or distant
 2671 photographs of spherical objects.
 2672 
 2673 Viewpoints outside the sphere treat the sphere as opaque by default,
 2674 though you can use the 'm' option to specify either the near or far
 2675 surface (relative to the origin).  Viewpoints below the surface treat
 2676 the sphere as transparent and undergo a mirror reversal for
 2677 consistency with projections that are special cases of the perspective
 2678 projection (e.g. t_gnomonic for r0=0 or t_stereographic for r0=-1).
 2679 
 2680 Magnification correction handles the extra edge distortion due to
 2681 higher angles between the focal plane and focused rays within the
 2682 optical system of your camera.  If you do not happen to know the
 2683 magnification of your camera, a simple rule of thumb is that the
 2684 magnification of a reflective telescope is roughly its focal length
 2685 (plate scale) divided by its physical length; and the magnification of 
 2686 a compound refractive telescope is roughly twice its physical length divided 
 2687 by its focal length.  Simple optical systems with a single optic have
 2688 magnification = 1.  Fisheye lenses have magnification < 1.
 2689 
 2690 This transformation was derived by direct geometrical calculation
 2691 rather than being translated from Voxland & Snyder.
 2692 
 2693 OPTIONS
 2694 
 2695 =over 3
 2696 
 2697 =item STANDARD POSITIONAL OPTIONS
 2698 
 2699 As always, the 'origin' field specifies the sub-camera point on the
 2700 sphere.
 2701 
 2702 The 'roll' option is the roll angle about the sub-camera point, for
 2703 consistency with the other projectons.
 2704 
 2705 =item p, ptg, pointing, Pointing (default (0,0,0))
 2706 
 2707 The pointing direction, in (horiz. offset, vert. offset, roll) of the
 2708 camera relative to the center of the sphere.  This is a spherical
 2709 coordinate system with the origin pointing directly at the sphere and
 2710 the pole pointing north in the pre-rolled coordinate system set by the
 2711 standard origin.  It's most useful for space-based images taken some distance
 2712 from the body in question (e.g. images of other planets or the Sun).
 2713 
 2714 Be careful not to confuse 'p' (pointing) with 'P' (P angle, a standard
 2715 synonym for roll).
 2716 
 2717 =item c, cam, camera, Camera (default undef) 
 2718 
 2719 Alternate way of specifying the camera pointing, using a spherical
 2720 coordinate system with poles at the zenith (positive) and nadir
 2721 (negative) -- this is useful for aerial photographs and such, where
 2722 the point of view is near the surface of the sphere.  You specify
 2723 (azimuth from N, altitude from horizontal, roll from vertical=up).  If
 2724 you specify pointing by this method, it overrides the 'pointing'
 2725 option, above.  This coordinate system is most useful for aerial photography
 2726 or low-orbit work, where the nadir is not necessarily the most interesting
 2727 part of the scene.
 2728 
 2729 =item r0, R0, radius, d, dist, distance [default 2.0] 
 2730 
 2731 The altitude of the point of view above the center of the sphere.
 2732 The default places the point of view 1 radius aboove the surface.
 2733 Do not confuse this with 'r', the standard origin roll angle!  Setting 
 2734 r0 < 1 gives a viewpoint inside the sphere.  In that case, the images are
 2735 mirror-reversed to preserve the chiralty of the perspective.  Setting 
 2736 r0=0 gives gnomonic projections; setting r0=-1 gives stereographic projections.
 2737 Setting r0 < -1 gives strange results.
 2738 
 2739 =item iu, im_unit, image_unit, Image_Unit (default 'degrees')
 2740 
 2741 This is the angular units in which the viewing camera is calibrated
 2742 at the center of the image.
 2743 
 2744 =item mag, magnification, Magnification (default 1.0)
 2745 
 2746 This is the magnification factor applied to the optics -- it affects the
 2747 amount of tangent-plane distortion within the telescope. 
 2748 1.0 yields the view from a simple optical system; higher values are 
 2749 telescopic, while lower values are wide-angle (fisheye).  Higher 
 2750 magnification leads to higher angles within the optical system, and more 
 2751 tangent-plane distortion at the edges of the image.  
 2752 The magnification is applied to the incident angles themselves, rather than
 2753 to their tangents (simple two-element telescopes magnify tan(theta) rather
 2754 than theta itself); this is appropriate because wide-field optics more
 2755 often conform to the equidistant azimuthal approximation than to the 
 2756 tangent plane approximation.  If you need more detailed control of 
 2757 the relationship between incident angle and focal-plane position, 
 2758 use mag=1.0 and compose the transform with something else to tweak the
 2759 angles.
 2760 
 2761 =item m, mask, Mask, h, hemisphere, Hemisphere [default 'near']
 2762 
 2763 'hemisphere' is by analogy to other cartography methods although the two 
 2764 regions to be selected are not really hemispheres.
 2765 
 2766 =item f, fov, field_of_view, Field_Of_View [default 60 degrees]
 2767 
 2768 The field of view of the telescope -- sets the crop radius on the
 2769 focal plane.  If you pass in a scalar, you get a circular crop.  If you
 2770 pass in a 2-element list ref, you get a rectilinear crop, with the
 2771 horizontal 'radius' and vertical 'radius' set separately. 
 2772 
 2773 =back
 2774 
 2775 EXAMPLES
 2776 
 2777 Model a camera looking at the Sun through a 10x telescope from Earth
 2778 (~230 solar radii from the Sun), with an 0.5 degree field of view and
 2779 a solar P (roll) angle of 30 degrees, in February (sub-Earth solar
 2780 latitude is 7 degrees south).  Convert a solar FITS image taken with
 2781 that camera to a FITS lon/lat map of the Sun with 20 pixels/degree
 2782 latitude:
 2783 
 2784   # Define map output header (no need if you don't want a FITS output map)
 2785   $maphdr = {NAXIS1=>7200,NAXIS2=>3600,            # Size of image
 2786          CTYPE1=>longitude,CTYPE2=>latitude,   # Type of axes
 2787          CUNIT1=>deg,CUNIT2=>deg,              # Unit of axes
 2788          CDELT1=>0.05,CDELT2=>0.05,            # Scale of axes
 2789          CRPIX1=>3601,CRPIX2=>1801,            # Center of map
 2790          CRVAL1=>0,CRVAL2=>0                   # (lon,lat) of center 
 2791          };
 2792 
 2793   # Set up the perspective transformation, and apply it.
 2794   $t = t_perspective(r0=>229,fov=>0.5,mag=>10,P=>30,B=>-7);
 2795   $map = $im->map( $t , $maphdr );
 2796 
 2797 Draw an aerial-view map of the Chesapeake Bay, as seen from a sounding
 2798 rocket at an altitude of 100km, looking NNE from ~200km south of
 2799 Washington (the radius of Earth is 6378 km; Washington D.C. is at
 2800 roughly 77W,38N).  Superimpose a linear coastline map on a photographic map.
 2801 
 2802   $x = graticule(1,0.1)->glue(1,earth_coast());
 2803   $t = t_perspective(r0=>6478/6378.0,fov=>60,cam=>[22.5,-20],o=>[-77,36])
 2804   $w = pgwin(size=>[10,6],J=>1);
 2805   $w->fits_imag(earth_image()->map($t,[800,500],{m=>linear}));
 2806   $w->hold;
 2807   $w->lines($x->apply($t),{xt=>'Degrees',yt=>'Degrees'});
 2808   $w->release;
 2809 
 2810 Model a 5x telescope looking at Betelgeuse with a 10 degree field of view
 2811 (since the telescope is looking at the Celestial sphere, r is 0 and this
 2812 is just an expensive modified-gnomonic projection).
 2813 
 2814   $t = t_perspective(r0=>0,fov=>10,mag=>5,o=>[88.79,7.41])
 2815 
 2816 =cut
 2817 
 2818 sub t_perspective {
 2819     my($me) = _new(@_,'Focal-Plane Perspective');
 2820     my $p = $me->{params};
 2821     
 2822     my $m= _opt($me->{options},
 2823         ['m','mask','Mask','h','hemi','hemisphere','Hemisphere'],
 2824         1);
 2825     $p->{m} = $m;
 2826     $p->{m} = 0 if($m=~m/^b/i);
 2827     $p->{m} = 1 if($m=~m/^n/i);
 2828     $p->{m} = 2 if($m=~m/^f/i);
 2829 
 2830     $p->{r0} = _opt($me->{options},
 2831             ['r0','R0','radius','Radius',
 2832              'd','dist','distance','Distance'],
 2833             2.0
 2834             );
 2835     
 2836     $p->{iu} = _opt($me->{options},
 2837            ['i','iu','image_unit','Image_Unit'],
 2838            'degrees');
 2839     
 2840     $p->{tconv} = _uconv($p->{iu});
 2841 
 2842     $p->{mag} = _opt($me->{options},
 2843              ['mag','magnification','Magnification'],
 2844              1.0);
 2845 
 2846     # Regular pointing pseudovector -- make sure there are exactly 3 elements
 2847     $p->{p} = (topdl(_opt($me->{options},
 2848             ['p','ptg','pointing','Pointing'],
 2849             [0,0,0])
 2850            )
 2851            * $p->{tconv}
 2852            )->append(zeroes(3))->slice("0:2");
 2853 
 2854     $p->{pmat} = _rotmat( (- $p->{p})->list );
 2855     
 2856     # Funky camera pointing pseudovector overrides normal pointing option 
 2857     $p->{c} = _opt($me->{options},
 2858            ['c','cam','camera','Camera'],
 2859            undef
 2860            );
 2861 
 2862     if(defined($p->{c})) {
 2863       $p->{c} = (topdl($p->{c}) * $p->{tconv})->append(zeroes(3))->slice("0:2");
 2864 
 2865       $p->{pmat} = ( _rotmat( 0,-$PI/2,0 ) x 
 2866              _rotmat( (-$p->{c})->list ) 
 2867              );
 2868     }
 2869 
 2870     # Reflect X axis if we're inside the sphere.
 2871     if($p->{r0}<1) {
 2872       $p->{pmat} = topdl([[-1,0,0],[0,1,0],[0,0,1]]) x $p->{pmat};
 2873     }
 2874 
 2875     $p->{f} = ( _opt($me->{options},
 2876              ['f','fov','field_of_view','Field_of_View'],
 2877              topdl($PI*2/3) / $p->{tconv} / $p->{mag} )
 2878         * $p->{tconv}
 2879         );
 2880     
 2881     $me->{otype} = ['Tan X','Tan Y'];
 2882     $me->{ounit} = [$p->{iu},$p->{iu}];
 2883 
 2884     # "Prefilter" -- subsidiary transform to convert the 
 2885     # spherical coordinates to 3-D coords in the viewer's 
 2886     # reference frame (Y,Z are more-or-less tangent-plane X and Y,
 2887     # and -X is the direction toward the planet, before rotation 
 2888     # to account for pointing).
 2889 
 2890     $me->{params}->{prefilt} = 
 2891       t_compose(
 2892                 # Offset for the camera pointing.
 2893         t_linear(m=>$p->{pmat},
 2894              d=>3),
 2895 
 2896                 # Rotate the sphere so the correct origin is at the 
 2897         # maximum-X point, then move the whole thing in the 
 2898         # -X direction  by r0.
 2899         t_linear(m=>(_rotmat($p->{o}->at(0),
 2900                      $p->{o}->at(1),
 2901                      $p->{roll}->at(0))
 2902                  ),
 2903              d=>3,
 2904              post=> topdl( [- $me->{params}->{r0},0,0] )
 2905              ),
 2906 
 2907                 # Put initial sci. coords into Cartesian space
 2908         t_unit_sphere(u=>'radian')  
 2909         );
 2910 
 2911     # Store the origin of the sphere -- useful for the inverse function
 2912     $me->{params}->{sph_origin} = (
 2913                    topdl([-$me->{params}->{r0},0,0]) x 
 2914                    $p->{pmat}
 2915                    )->slice(":,(0)");
 2916 
 2917     #
 2918     # Finally, the meat -- the forward function!
 2919     #
 2920     $me->{func} = sub {
 2921       my($d,$o) = @_;
 2922 
 2923       my($out) = $d->is_inplace ? $d : $d->copy;
 2924       $out->slice("0:1") *= $o->{conv};
 2925 
 2926       # If we're outside the sphere, do hemisphere filtering
 2927       my $idx;
 2928       if(abs($o->{r0}) < 1 ) {
 2929     $idx = null;
 2930       } else {
 2931     # Great-circle distance to origin
 2932     my($cos_c) = ( sin($o->{o}->slice("(1)")) * sin($out->slice("(1)"))
 2933              +
 2934              cos($o->{o}->slice("(1)")) * cos($out->slice("(1)")) *
 2935              cos($out->slice("(0)") - $o->{o}->slice("(0)"))
 2936              );
 2937 
 2938     my($thresh) = (1.0/$o->{r0});
 2939     
 2940     if($o->{m}==1) {
 2941       $idx = whichND($cos_c < $thresh);
 2942     } elsif($o->{m}==2) {
 2943       $idx = whichND($cos_c > $thresh);
 2944     } else {
 2945       $idx = null;
 2946     }
 2947       } 
 2948 
 2949       ### Transform everything -- just chuck out the bad points at the end.
 2950 
 2951       ## convert to 3-D viewer coordinates (there's a dimension change!)
 2952       my $dc = $out->apply($o->{prefilt});
 2953 
 2954       ## Apply the tangent-plane transform, and scale by the magnification.
 2955       my $dcyz = $dc->slice("1:2");
 2956 
 2957       my $r = ( $dcyz * $dcyz ) -> sumover -> sqrt ;     
 2958       my $rscale;
 2959       if( $o->{mag} == 1.0 ) {
 2960       $rscale = - 1.0 / $dc->slice("(0)");
 2961       } else {
 2962       print "(using magnification...)\n" if $PDL::verbose;
 2963       $rscale = - tan( $o->{mag} * atan( $r / $dc->slice("(0)") ) ) / $r;
 2964       }
 2965       $r *= $rscale;
 2966       $out->slice("0:1") .= $dcyz * $rscale->dummy(0,1);
 2967 
 2968       # Chuck points that are outside the FOV: glue those points
 2969       # onto the removal list.   The conditional works around a bug 
 2970       # in 2.3.4cvs and earlier: null ndarrays make append() crash.
 2971       my $w;
 2972       if(ref $o->{f} eq 'ARRAY') {
 2973     $w = whichND( ( abs($dcyz->slice("(0)")) > $o->{f}->[0] ) |
 2974               ( abs($dcyz->slice("(1)")) > $o->{f}->[1] ) |
 2975               ($r < 0)
 2976               );
 2977       } else {
 2978     $w = whichND( ($r > $o->{f}) | ($r < 0) );
 2979       }
 2980     
 2981       $idx = ($idx->nelem) ?  $idx->glue(1,$w)  : $w
 2982     if($w->nelem);
 2983 
 2984       if($idx->nelem) {
 2985     $out->slice("(0)")->range($idx) .= $o->{bad};
 2986     $out->slice("(1)")->range($idx) .= $o->{bad};
 2987       }
 2988 
 2989       ## Scale by the output conversion factor
 2990       $out->slice("0:1") /= $o->{tconv};
 2991 
 2992       $out;
 2993     };
 2994 
 2995     
 2996     #
 2997     # Inverse function
 2998     #
 2999     $me->{inv} = sub {
 3000     my($d,$o) = @_;
 3001 
 3002     my($out) = $d->is_inplace ? $d : $d->copy;
 3003     $out->slice("0:1") *= $o->{tconv};
 3004 
 3005     my $oyz = $out->slice("0:1") ;
 3006 
 3007     ## Inverse-magnify if required
 3008     if($o->{mag} != 1.0) {
 3009       my $r = ($oyz * $oyz)->sumover->sqrt;
 3010       my $scale = tan( atan( $r ) / $o->{mag} ) / $r;
 3011       $out->slice("0:1") *= $scale->dummy(0,1);
 3012     }
 3013     
 3014     ## Solve for the X coordinate of the surface.  
 3015     ## This is a quadratic in the tangent-plane coordinates;
 3016     ## so here we just figure out the coefficients and plug into
 3017     ## the quadratic formula.  $y here is actually -B/2.
 3018     my $a1 = ($oyz * $oyz)->sumover + 1;
 3019     my $y = ( $o->{sph_origin}->slice("(0)")
 3020           - ($o->{sph_origin}->slice("1:2") * $oyz)->sumover
 3021           );
 3022     my $c = topdl($o->{r0}*$o->{r0} - 1);
 3023     
 3024     my $x;
 3025     if($o->{m} == 2) { 
 3026         # Exceptional case: mask asks for the far hemisphere
 3027         $x = - ( $y - sqrt($y*$y - $a1 * $c) ) / $a1;
 3028     } else {
 3029         # normal case: mask asks for the near hemisphere
 3030         $x =   - ( $y + sqrt($y*$y - $a1 * $c) ) / $a1;
 3031     }
 3032 
 3033     ## Assemble the 3-space coordinates of the points
 3034     my $int = $out->slice("0")->append($out);
 3035     $int->sever;
 3036     $int->slice("(0)") .= -1.0;
 3037     $int->slice("0:2") *= $x->dummy(0,3);
 3038 
 3039     ## convert back to (lon,lat) coordinates...
 3040     $out .= $int->invert($o->{prefilt});    
 3041 
 3042     # If we're outside the sphere, do hemisphere filtering
 3043     my $idx;
 3044     if(abs($o->{r0}) < 1 ) {
 3045         $idx = null;
 3046     } else {
 3047         # Great-circle distance to origin
 3048         my($cos_c) = ( sin($o->{o}->slice("(1)")) * sin($out->slice("(1)"))
 3049                +
 3050                cos($o->{o}->slice("(1)")) * cos($out->slice("(1)")) *
 3051                cos($out->slice("(0)") - $o->{o}->slice("(0)"))
 3052                );
 3053         
 3054         my($thresh) = (1.0/$o->{r0});
 3055         
 3056         if($o->{m}==1) {
 3057         $idx = whichND($cos_c < $thresh);
 3058         } elsif($o->{m}==2) {
 3059         $idx = whichND($cos_c > $thresh);
 3060         }
 3061         else {
 3062         $idx = null;
 3063         }
 3064     }   
 3065 
 3066     ## Convert to the units the user requested
 3067     $out->slice("0:1") /= $o->{conv};
 3068 
 3069     ## Mark bad values
 3070     if($idx->nelem) {
 3071         $out->slice("(0)")->range($idx) .= $o->{bad};
 3072         $out->slice("(1)")->range($idx) .= $o->{bad};
 3073     }
 3074 
 3075     $out;
 3076 
 3077     };
 3078           
 3079     $me;
 3080   }
 3081 
 3082 1;