1Cartography(3) User Contributed Perl Documentation Cartography(3)
2
6 PDL::Transform::Cartography - Useful cartographic projections
7
9 # make a Mercator map of Earth
10 use PDL::Transform::Cartography;
11 $a = earth_coast();
12 $a = graticule(10,2)->glue(1,$a);
13 $t = t_mercator;
14 $w = pgwin(xs);
15 $w->lines($t->apply($a)->clean_lines());
16
18 PDL::Transform::Cartography includes a variety of useful cartographic
19 and observing projections (mappings of the surface of a sphere),
20 including reprojected observer coordinates. See PDL::Transform for
21 more information about image transforms in general.
22 Cartographic transformations are used for projecting not just terres‐
24 trial maps, but also any nearly spherical surface including the Sun,
25 the Celestial sphere, various moons and planets, distant stars, etc.
26 They also are useful for interpreting scientific images, which are
27 themselves generally projections of a sphere onto a flat focal plane
28 (e.g. the t_gnomonic projection).
29 Unless otherwise noted, all the transformations in this file convert
31 from (theta,phi) coordinates on the unit sphere (e.g. (lon,lat) on a
32 planet or (RA,dec) on the celestial sphere) into some sort of projected
33 coordinates, and have inverse transformations that convert back to
34 (theta,phi). This is equivalent to working from the equidistant cylin‐
35 drical (or "plate caree") projection, if you are a cartography wonk.
36 The projected coordinates are generally in units of body radii (radi‐
38 ans), so that multiplying the output by the scale of the map yields
39 physical units that are correct wherever the scale is correct for that
40 projection. For example, areas should be correct everywhere in the
41 authalic projections; and linear scales are correct along meridians in
42 the equidistant projections and along the standard parallels in all the
43 projections.
44 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 "unit
49 sphere" 3-D projection.
50 Extra dimensions tacked on to each point to be transformed are, in gen‐
52 eral, ignored. That is so that you can add on an extra index to keep
53 track of pen color. For example, earth_coast returns a 3x<n> piddle
54 containing (lon, lat, pen) at each list location. Transforming the
55 vector list retains the pen value as the first index after the dimen‐
56 sional directions.
57
59 Unless otherwise noted, the transformations and miscellaneous informa‐
60 tion in this section are taken from Snyder & Voxland 1989: "An Album of
61 Map Projections", US Geological Survey Professional Paper 1453, US
62 Printing Office (Denver); and from Snyder 1987: "Map Projections - A
63 Working Manual", US Geological Survey Professional Paper 1395, US
64 Printing Office (Denver, USA). You can obtain your own copy of both by
65 contacting the U.S. Geological Survey, Federal Center, Box 25425, Den‐
66 ver, CO 80225 USA.
67 The mathematics of cartography have a long history, and the details are
69 far trickier than the broad overview. For terrestrial (and, in gen‐
70 eral, planetary) cartography, the best reference datum is not a sphere
71 but an oblate ellipsoid due to centrifugal force from the planet's
72 rotation. Furthermore, because all rocky planets, including Earth,
73 have randomly placed mass concentrations that affect the gravitational
74 field, the reference gravitational isosurface (sea level on Earth) is
75 even more complex than an ellipsoid and, in general, different ellip‐
76 soids have been used for different locations at the same time and for
77 the same location at different times.
78 The transformations in this package use a spherical datum and hence
80 include global distortion at about the 0.5% level for terrestrial maps
81 (Earth's oblateness is ~1/300). This is roughly equal to the dimen‐
82 sional precision of physical maps printed on paper (due to stretching
83 and warping of the paper) but is significant at larger scales (e.g. for
84 regional maps). If you need more precision than that, you will want to
85 implement and use the ellipsoidal transformations from Snyder 1987 or
86 another reference work on geodesy. A good name for that package would
87 be "...::Cartography::Geodetic".
88
90 Cartographic transformations are useful for interpretation of scien‐
91 tific images, as all cameras produce projections of the celestial
92 sphere onto the focal plane of the camera. A simple (single-element)
93 optical system with a planar focal plane generates gnomonic images --
94 that is to say, gnomonic projections of a portion of the celestial
95 sphere near the paraxial direction. This is the projection that most
96 consumer grade cameras produce.
97 Magnification in an optical system changes the angle of incidence of
99 the rays on the focal plane for a given angle of incidence at the aper‐
100 ture. For example, a 10x telescope with a 2 degree field of view
101 exhibits the same gnomonic distortion as a simple optical system with a
102 20 degree field of view. Wide-angle optics typically have magnifica‐
103 tion less than 1 ('fisheye lenses'), reducing the gnomonic distortion
104 considerably but introducing "equidistant azimuthal" distortion --
105 there's no such thing as a free lunch!
106 Because many solar-system objects are spherical, PDL::Transform::Car‐
108 tography includes perspective projections for producing maps of spheri‐
109 cal bodies from perspective views. Those projections are "t_vertical"
110 and "t_perspective". They map between (lat,lon) on the spherical body
111 and planar projected coordinates at the viewpoint. "t_vertical" is the
112 vertical perspective projection given by Snyder, but "t_perspective" is
113 a fully general perspective projection that also handles magnification
114 correction.
115
117 Oblique projections rotate the sphere (and graticule) to an arbitrary
118 angle before generating the projection; transverse projections rotate
119 the sphere exactly 90 degrees before generating the projection.
120 Most of the projections accept the following standard options, useful
122 for making transverse and oblique projection maps.
123 o, origin, Origin [default (0,0,0)]
125 The origin of the oblique map coordinate system, in (old-theta,
126 old-phi) coordinates.
127 r, roll, Roll [default 0.0]
129 The roll angle of the sphere about the origin, measured CW from (N =
130 up) for reasonable values of phi and CW from (S = up) for unreason‐
131 able values of phi. This is equivalent to observer roll angle CCW
132 from the same direction.
133 u, unit, Unit [default 'degree']
135 This is the name of the angular unit to use in the lon/lat coordi‐
136 nate system.
137 b, B
139 The "B" angle of the body -- used for extraterrestrial maps. Set‐
140 ting this parameter is exactly equivalent to setting the phi compo‐
141 nent of the origin, and in fact overrides it.
142 l,L
144 The longitude of the central meridian as observed -- used for
145 extraterrestrial maps. Setting this parameter is exactly equivalent
146 to setting the theta component of the origin, and in fact overrides
147 it.
148 p,P
150 The "P" (or position) angle of the body -- used for extraterrestrial
151 maps. This parameter is a synonym for the roll angle, above.
152 bad, Bad, missing, Missing [default nan]
154 This is the value that missing points get. Mainly useful for the
155 inverse transforms. (This should work fine if set to BAD, if you
156 have bad-value support compiled in). The default nan is asin(1.2),
157 calculated at load time.
158
160 Draw a Mercator map of the world on-screen:
161 $w = pgwin(xs);
163 $w->lines(earth_coast->apply(t_mercator)->clean_lines);
164 Here, "earth_coast()" returns a 3xn piddle containing (lon, lat, pen)
166 values for the included world coastal outline; "t_mercator" converts
167 the values to projected Mercator coordinates, and "clean_lines" breaks
168 lines that cross the 180th meridian.
169 Draw a Mercator map of the world, with lon/lat at 10 degree intervals:
171 $w = pgwin(xs)
173 $a = earth_coast()->glue(1,graticule(10,1));
174 $w->lines($a->apply(t_mercator)->clean_lines);
175 This works just the same as the first example, except that a map
177 graticule has been applied with interline spacing of 10 degrees lon/lat
178 and inter-vertex spacing of 1 degree (so that each meridian contains
179 181 points, and each parallel contains 361 points).
180
182 Currently angular conversions are rather simpleminded. A list of com‐
183 mon conversions is present in the main constructor, which inserts a
184 conversion constant to radians into the {params} field of the new
185 transform. Something like Math::Convert::Units should be used instead
186 to generate the conversion constant.
187 A cleaner higher-level interface is probably needed (see the examples);
189 for example, earth_coast could return a graticule if asked, instead of
190 needing one to be glued on.
191 The class structure is somewhat messy because of the varying needs of
193 the different transformations. PDL::Transform::Cartography is a base
194 class that interprets the origin options and sets up the basic machin‐
195 ery of the Transform. The conic projections have their own subclass,
196 PDL::Transform::Conic, that interprets the standard parallels. Since
197 the cylindrical and azimuthal projections are pretty simple, they are
198 not subclassed.
199 The perl 5.6.1 compiler is quite slow at adding new classes to the
201 structure, so it does not makes sense to subclass new transformations
202 merely for the sake of pedantry.
203
205 Copyright 2002, Craig DeForest (deforest@boulder.swri.edu). This mod‐
206 ule may be modified and distributed under the same terms as PDL itself.
207 The module comes with NO WARRANTY.
208 The included digital world map is derived from the 1987 CIA World Map,
210 translated to ASCII in 1988 by Joe Dellinger (geojoe@freeusp.org) and
211 simplified in 1995 by Kirk Johnson (tuna@indra.com) for the program
212 XEarth. The map comes with NO WARRANTY. An ASCII version of the map,
213 and a sample PDL function to read it, may be found in the Demos subdi‐
214 rectory of the PDL source distribution.
215
217 The module exports both transform constructors ('t_<foo>') and some
218 auxiliary functions (no leading 't_').
219 graticule
221 $lonlatp = graticule(<grid-spacing>,<line-segment-size>);
223 (Cartography) PDL constructor - generate a lat/lon grid.
225 Returns a grid of meridians and parallels as a list of vectors suitable
227 for sending to PDL::Graphics::PGPLOT::Window::lines for plotting. The
228 grid is in degrees in (theta, phi) coordinates -- this is (E lon, N
229 lat) for terrestrial grids or (RA, dec) for celestial ones. You must
230 then transform the graticule in the same way that you transform the
231 map.
232 You can attach the graticule to a vector map using the syntax:
234 $out = graticule(10,2)->glue(1,$map);
236 In array context you get back a 2-element list containing a piddle of
238 the (theta,phi) pairs and a piddle of the pen values (1 or 0) suitable
239 for calling PDL::Graphics::PGPLOT::Window::lines. In scalar context
240 the two elements are combined into a single piddle.
241 The pen values associated with the graticule are negative, which will
243 cause PDL::Graphics::PGPLOT::Window::lines to plot them as hairlines.
244 earth_coast
246 $a = earth_coast()
248 (Cartography) PDL constructor - coastline map of Earth
250 Returns a vector coastline map based on the 1987 CIA World Coastline
252 database (see author information). The vector coastline data are in
253 plate caree format so they can be converted to other projections via
254 the apply method and cartographic transforms, and are suitable for
255 plotting with the lines method in the PGPLOT output library: the first
256 dimension is (X,Y,pen) with breaks having a pen value of 0 and hair‐
257 lines having negative pen values. The second dimension threads over
258 all the points in the data set.
259 The vector map includes lines that pass through the antipodean merid‐
261 ian, so if you want to plot it without reprojecting, you should run it
262 through clean_lines first:
263 $w = pgwin();
265 $w->lines(earth_coast->clean_lines); # plot plate caree map of world
266 $w->lines(earth_coast->apply(t_gnomonic))# plot gnomonic map of world
267 "earth_coast" is just a quick-and-dirty way of loading the file
269 "earth_coast.vec.fits" that is part of the normal installation tree.
270 earth_image
272 $rgb = earth_image()
274 (Cartography) PDL constructor - RGB pixel map of Earth
276 Returns an RGB image of Earth based on data from the MODIS instrument
278 on the NASA EOS/Terra satellite. (You can get a full-resolution image
279 from <http://earthobservatory.nasa.gov/Newsroom/BlueMarble/>). The
280 image is a plate caree map, so you can convert it to other projections
281 via the map method and cartographic transforms.
282 This is just a quick-and-dirty way of loading the earth-image files
284 that are distributed along with PDL.
285 clean_lines
287 $a = clean_lines(t_mercator->apply(scalar(earth_coast())));
289 $a = clean_lines($line_pen, [threshold]);
290 $a = $lines->clean_lines;
291 (Cartography) PDL method - remove projection irregularities
293 "clean_lines" massages vector data to remove jumps due to singularities
295 in the transform.
296 In the first (scalar) form, $line_pen contains both (X,Y) points and
298 pen values suitable to be fed to lines: in the second (list) form,
299 $lines contains the (X,Y) points and $pen contains the pen values.
300 "clean_lines" assumes that all the outline polylines are local -- that
302 is to say, there are no large jumps. Any jumps larger than a threshold
303 size are broken by setting the appropriate pen values to 0.
304 The "threshold" parameter sets the relative size of the largest jump,
306 relative to the map range (as determined by a min/max operation). The
307 default size is 0.1.
308 NOTES
310 This almost never catches stuff near the apex of cylindrical maps,
312 because the anomalous vectors get arbitrarily small. This could be
313 improved somewhat by looking at individual runs of the pen and using a
314 relative length scale that is calibrated to the rest of each run. it
315 is probably not worth the computational overhead.
316 t_unit_sphere
318 $t = t_unit_sphere(<options>);
320 (Cartography) 3-D globe projection (conformal; authalic)
322 This is similar to the inverse of t_spherical, but the inverse trans‐
324 form projects 3-D coordinates onto the unit sphere, yielding only a 2-D
325 (lon/lat) output. Similarly, the forward transform deprojects 2-D
326 (lon/lat) coordinates onto the surface of a unit sphere.
327 The cartesian system has its Z axis pointing through the pole of the
329 (lon,lat) system, and its X axis pointing through the equator at the
330 prime meridian.
331 Unit sphere mapping is unusual in that it is both conformal and
333 authalic. That is possible because it properly embeds the sphere in
334 3-space, as a notional globe.
335 This is handy as an intermediate step in lots of transforms, as Carte‐
337 sian 3-space is cleaner to work with than spherical 2-space.
338 Higher dimensional indices are preserved, so that "rider" indices (such
340 as pen value) are propagated.
341 There is no oblique transform for t_unit_sphere, largely because it's
343 so easy to rotate the output using t_linear once it's out into Carte‐
344 sian space. In fact, the other projections implement oblique trans‐
345 forms by wrapping t_linear with t_unit_sphere.
346 OPTIONS:
348 radius, Radius (default 1.0)
350 The radius of the sphere, for the inverse transform. (Radius is
351 ignored in the forward transform). Defaults to 1.0 so that the
352 resulting Cartesian coordinates are in units of "body radii".
353 t_rot_sphere
355 $t = t_rot_sphere({origin=>[<theta>,<phi>],roll=>[<roll>]});
357 (Cartography) Generate oblique projections
359 You feed in the origin in (theta,phi) and a roll angle, and you get
361 back out (theta', phi') coordinates. This is useful for making oblique
362 or transverse projections: just compose t_rot_sphere with your
363 favorite projection and you get an oblique one.
364 Most of the projections automagically compose themselves with
366 t_rot_sphere if you feed in an origin or roll angle.
367 t_rot_sphere converts the base plate caree projection (straight lon,
369 straight lat) to a Cassini projection.
370 OPTIONS
372 STANDARD POSITIONAL OPTIONS
374 t_orthographic
376 $t = t_orthographic(<options>);
378 (Cartography) Ortho. projection (azimuthal; perspective)
380 This is a perspective view as seen from infinite distance. You can
382 specify the sub-viewer point in (lon,lat) coordinates, and a rotation
383 angle of the map CW from (north=up). This is equivalent to specify
384 viewer roll angle CCW from (north=up).
385 t_orthographic is a convenience interface to t_unit_sphere -- it is
387 implemented as a composition of a t_unit_sphere call, a rotation, and a
388 slice.
389 [*] In the default case where the near hemisphere is mapped, the
391 inverse exists. There is no single inverse for the whole-sphere case,
392 so the inverse transform superimposes everything on a single hemi‐
393 sphere. If you want an invertible 3-D transform, you want
394 t_unit_sphere.
395 OPTIONS
397 STANDARD POSITIONAL OPTIONS
399 m, mask, Mask, h, hemisphere, Hemisphere [default 'near']
400 The hemisphere to keep in the projection (see PDL::Transform::Car‐
401 tography).
402 NOTES
404 Alone of the various projections, this one does not use t_rot_sphere to
406 handle the standard options, because the cartesian coordinates of the
407 rotated sphere are already correctly projected -- t_rot_sphere would
408 put them back into (theta', phi') coordinates.
409 t_caree
411 $t = t_caree(<options>);
413 (Cartography) Plate Caree projection (cylindrical; equidistant)
415 This is the simple Plate Caree projection -- also called a "lat/lon
417 plot". The horizontal axis is theta; the vertical axis is phi. This
418 is a no-op if the angular unit is radians; it is a simple scale other‐
419 wise.
420 OPTIONS
422 STANDARD POSITIONAL OPTIONS
424 s, std, standard, Standard (default 0)
425 The standard parallel where the transformation is conformal. Con‐
426 formality is achieved by shrinking of the horizontal scale to match
427 the vertical scale (which is correct everywhere).
428 t_mercator
430 $t = t_mercator(<options>);
432 (Cartography) Mercator projection (cylindrical; conformal)
434 This is perhaps the most famous of all map projections: meridians are
436 mapped to parallel vertical lines and parallels are unevenly spaced
437 horizontal lines. The poles are shifted to +/- infinity. The output
438 values are in units of globe-radii for easy conversion to kilometers;
439 hence the horizontal extent is -pi to pi.
440 You can get oblique Mercator projections by specifying the "origin" or
442 "roll" options; this is implemented via t_rot_sphere.
443 OPTIONS
445 STANDARD POSITIONAL OPTIONS
447 c, clip, Clip (default 75 [degrees])
448 The north/south clipping boundary of the transformation. Because
449 the poles are displaced to infinity, many applications require a
450 clipping boundary. The value is in whatever angular unit you set
451 with the standard 'units' option. The default roughly matches
452 interesting landforms on Earth. For no clipping at all, set b=>0.
453 For asymmetric clipping, use a 2-element list ref or piddle.
454 s, std, Standard (default 0)
456 This is the parallel at which the map has correct scale. The scale
457 is also correct at the parallel of opposite sign.
458 t_utm
460 $t = t_utm(<zone>,<options>);
462 (Cartography) Universal Transverse Mercator projection (cylindrical)
464 This is the internationally used UTM projection, with 2 subzones
466 (North/South). The UTM zones are parametrized individually, so if you
467 want a Zone 30 map you should use "t_utm(30)". By default you get the
468 northern subzone, so that locations in the southern hemisphere get neg‐
469 ative Y coordinates. If you select the southern subzone (with the
470 "subzone=>-1" option), you get offset southern UTM coordinates.
471 The 20-subzone military system is not yet supported. If/when it is
473 implemented, you will be able to enter "subzone=>[a-t]" to select a N/S
474 subzone.
475 Note that UTM is really a family of transverse Mercator projections
477 with different central meridia. Each zone properly extends for six
478 degrees of longitude on either side of its appropriate central merid‐
479 ian, with Zone 1 being centered at -177 degrees longitude (177 west).
480 Properly speaking, the zones only extend from 80 degrees south to 84
481 degrees north; but this implementation lets you go all the way to 90
482 degrees. The default UTM coordinates are meters. The origin for each
483 zone is on the equator, at an easting of -500,000 meters.
484 The default output units are meters, assuming that you are wanting a
486 map of the Earth. This will break for bodies other than Earth (which
487 have different radii and hence different conversions between lat/lon
488 angle and meters).
489 The standard UTM projection has a slight reduction in scale at the
491 prime meridian of each zone: the transverse Mercator projection's stan‐
492 dard "parallels" are 180km e/w of the central meridian. However, many
493 Europeans prefer the "Gauss-Kruger" system, which is virtually identi‐
494 cal to UTM but with a normal tangent Mercator (standard parallel on the
495 prime meridian). To get this behavior, set "gk=>1".
496 Like the rest of the PDL::Transform::Cartography package, t_utm uses a
498 spherical datum rather than the "official" ellipsoidal datums for the
499 UTM system.
500 This implementation was derived from the rather nice description by
502 Denis J. Dean, located on the web at:
503 http://www.cnr.colostate.edu/class_info/nr502/lg3/datums_coordi‐
504 nates/utm.html
505 OPTIONS
507 STANDARD OPTIONS
509 (No positional options -- Origin and Roll are ignored)
510 ou, ounit, OutputUnit (default 'meters')
512 (This is likely to become a standard option in a future release) The
513 unit of the output map. By default, this is 'meters' for UTM, but
514 you may specify 'deg' or 'km' or even (heaven help us) 'miles' if
515 you prefer.
516 sz, subzone, SubZone (default 1)
518 Set this to -1 for the southern hemisphere subzone. Ultimately you
519 should be able to set it to a letter to get the corresponding mili‐
520 tary subzone, but that's too much effort for now.
521 gk, gausskruger (default 0)
523 Set this to 1 to get the (European-style) tangent-plane Mercator
524 with standard parallel on the prime meridian. The default of 0
525 places the standard parallels 180km east/west of the prime meridian,
526 yielding better average scale across the zone. Setting gk=>1 makes
527 the scale exactly 1.0 at the central meridian, and >1.0 everywhere
528 else on the projection. The difference in scale is about 0.3%.
529 t_sin_lat
531 $t = t_sin_lat(<options>);
533 (Cartography) Cyl. equal-area projection (cyl.; authalic)
535 This projection is commonly used in solar Carrington plots; but not
537 much for terrestrial mapping.
538 OPTIONS
540 STANDARD POSITIONAL OPTIONS
542 s,std, Standard (default 0)
543 This is the parallel at which the map is conformal. It is also con‐
544 formal at the parallel of opposite sign. The conformality is
545 achieved by matched vertical stretching and horizontal squishing (to
546 achieve constant area).
547 t_sinusoidal
549 $t = t_sinusoidal(<options>);
551 (Cartography) Sinusoidal projection (authalic)
553 Sinusoidal projection preserves the latitude scale but scales longitude
555 according to sin(lat); in this respect it is the companion to
556 t_sin_lat, which is also authalic but preserves the longitude scale
557 instead.
558 OPTIONS
560 STANDARD POSITIONAL OPTIONS
562 t_conic
564 $t = t_conic(<options>)
566 (Cartography) Simple conic projection (conic; equidistant)
568 This is the simplest conic projection, with parallels mapped to
570 equidistant concentric circles. It is neither authalic nor conformal.
571 This transformation is also referred to as the "Modified Transverse
572 Mercator" projection in several maps of Alaska published by the USGS;
573 and the American State of New Mexico re-invented the projection in 1936
574 for an official map of that State.
575 OPTIONS
577 STANDARD POSITIONAL OPTIONS
579 s, std, Standard (default 29.5, 45.5)
580 The locations of the standard parallel(s) (where the cone intersects
581 the surface of the sphere). If you specify only one then the other
582 is taken to be the nearest pole. If you specify both of them to be
583 one pole then you get an equidistant azimuthal map. If you specify
584 both of them to be opposite and equidistant from the equator you get
585 a Plate Caree projection.
586 t_albers
588 $t = t_albers(<options>)
590 (Cartography) Albers conic projection (conic; authalic)
592 This is the standard projection used by the US Geological Survey for
594 sectionals of the 50 contiguous United States of America.
595 The projection reduces to the Lambert equal-area conic (infrequently
597 used and not to be confused with the Lambert conformal conic, t_lam‐
598 bert!) if the pole is used as one of the two standard parallels.
599 Notionally, this is a conic projection onto a cone that intersects the
601 sphere at the two standard parallels; it works best when the two paral‐
602 lels straddle the region of interest.
603 OPTIONS
605 STANDARD POSITIONAL OPTIONS
607 s, std, standard, Standard (default (29.5,45.5))
608 The locations of the standard parallel(s). If you specify only one
609 then the other is taken to be the nearest pole and a Lambert Equal-
610 Area Conic map results. If you specify both standard parallels to
611 be the same pole, then the projection reduces to the Lambert
612 Azimuthal Equal-Area map as aq special case. (Note that t_lambert
613 is Lambert's Conformal Conic, the most commonly used of Lambert's
614 projections.)
615 The default values for the standard parallels are those chosen by
617 Adams for maps of the lower 48 US states: (29.5,45.5). The USGS
618 recommends (55,65) for maps of Alaska and (8,18) for maps of Hawaii
619 -- these latter are chosen to also include the Canal Zone and
620 Philippine Islands farther south, which is why both of those paral‐
621 lels are south of the Hawaiian islands.
622 The transformation reduces to the cylindrical equal-area (sin-lat)
624 transformation in the case where the standard parallels are opposite
625 and equidistant from the equator, and in fact this is implemented by
626 a call to t_sin_lat.
627 t_lambert
629 $t = t_lambert(<options>);
631 (Cartography) Lambert conic projection (conic; conformal)
633 Lambert conformal conic projection is widely used in aeronautical
635 charts and state base maps published by the USA's FAA and USGS. It's
636 especially useful for mid-latitude charts. In particular, straight
637 lines approximate (but are not exactly) great circle routes of up to ~2
638 radians.
639 The default standard parallels are 33 and 45 to match the USGS state
641 1:500,000 base maps of the United States. At scales of 1:500,000 and
642 larger, discrepancies between the spherical and ellipsoidal projections
643 become important; use care with this projection on spheres.
644 OPTIONS
646 STANDARD POSITIONAL OPTIONS
648 s, std, standard, Standard (default (33,45))
649 The locations of the standard parallel(s) for the conic projection.
650 The transform reduces to the Mercator projection in the case where
651 the standard parallels are opposite and equidistant from the equa‐
652 tor, and in fact this is implemented by a call to t_mercator.
653 c, clip, Clip (default [-75,75])
655 Because the transform is conformal, the distant pole is displaced to
656 infinity. Many applications require a clipping boundary. The value
657 is in whatever angular unit you set with the standard 'unit' option.
658 For consistency with t_mercator, clipping works the same way even
659 though in most cases only one pole needs it. Set this to 0 for no
660 clipping at all.
661 t_stereographic
663 $t = t_stereographic(<options>);
665 (Cartography) Stereographic projection (az.; conf.; persp.)
667 The stereographic projection is a true perspective (planar) projection
669 from a point on the spherical surface opposite the origin of the map.
670 OPTIONS
672 STANDARD POSITIONAL OPTIONS
674 c, clip, Clip (default 120)
675 This is the angular distance from the center to the edge of the pro‐
676 jected map. The default 120 degrees gives you most of the opposite
677 hemisphere but avoids the hugely distorted part near the antipodes.
678 t_gnomonic
680 $t = t_gnomonic(<options>);
682 (Cartography) Gnomonic (focal-plane) projection (az.; persp.)
684 The gnomonic projection projects a hemisphere onto a tangent plane. It
686 is useful in cartography for the property that straight lines are great
687 circles; and it is useful in scientific imaging because it is the pro‐
688 jection generated by a simple optical system with a flat focal plane.
689 OPTIONS
691 STANDARD POSITIONAL OPTIONS
693 c, clip, Clip (default 75)
694 This is the angular distance from the center to the edge of the pro‐
695 jected map. The default 75 degrees gives you most of the hemisphere
696 but avoids the hugely distorted part near the horizon.
697 t_az_eqd
699 $t = t_az_eqd(<options>);
701 (Cartography) Azimuthal equidistant projection (az.; equi.)
703 Basic azimuthal projection preserving length along radial lines from
705 the origin (meridians, in the original polar aspect). Hence, both
706 azimuth and distance are correct for journeys beginning at the origin.
707 Applied to the celestial sphere, this is the projection made by fisheye
709 lenses; it is also the projection into which "t_vertical" puts perspec‐
710 tive views.
711 The projected plane scale is normally taken to be planetary radii; this
713 is useful for cartographers but not so useful for scientific observers.
714 Setting the 't=>1' option causes the output scale to shift to camera
715 angular coordinates (the angular unit is determined by the standard
716 'Units' option; default is degrees).
717 OPTIONS
719 STANDARD POSITIONAL OPTIONS
721 c, clip, Clip (default 180 degrees)
722 The largest angle relative to the origin. Default is the whole
723 sphere.
724 t_az_eqa
726 $t = t_az_eqa(<options>);
728 (Cartography) Azimuthal equal-area projection (az.; auth.)
730 OPTIONS
732 STANDARD POSITIONAL OPTIONS
734 c, clip, Clip (default 180 degrees)
735 The largest angle relative to the origin. Default is the whole
736 sphere.
737 t_aitoff
739 t_hammer
741 (Cartography) Hammer/Aitoff elliptical projection (az.; auth.)
743 The Hammer/Aitoff projection is often used to display the Celestial
745 sphere. It is mathematically related to the Lambert Azimuthal Equal-
746 Area projection (t_az_eqa), and maps the sphere to an ellipse of unit
747 eccentricity, with vertical radius sqrt(2) and horizontal radius of 2
748 sqrt(2).
749 OPTIONS
751 STANDARD POSITIONAL OPTIONS
753 t_vertical
755 $t = t_vertical(<options>);
757 (Cartography) Vertical perspective projection (az.; persp.)
759 Vertical perspective projection is a generalization of gnomonic and
761 stereographic projection, and a special case of perspective projection.
762 It is a projection from the sphere onto a focal plane at the camera
763 location.
764 OPTIONS
766 STANDARD POSITIONAL OPTIONS
768 m, mask, Mask, h, hemisphere, Hemisphere [default 'near']
769 The hemisphere to keep in the projection (see PDL::Transform::Car‐
770 tography).
771 r0, R0, radius, d, dist, distance [default 2.0]
773 The altitude of the focal plane above the center of the sphere. The
774 default places the point of view one radius above the surface.
775 t, telescope, Telescope, cam, Camera (default '')
777 If this is set, then the central scale is in telescope or camera
778 angular units rather than in planetary radii. The angular units are
779 parsed as with the normal 'u' option for the lon/lat specification.
780 If you specify a non-string value (such as 1) then you get tele‐
781 scope-frame radians, suitable for working on with other transforma‐
782 tions.
783 f, fish, fisheye (default '')
785 If this is set then the output is in azimuthal equidistant coordi‐
786 nates instead of in tangent-plane coordinates. This is a conve‐
787 nience function for '(t_az_eqd) x !(t_gnomonic) x (t_vertical)'.
788 t_perspective
790 $t = t_perspective(<options>);
792 (Cartography) Arbitrary perspective projection
794 Perspective projection onto a focal plane from an arbitrary location
796 within or without the sphere, with an arbitary central look direction,
797 and with correction for magnification within the optical system.
798 In the forward direction, t_perspective generates perspective views of
800 a sphere given (lon/lat) mapping or vector information. In the reverse
801 direction, t_perspective produces (lon/lat) maps from aerial or distant
802 photographs of spherical objects.
803 Viewpoints outside the sphere treat the sphere as opaque by default,
805 though you can use the 'm' option to specify either the near or far
806 surface (relative to the origin). Viewpoints below the surface treat
807 the sphere as transparent and undergo a mirror reversal for consistency
808 with projections that are special cases of the perspective projection
809 (e.g. t_gnomonic for r0=0 or t_stereographic for r0=-1).
810 Magnification correction handles the extra edge distortion due to
812 higher angles between the focal plane and focused rays within the opti‐
813 cal system of your camera. If you do not happen to know the magnifica‐
814 tion of your camera, a simple rule of thumb is that the magnification
815 of a reflective telescope is roughly its focal length (plate scale)
816 divided by its physical length; and the magnification of a compound
817 refractive telescope is roughly twice its physical length divided by
818 its focal length. Simple optical sytems with a single optic have mag‐
819 nification = 1. Fisheye lenses have magnification < 1.
820 This transformation was derived by direct geometrical calculation
822 rather than being translated from Voxland & Snyder.
823 OPTIONS
825 STANDARD POSITIONAL OPTIONS
827 As always, the 'origin' field specifies the sub-camera point on the
828 sphere.
829 The 'roll' option is the roll angle about the sub-camera point, for
831 consistency with the other projectons.
832 p, ptg, pointing, Pointing (default (0,0,0))
834 The pointing direction, in (horiz. offset, vert. offset, roll) of
835 the camera relative to the center of the sphere. This is a spheri‐
836 cal coordinate system with the origin pointing directly at the
837 sphere and the pole pointing north in the pre-rolled coordinate sys‐
838 tem set by the standard origin. It's most useful for space-based
839 images taken some distance from the body in question (e.g. images of
840 other planets or the Sun).
841 Be careful not to confuse 'p' (pointing) with 'P' (P angle, a stan‐
843 dard synonym for roll).
844 c, cam, camera, Camera (default undef)
846 Alternate way of specifying the camera pointing, using a spherical
847 coordinate system with poles at the zenith (positive) and nadir
848 (negative) -- this is useful for aerial photographs and such, where
849 the point of view is near the surface of the sphere. You specify
850 (azimuth from N, altitude from horizontal, roll from vertical=up).
851 If you specify pointing by this method, it overrides the 'pointing'
852 option, above. This coordinate system is most useful for aerial
853 photography or low-orbit work, where the nadir is not necessarily
854 the most interesting part of the scene.
855 r0, R0, radius, d, dist, distance [default 2.0]
857 The altitude of the point of view above the center of the sphere.
858 The default places the point of view 1 radius aboove the surface.
859 Do not confuse this with 'r', the standard origin roll angle! Set‐
860 ting r0 < 1 gives a viewpoint inside the sphere. In that case, the
861 images are mirror-reversed to preserve the chiralty of the perspec‐
862 tive. Setting r0=0 gives gnomonic projections; setting r0=-1 gives
863 stereographic projections. Setting r0 < -1 gives strange results.
864 iu, im_unit, image_unit, Image_Unit (default 'degrees')
866 This is the angular units in which the viewing camera is calibrated
867 at the center of the image.
868 mag, magnification, Magnification (default 1.0)
870 This is the magnification factor applied to the optics -- it affects
871 the amount of tangent-plane distortion within the telescope. 1.0
872 yields the view from a simple optical system; higher values are
873 telescopic, while lower values are wide-angle (fisheye). Higher
874 magnification leads to higher angles within the optical system, and
875 more tangent-plane distortion at the edges of the image. The magni‐
876 fication is applied to the incident angles themselves, rather than
877 to their tangents (simple two-element telescopes magnify tan(theta)
878 rather than theta itself); this is appropriate because wide-field
879 optics more often conform to the equidistant azimuthal approximation
880 than to the tangent plane approximation. If you need more detailed
881 control of the relationship between incident angle and focal-plane
882 position, use mag=1.0 and compose the transform with something else
883 to tweak the angles.
884 m, mask, Mask, h, hemisphere, Hemisphere [default 'near']
886 'hemisphere' is by analogy to other cartography methods although the
887 two regions to be selected are not really hemispheres.
888 f, fov, field_of_view, Field_Of_View [default 60 degrees]
890 The field of view of the telescope -- sets the crop radius on the
891 focal plane. If you pass in a scalar, you get a circular crop. If
892 you pass in a 2-element list ref, you get a rectilinear crop, with
893 the horizontal 'radius' and vertical 'radius' set separately.
894 EXAMPLES
896 Model a camera looking at the Sun through a 10x telescope from Earth
898 (~230 solar radii from the Sun), with an 0.5 degree field of view and a
899 solar P (roll) angle of 30 degrees, in February (sub-Earth solar lati‐
900 tude is 7 degrees south). Convert a solar FITS image taken with that
901 camera to a FITS lon/lat map of the Sun with 20 pixels/degree latitude:
902 # Define map output header (no need if you don't want a FITS output map)
904 $maphdr = {NAXIS1=>7200,NAXIS2=>3600, # Size of image
905 CTYPE1=>longitude,CTYPE2=>latitude, # Type of axes
906 CUNIT1=>deg,CUNIT2=>deg, # Unit of axes
907 CDELT1=>0.05,CDELT2=>0.05, # Scale of axes
908 CRPIX1=>3601,CRPIX2=>1801, # Center of map
909 CRVAL1=>0,CRVAL2=>0 # (lon,lat) of center
910 };
911 # Set up the perspective transformation, and apply it.
913 $t = t_perspective(r0=>229,fov=>0.5,mag=>10,P=>30,B=>-7);
914 $map = $im->map( $t , $maphdr );
915 Draw an aerial-view map of the Chesapeake Bay, as seen from a sounding
917 rocket at an altitude of 100km, looking NNE from ~200km south of Wash‐
918 ington (the radius of Earth is 6378 km; Washington D.C. is at roughly
919 77W,38N). Superimpose a linear coastline map on a photographic map.
920 $a = graticule(1,0.1)->glue(1,earth_coast());
922 $t = t_perspective(r0=>6478/6378.0,fov=>60,cam=>[22.5,-20],o=>[-77,36])
923 $w = pgwin(size=>[10,6],J=>1);
924 $w->fits_imag(earth_image()->map($t,[800,500],{m=>linear}));
925 $w->hold;
926 $w->lines($a->apply($t),{xt=>'Degrees',yt=>'Degrees'});
927 $w->release;
928 Model a 5x telescope looking at Betelgeuse with a 10 degree field of
930 view (since the telescope is looking at the Celestial sphere, r is 0
931 and this is just an expensive modified-gnomonic projection).
932 $t = t_perspective(r0=>0,fov=>10,mag=>5,o=>[88.79,7.41])
934perl v5.8.8 2004-07-21 Cartography(3)