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 $x = earth_coast();
12 $x = graticule(10,2)->glue(1,$x);
13 $t = t_mercator;
14 $w = pgwin(xs);
15 $w->lines($t->apply($x)->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
24 terrestrial maps, but also any nearly spherical surface including the
25 Sun, the Celestial sphere, various moons and planets, distant stars,
26 etc. They also are useful for interpreting scientific images, which
27 are themselves generally projections of a sphere onto a flat focal
28 plane (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
35 cylindrical (or "plate carree") projection, if you are a cartography
36 wonk.
37 The projected coordinates are generally in units of body radii
39 (radians), so that multiplying the output by the scale of the map
40 yields physical units that are correct wherever the scale is correct
41 for that projection. For example, areas should be correct everywhere
42 in the authalic projections; and linear scales are correct along
43 meridians in the equidistant projections and along the standard
44 parallels in all the projections.
45 The transformations that are authalic (equal-area), conformal (equal-
47 angle), azimuthal (circularly symmetric), or perspective (true
48 perspective on a focal plane from some viewpoint) are marked. The
49 first two categories are mutually exclusive for all but the "unit
50 sphere" 3-D projection.
51 Extra dimensions tacked on to each point to be transformed are, in
53 general, ignored. That is so that you can add on an extra index to
54 keep track of pen color. For example, "earth_coast" returns a 3x<n>
55 ndarray containing (lon, lat, pen) at each list location. Transforming
56 the vector list retains the pen value as the first index after the
57 dimensional directions.
58
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 The mathematics of cartography have a long history, and the details are
70 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 same
78 time and for the same location at different times.
79 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 that,
86 you will want to implement and use the ellipsoidal transformations from
87 Snyder 1987 or another reference work on geodesy. A good name for that
88 package would be "...::Cartography::Geodetic".
89
91 Cartographic transformations are useful for interpretation of
92 scientific images, as all cameras produce projections of the celestial
93 sphere onto the focal plane of the camera. A simple (single-element)
94 optical system with a planar focal plane generates gnomonic images --
95 that is to say, gnomonic projections of a portion of the celestial
96 sphere near the paraxial direction. This is the projection that most
97 consumer grade cameras produce.
98 Magnification in an optical system changes the angle of incidence of
100 the rays on the focal plane for a given angle of incidence at the
101 aperture. For example, a 10x telescope with a 2 degree field of view
102 exhibits the same gnomonic distortion as a simple optical system with a
103 20 degree field of view. Wide-angle optics typically have
104 magnification less than 1 ('fisheye lenses'), reducing the gnomonic
105 distortion considerably but introducing "equidistant azimuthal"
106 distortion -- there's no such thing as a free lunch!
107 Because many solar-system objects are spherical,
109 PDL::Transform::Cartography includes perspective projections for
110 producing maps of spherical bodies from perspective views. Those
111 projections are "t_vertical" and "t_perspective". They map between
112 (lat,lon) on the spherical body and planar projected coordinates at the
113 viewpoint. "t_vertical" is the vertical perspective projection given
114 by Snyder, but "t_perspective" is a fully general perspective
115 projection that also handles magnification correction.
116
118 Oblique projections rotate the sphere (and graticule) to an arbitrary
119 angle before generating the projection; transverse projections rotate
120 the sphere exactly 90 degrees before generating the projection.
121 Most of the projections accept the following standard options, useful
123 for making transverse and oblique projection maps.
124 o, origin, Origin [default (0,0,0)]
126 The origin of the oblique map coordinate system, in (old-theta, old-
127 phi) coordinates.
128 r, roll, Roll [default 0.0]
130 The roll angle of the sphere about the origin, measured CW from (N =
131 up) for reasonable values of phi and CW from (S = up) for
132 unreasonable values of phi. This is equivalent to observer roll
133 angle CCW from the same direction.
134 u, unit, Unit [default 'degree']
136 This is the name of the angular unit to use in the lon/lat
137 coordinate system.
138 b, B
140 The "B" angle of the body -- used for extraterrestrial maps.
141 Setting this parameter is exactly equivalent to setting the phi
142 component of the origin, and in fact overrides it.
143 l,L
145 The longitude of the central meridian as observed -- used for
146 extraterrestrial maps. Setting this parameter is exactly equivalent
147 to setting the theta component of the origin, and in fact overrides
148 it.
149 p,P
151 The "P" (or position) angle of the body -- used for extraterrestrial
152 maps. This parameter is a synonym for the roll angle, above.
153 bad, Bad, missing, Missing [default nan]
155 This is the value that missing points get. Mainly useful for the
156 inverse transforms. (This should work fine if set to BAD, if you
157 have bad-value support compiled in). The default nan is asin(1.2),
158 calculated at load time.
159
161 Draw a Mercator map of the world on-screen:
162 $w = pgwin(xs);
164 $w->lines(earth_coast->apply(t_mercator)->clean_lines);
165 Here, earth_coast() returns a 3xn ndarray containing (lon, lat, pen)
167 values for the included world coastal outline; "t_mercator" converts
168 the values to projected Mercator coordinates, and "clean_lines" breaks
169 lines that cross the 180th meridian.
170 Draw a Mercator map of the world, with lon/lat at 10 degree intervals:
172 $w = pgwin(xs)
174 $x = earth_coast()->glue(1,graticule(10,1));
175 $w->lines($x->apply(t_mercator)->clean_lines);
176 This works just the same as the first example, except that a map
178 graticule has been applied with interline spacing of 10 degrees lon/lat
179 and inter-vertex spacing of 1 degree (so that each meridian contains
180 181 points, and each parallel contains 361 points).
181
183 Currently angular conversions are rather simpleminded. A list of
184 common conversions is present in the main constructor, which inserts a
185 conversion constant to radians into the {params} field of the new
186 transform. Something like Math::Convert::Units should be used instead
187 to generate the conversion constant.
188 A cleaner higher-level interface is probably needed (see the examples);
190 for example, earth_coast could return a graticule if asked, instead of
191 needing one to be glued on.
192 The class structure is somewhat messy because of the varying needs of
194 the different transformations. PDL::Transform::Cartography is a base
195 class that interprets the origin options and sets up the basic
196 machinery of the Transform. The conic projections have their own
197 subclass, PDL::Transform::Conic, that interprets the standard
198 parallels. Since the cylindrical and azimuthal projections are pretty
199 simple, they are not subclassed.
200 The perl 5.6.1 compiler is quite slow at adding new classes to the
202 structure, so it does not makes sense to subclass new transformations
203 merely for the sake of pedantry.
204
206 Copyright 2002, Craig DeForest (deforest@boulder.swri.edu). This
207 module may be modified and distributed under the same terms as PDL
208 itself. The module comes with NO WARRANTY.
209 The included digital world map is derived from the 1987 CIA World Map,
211 translated to ASCII in 1988 by Joe Dellinger (geojoe@freeusp.org) and
212 simplified in 1995 by Kirk Johnson (tuna@indra.com) for the program
213 XEarth. The map comes with NO WARRANTY. An ASCII version of the map,
214 and a sample PDL function to read it, may be found in the Demos
215 subdirectory of the PDL source distribution.
216
218 The module exports both transform constructors ('t_<foo>') and some
219 auxiliary functions (no leading 't_').
220 graticule
222 $lonlatp = graticule(<grid-spacing>,<line-segment-size>);
223 $lonlatp = graticule(<grid-spacing>,<line-segment-size>,1);
224 (Cartography) PDL constructor - generate a lat/lon grid.
226 Returns a grid of meridians and parallels as a list of vectors suitable
228 for sending to PDL::Graphics::PGPLOT::Window::lines for plotting. The
229 grid is in degrees in (theta, phi) coordinates -- this is (E lon, N
230 lat) for terrestrial grids or (RA, dec) for celestial ones. You must
231 then transform the graticule in the same way that you transform the
232 map.
233 You can attach the graticule to a vector map using the syntax:
235 $out = graticule(10,2)->glue(1,$map);
237 In array context you get back a 2-element list containing an ndarray of
239 the (theta,phi) pairs and an ndarray of the pen values (1 or 0)
240 suitable for calling PDL::Graphics::PGPLOT::Window::lines. In scalar
241 context the two elements are combined into a single ndarray.
242 The pen values associated with the graticule are negative, which will
244 cause PDL::Graphics::PGPLOT::Window::lines to plot them as hairlines.
245 If a third argument is given, it is a hash of options, which can be:
247 nan - if true, use two columns instead of three, and separate lines
249 with a 'nan' break
250 lonpos - if true, all reported longitudes are positive (0 to 360)
251 instead of (-180 to 180).
252 dup - if true, the meridian at the far boundary is duplicated.
253 earth_coast
255 $x = earth_coast()
256 (Cartography) PDL constructor - coastline map of Earth
258 Returns a vector coastline map based on the 1987 CIA World Coastline
260 database (see author information). The vector coastline data are in
261 plate carree format so they can be converted to other projections via
262 the apply method and cartographic transforms, and are suitable for
263 plotting with the lines method in the PGPLOT output library: the first
264 dimension is (X,Y,pen) with breaks having a pen value of 0 and
265 hairlines having negative pen values. The second dimension broadcasts
266 over all the points in the data set.
267 The vector map includes lines that pass through the antipodean
269 meridian, so if you want to plot it without reprojecting, you should
270 run it through "clean_lines" first:
271 $w = pgwin();
273 $w->lines(earth_coast->clean_lines); # plot plate carree map of world
274 $w->lines(earth_coast->apply(t_gnomonic))# plot gnomonic map of world
275 "earth_coast" is just a quick-and-dirty way of loading the file
277 "earth_coast.vec.fits" that is part of the normal installation tree.
278 earth_image
280 $rgb = earth_image()
281 (Cartography) PDL constructor - RGB pixel map of Earth
283 Returns an RGB image of Earth based on data from the MODIS instrument
285 on the NASA EOS/Terra satellite. (You can get a full-resolution image
286 from <http://earthobservatory.nasa.gov/Newsroom/BlueMarble/>). The
287 image is a plate carree map, so you can convert it to other projections
288 via the map method and cartographic transforms.
289 This is just a quick-and-dirty way of loading the earth-image files
291 that are distributed along with PDL.
292 earth_shape
294 $fits_shape = earth_shape()
295 (Cartography) PDL constructor - height map of Earth
297 Returns a height map of Earth based on data from the General
299 Bathymetric Chart of the Oceans (GEBCO) produced by the British
300 Oceanographic Data Centre. (You can get a full-resolution image from
301 <http://visibleearth.nasa.gov/view.php?id=73934>). The image is a
302 plate carree map, so you can convert it to other projections via the
303 map method and cartographic transforms. The data is from 8-bit
304 grayscale (so only 256 levels), but is returned in a similar format to
305 "earth_image". The range represents a span of 6400m, so Everest and the
306 Marianas Trench are not accurately represented.
307 To turn this into a "float", ("lonlatradius,x,y") with "x" and "y" in
309 radians, and the radius as a "float" as a proportion of the Earth's
310 mean radius, use "t_raster2float". The Earth is treated here as a
311 perfect sphere with sea level at radius 6,371km.
312 Value Hex value Float From centre in km Float as radius
314 Base 00 0.0 6370.69873km 0.99995
315 Sea level 0C 0.04705 6371km 1.0
316 Highest FF 1.0 6377.09863km 1.00096
317 Code:
319 $shape = earth_shape();
321 $floats = t_raster2float()->apply($shape->mv(2,0));
322 $lonlatradius = $floats->slice('0:2'); # r g b all same
323 $lonlatradius->slice('(2)') *= float((6377.09863 - 6370.69873) / 6371);
324 $lonlatradius->slice('(2)') += float(6370.69873 / 6371);
325 clean_lines
327 $x = clean_lines(t_mercator->apply(scalar(earth_coast())));
328 $x = clean_lines($line_pen, [threshold]);
329 $x = $lines->clean_lines;
330 (Cartography) PDL method - remove projection irregularities
332 "clean_lines" massages vector data to remove jumps due to singularities
334 in the transform.
335 In the first (scalar) form, $line_pen contains both (X,Y) points and
337 pen values suitable to be fed to lines: in the second (list) form,
338 $lines contains the (X,Y) points and $pen contains the pen values.
339 "clean_lines" assumes that all the outline polylines are local -- that
341 is to say, there are no large jumps. Any jumps larger than a threshold
342 size are broken by setting the appropriate pen values to 0.
343 The "threshold" parameter sets the relative size of the largest jump,
345 relative to the map range (as determined by a min/max operation). The
346 default size is 0.1.
347 NOTES
349 This almost never catches stuff near the apex of cylindrical maps,
351 because the anomalous vectors get arbitrarily small. This could be
352 improved somewhat by looking at individual runs of the pen and using a
353 relative length scale that is calibrated to the rest of each run. it
354 is probably not worth the computational overhead.
355 t_raster2float
357 $t = t_raster2float();
358 (Cartography) Convert a raster (3,x,y) to "float" (lonlatrgb,x,y)
360 Assumes "bytes" input, and radians and "float" output, with the first 2
362 coordinates suitable for use as plate carree.
363 t_raster2fits
365 $t = t_raster2fits();
366 (Cartography) Convert a raster (3,x,y) to FITS plate carree (x,y,3)
368 Adds suitable "hdr". Assumes degrees. Used by "earth_image".
370 t_unit_sphere
372 $t = t_unit_sphere(<options>);
373 (Cartography) 3-D globe projection (conformal; authalic)
375 This is similar to the inverse of t_spherical, but the inverse
377 transform projects 3-D coordinates onto the unit sphere, yielding only
378 a 2-D (lon/lat) output. Similarly, the forward transform deprojects
379 2-D (lon/lat) coordinates onto the surface of a unit sphere.
380 The cartesian system has its Z axis pointing through the pole of the
382 (lon,lat) system, and its X axis pointing through the equator at the
383 prime meridian.
384 Unit sphere mapping is unusual in that it is both conformal and
386 authalic. That is possible because it properly embeds the sphere in
387 3-space, as a notional globe.
388 This is handy as an intermediate step in lots of transforms, as
390 Cartesian 3-space is cleaner to work with than spherical 2-space.
391 Higher dimensional indices are preserved, so that "rider" indices (such
393 as pen value) are propagated.
394 There is no oblique transform for t_unit_sphere, largely because it's
396 so easy to rotate the output using t_linear once it's out into
397 Cartesian space. In fact, the other projections implement oblique
398 transforms by wrapping t_linear with "t_unit_sphere".
399 OPTIONS:
401 radius, Radius (default 1.0)
403 The radius of the sphere, for the inverse transform. (Radius is
404 ignored in the forward transform). Defaults to 1.0 so that the
405 resulting Cartesian coordinates are in units of "body radii".
406 t_rot_sphere
408 $t = t_rot_sphere({origin=>[<theta>,<phi>],roll=>[<roll>]});
409 (Cartography) Generate oblique projections
411 You feed in the origin in (theta,phi) and a roll angle, and you get
413 back out (theta', phi') coordinates. This is useful for making oblique
414 or transverse projections: just compose t_rot_sphere with your
415 favorite projection and you get an oblique one.
416 Most of the projections automagically compose themselves with
418 t_rot_sphere if you feed in an origin or roll angle.
419 t_rot_sphere converts the base plate carree projection (straight lon,
421 straight lat) to a Cassini projection.
422 OPTIONS
424 STANDARD POSITIONAL OPTIONS
426 t_orthographic
428 $t = t_orthographic(<options>);
429 (Cartography) Ortho. projection (azimuthal; perspective)
431 This is a perspective view as seen from infinite distance. You can
433 specify the sub-viewer point in (lon,lat) coordinates, and a rotation
434 angle of the map CW from (north=up). This is equivalent to specify
435 viewer roll angle CCW from (north=up).
436 t_orthographic is a convenience interface to t_unit_sphere -- it is
438 implemented as a composition of a t_unit_sphere call, a rotation, and a
439 slice.
440 [*] In the default case where the near hemisphere is mapped, the
442 inverse exists. There is no single inverse for the whole-sphere case,
443 so the inverse transform superimposes everything on a single
444 hemisphere. If you want an invertible 3-D transform, you want
445 "t_unit_sphere".
446 OPTIONS
448 STANDARD POSITIONAL OPTIONS
450 m, mask, Mask, h, hemisphere, Hemisphere [default 'near']
451 The hemisphere to keep in the projection (see
452 PDL::Transform::Cartography).
453 NOTES
455 Alone of the various projections, this one does not use "t_rot_sphere"
457 to handle the standard options, because the cartesian coordinates of
458 the rotated sphere are already correctly projected -- t_rot_sphere
459 would put them back into (theta', phi') coordinates.
460 t_carree
462 $t = t_carree(<options>);
463 (Cartography) Plate Carree projection (cylindrical; equidistant)
465 This is the simple Plate Carree projection -- also called a "lat/lon
467 plot". The horizontal axis is theta; the vertical axis is phi. This
468 is a no-op if the angular unit is radians; it is a simple scale
469 otherwise.
470 OPTIONS
472 STANDARD POSITIONAL OPTIONS
474 s, std, standard, Standard (default 0)
475 The standard parallel where the transformation is conformal.
476 Conformality is achieved by shrinking of the horizontal scale to
477 match the vertical scale (which is correct everywhere).
478 t_mercator
480 $t = t_mercator(<options>);
481 (Cartography) Mercator projection (cylindrical; conformal)
483 This is perhaps the most famous of all map projections: meridians are
485 mapped to parallel vertical lines and parallels are unevenly spaced
486 horizontal lines. The poles are shifted to +/- infinity. The output
487 values are in units of globe-radii for easy conversion to kilometers;
488 hence the horizontal extent is -pi to pi.
489 You can get oblique Mercator projections by specifying the "origin" or
491 "roll" options; this is implemented via "t_rot_sphere".
492 OPTIONS
494 STANDARD POSITIONAL OPTIONS
496 c, clip, Clip (default 75 [degrees])
497 The north/south clipping boundary of the transformation. Because
498 the poles are displaced to infinity, many applications require a
499 clipping boundary. The value is in whatever angular unit you set
500 with the standard 'units' option. The default roughly matches
501 interesting landforms on Earth. For no clipping at all, set b=>0.
502 For asymmetric clipping, use a 2-element list ref or ndarray.
503 s, std, Standard (default 0)
505 This is the parallel at which the map has correct scale. The scale
506 is also correct at the parallel of opposite sign.
507 t_utm
509 $t = t_utm(<zone>,<options>);
510 (Cartography) Universal Transverse Mercator projection (cylindrical)
512 This is the internationally used UTM projection, with 2 subzones
514 (North/South). The UTM zones are parametrized individually, so if you
515 want a Zone 30 map you should use t_utm(30). By default you get the
516 northern subzone, so that locations in the southern hemisphere get
517 negative Y coordinates. If you select the southern subzone (with the
518 "subzone=>-1" option), you get offset southern UTM coordinates.
519 The 20-subzone military system is not yet supported. If/when it is
521 implemented, you will be able to enter "subzone=>[a-t]" to select a N/S
522 subzone.
523 Note that UTM is really a family of transverse Mercator projections
525 with different central meridia. Each zone properly extends for six
526 degrees of longitude on either side of its appropriate central
527 meridian, with Zone 1 being centered at -177 degrees longitude (177
528 west). Properly speaking, the zones only extend from 80 degrees south
529 to 84 degrees north; but this implementation lets you go all the way to
530 90 degrees. The default UTM coordinates are meters. The origin for
531 each zone is on the equator, at an easting of -500,000 meters.
532 The default output units are meters, assuming that you are wanting a
534 map of the Earth. This will break for bodies other than Earth (which
535 have different radii and hence different conversions between lat/lon
536 angle and meters).
537 The standard UTM projection has a slight reduction in scale at the
539 prime meridian of each zone: the transverse Mercator projection's
540 standard "parallels" are 180km e/w of the central meridian. However,
541 many Europeans prefer the "Gauss-Kruger" system, which is virtually
542 identical to UTM but with a normal tangent Mercator (standard parallel
543 on the prime meridian). To get this behavior, set "gk=>1".
544 Like the rest of the PDL::Transform::Cartography package, t_utm uses a
546 spherical datum rather than the "official" ellipsoidal datums for the
547 UTM system.
548 This implementation was derived from the rather nice description by
550 Denis J. Dean, located on the web at:
551 http://www.cnr.colostate.edu/class_info/nr502/lg3/datums_coordinates/utm.html
552 OPTIONS
554 STANDARD OPTIONS
556 (No positional options -- Origin and Roll are ignored)
557 ou, ounit, OutputUnit (default 'meters')
559 (This is likely to become a standard option in a future release) The
560 unit of the output map. By default, this is 'meters' for UTM, but
561 you may specify 'deg' or 'km' or even (heaven help us) 'miles' if
562 you prefer.
563 sz, subzone, SubZone (default 1)
565 Set this to -1 for the southern hemisphere subzone. Ultimately you
566 should be able to set it to a letter to get the corresponding
567 military subzone, but that's too much effort for now.
568 gk, gausskruger (default 0)
570 Set this to 1 to get the (European-style) tangent-plane Mercator
571 with standard parallel on the prime meridian. The default of 0
572 places the standard parallels 180km east/west of the prime meridian,
573 yielding better average scale across the zone. Setting gk=>1 makes
574 the scale exactly 1.0 at the central meridian, and >1.0 everywhere
575 else on the projection. The difference in scale is about 0.3%.
576 t_sin_lat
578 $t = t_sin_lat(<options>);
579 (Cartography) Cyl. equal-area projection (cyl.; authalic)
581 This projection is commonly used in solar Carrington plots; but not
583 much for terrestrial mapping.
584 OPTIONS
586 STANDARD POSITIONAL OPTIONS
588 s,std, Standard (default 0)
589 This is the parallel at which the map is conformal. It is also
590 conformal at the parallel of opposite sign. The conformality is
591 achieved by matched vertical stretching and horizontal squishing (to
592 achieve constant area).
593 t_sinusoidal
595 $t = t_sinusoidal(<options>);
596 (Cartography) Sinusoidal projection (authalic)
598 Sinusoidal projection preserves the latitude scale but scales longitude
600 according to sin(lat); in this respect it is the companion to
601 "t_sin_lat", which is also authalic but preserves the longitude scale
602 instead.
603 OPTIONS
605 STANDARD POSITIONAL OPTIONS
607 t_conic
609 $t = t_conic(<options>)
610 (Cartography) Simple conic projection (conic; equidistant)
612 This is the simplest conic projection, with parallels mapped to
614 equidistant concentric circles. It is neither authalic nor conformal.
615 This transformation is also referred to as the "Modified Transverse
616 Mercator" projection in several maps of Alaska published by the USGS;
617 and the American State of New Mexico re-invented the projection in 1936
618 for an official map of that State.
619 OPTIONS
621 STANDARD POSITIONAL OPTIONS
623 s, std, Standard (default 29.5, 45.5)
624 The locations of the standard parallel(s) (where the cone intersects
625 the surface of the sphere). If you specify only one then the other
626 is taken to be the nearest pole. If you specify both of them to be
627 one pole then you get an equidistant azimuthal map. If you specify
628 both of them to be opposite and equidistant from the equator you get
629 a Plate Carree projection.
630 t_albers
632 $t = t_albers(<options>)
633 (Cartography) Albers conic projection (conic; authalic)
635 This is the standard projection used by the US Geological Survey for
637 sectionals of the 50 contiguous United States of America.
638 The projection reduces to the Lambert equal-area conic (infrequently
640 used and not to be confused with the Lambert conformal conic,
641 "t_lambert"!) if the pole is used as one of the two standard
642 parallels.
643 Notionally, this is a conic projection onto a cone that intersects the
645 sphere at the two standard parallels; it works best when the two
646 parallels straddle the region of interest.
647 OPTIONS
649 STANDARD POSITIONAL OPTIONS
651 s, std, standard, Standard (default (29.5,45.5))
652 The locations of the standard parallel(s). If you specify only one
653 then the other is taken to be the nearest pole and a Lambert Equal-
654 Area Conic map results. If you specify both standard parallels to
655 be the same pole, then the projection reduces to the Lambert
656 Azimuthal Equal-Area map as aq special case. (Note that "t_lambert"
657 is Lambert's Conformal Conic, the most commonly used of Lambert's
658 projections.)
659 The default values for the standard parallels are those chosen by
661 Adams for maps of the lower 48 US states: (29.5,45.5). The USGS
662 recommends (55,65) for maps of Alaska and (8,18) for maps of Hawaii
663 -- these latter are chosen to also include the Canal Zone and
664 Philippine Islands farther south, which is why both of those
665 parallels are south of the Hawaiian islands.
666 The transformation reduces to the cylindrical equal-area (sin-lat)
668 transformation in the case where the standard parallels are opposite
669 and equidistant from the equator, and in fact this is implemented by
670 a call to t_sin_lat.
671 t_lambert
673 $t = t_lambert(<options>);
674 (Cartography) Lambert conic projection (conic; conformal)
676 Lambert conformal conic projection is widely used in aeronautical
678 charts and state base maps published by the USA's FAA and USGS. It's
679 especially useful for mid-latitude charts. In particular, straight
680 lines approximate (but are not exactly) great circle routes of up to ~2
681 radians.
682 The default standard parallels are 33 and 45 to match the USGS state
684 1:500,000 base maps of the United States. At scales of 1:500,000 and
685 larger, discrepancies between the spherical and ellipsoidal projections
686 become important; use care with this projection on spheres.
687 OPTIONS
689 STANDARD POSITIONAL OPTIONS
691 s, std, standard, Standard (default (33,45))
692 The locations of the standard parallel(s) for the conic projection.
693 The transform reduces to the Mercator projection in the case where
694 the standard parallels are opposite and equidistant from the
695 equator, and in fact this is implemented by a call to t_mercator.
696 c, clip, Clip (default [-75,75])
698 Because the transform is conformal, the distant pole is displaced to
699 infinity. Many applications require a clipping boundary. The value
700 is in whatever angular unit you set with the standard 'unit' option.
701 For consistency with "t_mercator", clipping works the same way even
702 though in most cases only one pole needs it. Set this to 0 for no
703 clipping at all.
704 t_stereographic
706 $t = t_stereographic(<options>);
707 (Cartography) Stereographic projection (az.; conf.; persp.)
709 The stereographic projection is a true perspective (planar) projection
711 from a point on the spherical surface opposite the origin of the map.
712 OPTIONS
714 STANDARD POSITIONAL OPTIONS
716 c, clip, Clip (default 120)
717 This is the angular distance from the center to the edge of the
718 projected map. The default 120 degrees gives you most of the
719 opposite hemisphere but avoids the hugely distorted part near the
720 antipodes.
721 t_gnomonic
723 $t = t_gnomonic(<options>);
724 (Cartography) Gnomonic (focal-plane) projection (az.; persp.)
726 The gnomonic projection projects a hemisphere onto a tangent plane. It
728 is useful in cartography for the property that straight lines are great
729 circles; and it is useful in scientific imaging because it is the
730 projection generated by a simple optical system with a flat focal
731 plane.
732 OPTIONS
734 STANDARD POSITIONAL OPTIONS
736 c, clip, Clip (default 75)
737 This is the angular distance from the center to the edge of the
738 projected map. The default 75 degrees gives you most of the
739 hemisphere but avoids the hugely distorted part near the horizon.
740 t_az_eqd
742 $t = t_az_eqd(<options>);
743 (Cartography) Azimuthal equidistant projection (az.; equi.)
745 Basic azimuthal projection preserving length along radial lines from
747 the origin (meridians, in the original polar aspect). Hence, both
748 azimuth and distance are correct for journeys beginning at the origin.
749 Applied to the celestial sphere, this is the projection made by fisheye
751 lenses; it is also the projection into which "t_vertical" puts
752 perspective views.
753 The projected plane scale is normally taken to be planetary radii; this
755 is useful for cartographers but not so useful for scientific observers.
756 Setting the 't=>1' option causes the output scale to shift to camera
757 angular coordinates (the angular unit is determined by the standard
758 'Units' option; default is degrees).
759 OPTIONS
761 STANDARD POSITIONAL OPTIONS
763 c, clip, Clip (default 180 degrees)
764 The largest angle relative to the origin. Default is the whole
765 sphere.
766 t_az_eqa
768 $t = t_az_eqa(<options>);
769 (Cartography) Azimuthal equal-area projection (az.; auth.)
771 OPTIONS
773 STANDARD POSITIONAL OPTIONS
775 c, clip, Clip (default 180 degrees)
776 The largest angle relative to the origin. Default is the whole
777 sphere.
778 t_aitoff
780 "t_aitoff" in an alias for "t_hammer"
781 t_hammer
783 (Cartography) Hammer/Aitoff elliptical projection (az.; auth.)
784 The Hammer/Aitoff projection is often used to display the Celestial
786 sphere. It is mathematically related to the Lambert Azimuthal Equal-
787 Area projection ("t_az_eqa"), and maps the sphere to an ellipse of unit
788 eccentricity, with vertical radius sqrt(2) and horizontal radius of 2
789 sqrt(2).
790 OPTIONS
792 STANDARD POSITIONAL OPTIONS
794 t_zenithal
796 Vertical projections are also called "zenithal", and "t_zenithal" is an
797 alias for "t_vertical".
798 t_vertical
800 $t = t_vertical(<options>);
801 (Cartography) Vertical perspective projection (az.; persp.)
803 Vertical perspective projection is a generalization of gnomonic and
805 stereographic projection, and a special case of perspective projection.
806 It is a projection from the sphere onto a tangent plane from a point at
807 the camera location.
808 OPTIONS
810 STANDARD POSITIONAL OPTIONS
812 m, mask, Mask, h, hemisphere, Hemisphere [default 'near']
813 The hemisphere to keep in the projection (see
814 PDL::Transform::Cartography).
815 r0, R0, radius, d, dist, distance [default 2.0]
817 The altitude of the focal plane above the center of the sphere. The
818 default places the point of view one radius above the surface.
819 t, telescope, Telescope, cam, Camera (default '')
821 If this is set, then the central scale is in telescope or camera
822 angular units rather than in planetary radii. The angular units are
823 parsed as with the normal 'u' option for the lon/lat specification.
824 If you specify a non-string value (such as 1) then you get
825 telescope-frame radians, suitable for working on with other
826 transformations.
827 f, fish, fisheye (default '')
829 If this is set then the output is in azimuthal equidistant
830 coordinates instead of in tangent-plane coordinates. This is a
831 convenience function for '(t_az_eqd) x !(t_gnomonic) x
832 (t_vertical)'.
833 t_perspective
835 $t = t_perspective(<options>);
836 (Cartography) Arbitrary perspective projection
838 Perspective projection onto a focal plane from an arbitrary location
840 within or without the sphere, with an arbitrary central look direction,
841 and with correction for magnification within the optical system.
842 In the forward direction, t_perspective generates perspective views of
844 a sphere given (lon/lat) mapping or vector information. In the reverse
845 direction, t_perspective produces (lon/lat) maps from aerial or distant
846 photographs of spherical objects.
847 Viewpoints outside the sphere treat the sphere as opaque by default,
849 though you can use the 'm' option to specify either the near or far
850 surface (relative to the origin). Viewpoints below the surface treat
851 the sphere as transparent and undergo a mirror reversal for consistency
852 with projections that are special cases of the perspective projection
853 (e.g. t_gnomonic for r0=0 or t_stereographic for r0=-1).
854 Magnification correction handles the extra edge distortion due to
856 higher angles between the focal plane and focused rays within the
857 optical system of your camera. If you do not happen to know the
858 magnification of your camera, a simple rule of thumb is that the
859 magnification of a reflective telescope is roughly its focal length
860 (plate scale) divided by its physical length; and the magnification of
861 a compound refractive telescope is roughly twice its physical length
862 divided by its focal length. Simple optical systems with a single
863 optic have magnification = 1. Fisheye lenses have magnification < 1.
864 This transformation was derived by direct geometrical calculation
866 rather than being translated from Voxland & Snyder.
867 OPTIONS
869 STANDARD POSITIONAL OPTIONS
871 As always, the 'origin' field specifies the sub-camera point on the
872 sphere.
873 The 'roll' option is the roll angle about the sub-camera point, for
875 consistency with the other projectons.
876 p, ptg, pointing, Pointing (default (0,0,0))
878 The pointing direction, in (horiz. offset, vert. offset, roll) of
879 the camera relative to the center of the sphere. This is a
880 spherical coordinate system with the origin pointing directly at the
881 sphere and the pole pointing north in the pre-rolled coordinate
882 system set by the standard origin. It's most useful for space-based
883 images taken some distance from the body in question (e.g. images of
884 other planets or the Sun).
885 Be careful not to confuse 'p' (pointing) with 'P' (P angle, a
887 standard synonym for roll).
888 c, cam, camera, Camera (default undef)
890 Alternate way of specifying the camera pointing, using a spherical
891 coordinate system with poles at the zenith (positive) and nadir
892 (negative) -- this is useful for aerial photographs and such, where
893 the point of view is near the surface of the sphere. You specify
894 (azimuth from N, altitude from horizontal, roll from vertical=up).
895 If you specify pointing by this method, it overrides the 'pointing'
896 option, above. This coordinate system is most useful for aerial
897 photography or low-orbit work, where the nadir is not necessarily
898 the most interesting part of the scene.
899 r0, R0, radius, d, dist, distance [default 2.0]
901 The altitude of the point of view above the center of the sphere.
902 The default places the point of view 1 radius aboove the surface.
903 Do not confuse this with 'r', the standard origin roll angle!
904 Setting r0 < 1 gives a viewpoint inside the sphere. In that case,
905 the images are mirror-reversed to preserve the chiralty of the
906 perspective. Setting r0=0 gives gnomonic projections; setting r0=-1
907 gives stereographic projections. Setting r0 < -1 gives strange
908 results.
909 iu, im_unit, image_unit, Image_Unit (default 'degrees')
911 This is the angular units in which the viewing camera is calibrated
912 at the center of the image.
913 mag, magnification, Magnification (default 1.0)
915 This is the magnification factor applied to the optics -- it affects
916 the amount of tangent-plane distortion within the telescope. 1.0
917 yields the view from a simple optical system; higher values are
918 telescopic, while lower values are wide-angle (fisheye). Higher
919 magnification leads to higher angles within the optical system, and
920 more tangent-plane distortion at the edges of the image. The
921 magnification is applied to the incident angles themselves, rather
922 than to their tangents (simple two-element telescopes magnify
923 tan(theta) rather than theta itself); this is appropriate because
924 wide-field optics more often conform to the equidistant azimuthal
925 approximation than to the tangent plane approximation. If you need
926 more detailed control of the relationship between incident angle and
927 focal-plane position, use mag=1.0 and compose the transform with
928 something else to tweak the angles.
929 m, mask, Mask, h, hemisphere, Hemisphere [default 'near']
931 'hemisphere' is by analogy to other cartography methods although the
932 two regions to be selected are not really hemispheres.
933 f, fov, field_of_view, Field_Of_View [default 60 degrees]
935 The field of view of the telescope -- sets the crop radius on the
936 focal plane. If you pass in a scalar, you get a circular crop. If
937 you pass in a 2-element list ref, you get a rectilinear crop, with
938 the horizontal 'radius' and vertical 'radius' set separately.
939 EXAMPLES
941 Model a camera looking at the Sun through a 10x telescope from Earth
943 (~230 solar radii from the Sun), with an 0.5 degree field of view and a
944 solar P (roll) angle of 30 degrees, in February (sub-Earth solar
945 latitude is 7 degrees south). Convert a solar FITS image taken with
946 that camera to a FITS lon/lat map of the Sun with 20 pixels/degree
947 latitude:
948 # Define map output header (no need if you don't want a FITS output map)
950 $maphdr = {NAXIS1=>7200,NAXIS2=>3600, # Size of image
951 CTYPE1=>longitude,CTYPE2=>latitude, # Type of axes
952 CUNIT1=>deg,CUNIT2=>deg, # Unit of axes
953 CDELT1=>0.05,CDELT2=>0.05, # Scale of axes
954 CRPIX1=>3601,CRPIX2=>1801, # Center of map
955 CRVAL1=>0,CRVAL2=>0 # (lon,lat) of center
956 };
957 # Set up the perspective transformation, and apply it.
959 $t = t_perspective(r0=>229,fov=>0.5,mag=>10,P=>30,B=>-7);
960 $map = $im->map( $t , $maphdr );
961 Draw an aerial-view map of the Chesapeake Bay, as seen from a sounding
963 rocket at an altitude of 100km, looking NNE from ~200km south of
964 Washington (the radius of Earth is 6378 km; Washington D.C. is at
965 roughly 77W,38N). Superimpose a linear coastline map on a photographic
966 map.
967 $x = graticule(1,0.1)->glue(1,earth_coast());
969 $t = t_perspective(r0=>6478/6378.0,fov=>60,cam=>[22.5,-20],o=>[-77,36])
970 $w = pgwin(size=>[10,6],J=>1);
971 $w->fits_imag(earth_image()->map($t,[800,500],{m=>linear}));
972 $w->hold;
973 $w->lines($x->apply($t),{xt=>'Degrees',yt=>'Degrees'});
974 $w->release;
975 Model a 5x telescope looking at Betelgeuse with a 10 degree field of
977 view (since the telescope is looking at the Celestial sphere, r is 0
978 and this is just an expensive modified-gnomonic projection).
979 $t = t_perspective(r0=>0,fov=>10,mag=>5,o=>[88.79,7.41])
981perl v5.36.0 2023-01-20 Cartography(3)