File Coverage

blib/lib/PDLA/Transform/Cartography.pm
Criterion Covered Total %
statement 142 829 17.1
branch 29 332 8.7
condition 3 48 6.2
subroutine 19 73 26.0
pod 22 23 95.6
total 215 1305 16.4


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