PDL::Transform::Cartography - Useful cartographic projections |
PDL::Transform::Cartography - Useful cartographic projections
# make a Mercator map of Earth use PDL::Transform::Cartography; $a = earth_coast(); $a = graticule(10,2)->glue(1,$a); $t = t_mercator; $w = pgwin(xs); $w->lines($t->apply($a)->clean_lines());
PDL::Transform::Cartography includes a variety of useful cartographic and observing projections (mappings of the surface of a sphere), including reprojected observer coordinates. See the PDL::Transform manpage for more information about image transforms in general.
Cartographic transformations are used for projecting not just terrestrial maps, but also any nearly spherical surface including the Sun, the Celestial sphere, various moons and planets, distant stars, etc. They also are useful for interpreting scientific images, which are themselves generally projections of a sphere onto a flat focal plane (e.g. the t_gnomonic projection).
Unless otherwise noted, all the transformations in this file convert from (theta,phi) coordinates on the unit sphere (e.g. (lon,lat) on a planet or (RA,dec) on the celestial sphere) into some sort of projected coordinates, and have inverse transformations that convert back to (theta,phi). This is equivalent to working from the equidistant cylindrical (or ``plate caree'') projection, if you are a cartography wonk.
The projected coordinates are generally in units of body radii (radians), so that multiplying the output by the scale of the map yields physical units that are correct wherever the scale is correct for that projection. For example, areas should be correct everywhere in the authalic projections; and linear scales are correct along meridians in the equidistant projections and along the standard parallels in all the projections.
The transformations that are authalic (equal-area), conformal (equal-angle), azimuthal (circularly symmetric), or perspective (true perspective on a focal plane from some viewpoint) are marked. The first two categories are mutually exclusive for all but the ``unit sphere'' 3-D projection.
Extra dimensions tacked on to each point to be transformed are, in general, ignored. That is so that you can add on an extra index to keep track of pen color. For example, earth_coast returns a 3x<n> piddle containing (lon, lat, pen) at each list location. Transforming the vector list retains the pen value as the first index after the dimensional directions.
Unless otherwise noted, the transformations and miscellaneous information in this section are taken from Snyder & Voxland 1989: ``An Album of Map Projections'', US Geological Survey Professional Paper 1453, US Printing Office (Denver); and from Snyder 1987: ``Map Projections - A Working Manual'', US Geological Survey Professional Paper 1395, US Printing Office (Denver, USA). You can obtain your own copy of both by contacting the U.S. Geological Survey, Federal Center, Box 25425, Denver, CO 80225 USA.
The mathematics of cartography have a long history, and the details are far trickier than the broad overview. For terrestrial (and, in general, planetary) cartography, the best reference datum is not a sphere but an oblate ellipsoid due to centrifugal force from the planet's rotation. Furthermore, because all rocky planets, including Earth, have randomly placed mass concentrations that affect the gravitational field, the reference gravitational isosurface (sea level on Earth) is even more complex than an ellipsoid and, in general, different ellipsoids have been used for different locations at the same time and for the same location at different times.
The transformations in this package use a spherical datum and hence
include global distortion at about the 0.5% level for terrestrial maps
(Earth's oblateness is ~1/300). This is roughly equal to the
dimensional precision of physical maps printed on paper (due to
stretching and warping of the paper) but is significant at larger
scales (e.g. for regional maps). If you need more precision than
that, you will want to implement and use the ellipsoidal
transformations from Snyder 1987 or another reference work on geodesy.
A good name for that package would be ...::Cartography::Geodetic
.
Cartographic transformations are useful for interpretation of scientific images, as all cameras produce projections of the celestial sphere onto the focal plane of the camera. A simple (single-element) optical system with a planar focal plane generates gnomonic images -- that is to say, gnomonic projections of a portion of the celestial sphere near the paraxial direction. This is the projection that most consumer grade cameras produce.
Magnification in an optical system changes the angle of incidence of the rays on the focal plane for a given angle of incidence at the aperture. For example, a 10x telescope with a 2 degree field of view exhibits the same gnomonic distortion as a simple optical system with a 20 degree field of view. Wide-angle optics typically have magnification less than 1 ('fisheye lenses'), reducing the gnomonic distortion considerably but introducing ``equidistant azimuthal'' distortion -- there's no such thing as a free lunch!
Because many solar-system objects are spherical, PDL::Transform::Cartography includes perspective projections for producing maps of spherical bodies from perspective views. Those projections are ``t_vertical'' and ``t_perspective''. They map between (lat,lon) on the spherical body and planar projected coordinates at the viewpoint. ``t_vertical'' is the vertical perspective projection given by Snyder, but ``t_perspective'' is a fully general perspective projection that also handles magnification correction.
Oblique projections rotate the sphere (and graticule) to an arbitrary angle before generating the projection; transverse projections rotate the sphere exactly 90 degrees before generating the projection.
Most of the projections accept the following standard options, useful for making transverse and oblique projection maps.
Draw a Mercator map of the world on-screen:
$w = pgwin(xs); $w->lines(earth_coast->apply(t_mercator)->clean_lines);
Here, earth_coast()
returns a 3xn piddle containing (lon, lat, pen)
values for the included world coastal outline; t_mercator
converts
the values to projected Mercator coordinates, and clean_lines
breaks
lines that cross the 180th meridian.
Draw a Mercator map of the world, with lon/lat at 10 degree intervals:
$w = pgwin(xs) $a = earth_coast()->glue(1,graticule(10,1)); $w->lines($a->apply(t_mercator)->clean_lines);
This works just the same as the first example, except that a map graticule has been applied with interline spacing of 10 degrees lon/lat and inter-vertex spacing of 1 degree (so that each meridian contains 181 points, and each parallel contains 361 points).
Currently angular conversions are rather simpleminded. A list of common conversions is present in the main constructor, which inserts a conversion constant to radians into the {params} field of the new transform. Something like Math::Convert::Units should be used instead to generate the conversion constant.
A cleaner higher-level interface is probably needed (see the examples); for example, earth_coast could return a graticule if asked, instead of needing one to be glued on.
The class structure is somewhat messy because of the varying needs of the different transformations. PDL::Transform::Cartography is a base class that interprets the origin options and sets up the basic machinery of the Transform. The conic projections have their own subclass, PDL::Transform::Conic, that interprets the standard parallels. Since the cylindrical and azimuthal projections are pretty simple, they are not subclassed.
The perl 5.6.1 compiler is quite slow at adding new classes to the structure, so it does not makes sense to subclass new transformations merely for the sake of pedantry.
Copyright 2002, Craig DeForest (deforest@boulder.swri.edu) This module may be modified and distributed under the same terms as PDL itself. The module comes with NO WARRANTY.
The included digital world map is derived from the 1987 CIA World Map, translated to ASCII in 1988 by Joe Dellinger (geojoe@freeusp.org) and simplified in 1995 by Kirk Johnson (tuna@indra.com) for the program XEarth. The map comes with NO WARRANTY. An ASCII version of the map, and a sample PDL function to read it, may be found in the Demos subdirectory of the PDL source distribution.
The module exports both transform constructors ('t_<foo>') and some auxiliary functions (no leading 't_').
$lonlatp = graticule(<grid-spacing>,<line-segment-size>);
(Cartography) PDL constructor - generate a lat/lon grid.
Returns a grid of meridians and parallels as a list of vectors suitable for sending to PDL::Graphics::PGPLOT::Window::lines for plotting. The grid is in degrees in (theta, phi) coordinates -- this is (E lon, N lat) for terrestrial grids or (RA, dec) for celestial ones. You must then transform the graticule in the same way that you transform the map.
You can attach the graticule to a vector map using the syntax:
$out = graticule(10,2)->glue(1,$map);
In array context you get back a 2-element list containing a piddle of the (theta,phi) pairs and a piddle of the pen values (1 or 0) suitable for calling PDL::Graphics::PGPLOT::Window::lines. In scalar context the two elements are combined into a single piddle.
The pen values associated with the graticule are negative, which will cause PDL::Graphics::PGPLOT::Window::lines to plot them as hairlines.
$a = earth_coast()
(Cartography) PDL constructor - coastline map of Earth
Returns a vector coastline map based on the 1987 CIA World Coastline database (see author information). The vector coastline data are in plate caree format so they can be converted to other projections via the apply method and cartographic transforms, and are suitable for plotting with the lines method in the PGPLOT output library: the first dimension is (X,Y,pen) with breaks having a pen value of 0 and hairlines having negative pen values. The second dimension threads over all the points in the data set.
The vector map includes lines that pass through the antipodean meridian, so if you want to plot it without reprojecting, you should run it through clean_lines first:
$w = pgwin(); $w->lines(earth_coast->clean_lines); # plot plate caree map of world $w->lines(earth_coast->apply(t_gnomonic))# plot gnomonic map of world
earth_coast
is just a quick-and-dirty way of loading the file
``earth_coast.vec.fits'' that is part of the normal installation tree.
$rgb = earth_image()
(Cartography) PDL constructor - RGB pixel map of Earth
Returns an RGB image of Earth based on data from the MODIS instrument on the NASA EOS/Terra satellite. (You can get a full-resolution image from http://earthobservatory.nasa.gov/Newsroom/BlueMarble/). The image is a plate caree map, so you can convert it to other projections via the map method and cartographic transforms.
This is just a quick-and-dirty way of loading the earth-image files that are distributed along with PDL.
$a = clean_lines(t_mercator->apply(scalar(earth_coast()))); $a = clean_lines($line_pen, [threshold]); $a = $lines->clean_lines;
(Cartography) PDL method - remove projection irregularities
clean_lines
massages vector data to remove jumps due to singularities
in the transform.
In the first (scalar) form, $line_pen
contains both (X,Y) points and pen
values suitable to be fed to lines;
in the second (list) form, $lines
contains the (X,Y) points and $pen
contains the pen values.
clean_lines
assumes that all the outline polylines are local --
that is to say, there are no large jumps. Any jumps larger than a
threshold size are broken by setting the appropriate pen values to 0.
The threshold
parameter sets the relative size of the largest jump, relative
to the map range (as determined by a min/max operation). The default size is
0.1.
NOTES
This almost never catches stuff near the apex of cylindrical maps, because the anomalous vectors get arbitrarily small. This could be improved somewhat by looking at individual runs of the pen and using a relative length scale that is calibrated to the rest of each run. it is probably not worth the computational overhead.
$t = t_unit_sphere(<options>);
(Cartography) 3-D globe projection (conformal; authalic)
This is similar to the inverse of t_spherical, but the inverse transform projects 3-D coordinates onto the unit sphere, yielding only a 2-D (lon/lat) output. Similarly, the forward transform deprojects 2-D (lon/lat) coordinates onto the surface of a unit sphere.
The cartesian system has its Z axis pointing through the pole of the (lon,lat) system, and its X axis pointing through the equator at the prime meridian.
Unit sphere mapping is unusual in that it is both conformal and authalic. That is possible because it properly embeds the sphere in 3-space, as a notional globe.
This is handy as an intermediate step in lots of transforms, as Cartesian 3-space is cleaner to work with than spherical 2-space.
Higher dimensional indices are preserved, so that ``rider'' indices (such as pen value) are propagated.
There is no oblique transform for t_unit_sphere, largely because it's so easy to rotate the output using t_linear once it's out into Cartesian space. In fact, the other projections implement oblique transforms by wrapping t_linear with t_unit_sphere.
OPTIONS:
$t = t_rot_sphere({origin=>[<theta>,<phi>],roll=>[<roll>]});
(Cartography) Generate oblique projections
You feed in the origin in (theta,phi) and a roll angle, and you get back out (theta', phi') coordinates. This is useful for making oblique or transverse projections: just compose t_rot_sphere with your favorite projection and you get an oblique one.
Most of the projections automagically compose themselves with t_rot_sphere if you feed in an origin or roll angle.
t_rot_sphere converts the base plate caree projection (straight lon, straight lat) to a Cassini projection.
OPTIONS
$t = t_orthographic(<options>);
(Cartography) Ortho. projection (azimuthal; perspective)
This is a perspective view as seen from infinite distance. You can specify the sub-viewer point in (lon,lat) coordinates, and a rotation angle of the map CW from (north=up). This is equivalent to specify viewer roll angle CCW from (north=up).
t_orthographic is a convenience interface to t_unit_sphere -- it is implemented as a composition of a t_unit_sphere call, a rotation, and a slice.
[*] In the default case where the near hemisphere is mapped, the inverse exists. There is no single inverse for the whole-sphere case, so the inverse transform superimposes everything on a single hemisphere. If you want an invertible 3-D transform, you want t_unit_sphere.
OPTIONS
NOTES
Alone of the various projections, this one does not use t_rot_sphere to handle the standard options, because the cartesian coordinates of the rotated sphere are already correctly projected -- t_rot_sphere would put them back into (theta', phi') coordinates.
$t = t_caree(<options>);
(Cartography) Plate Caree projection (cylindrical; equidistant)
This is the simple Plate Caree projection -- also called a ``lat/lon plot''. The horizontal axis is theta; the vertical axis is phi. This is a no-op if the angular unit is radians; it is a simple scale otherwise.
OPTIONS
$t = t_mercator(<options>);
(Cartography) Mercator projection (cylindrical; conformal)
This is perhaps the most famous of all map projections: meridians are mapped to parallel vertical lines and parallels are unevenly spaced horizontal lines. The poles are shifted to +/- infinity. The output values are in units of globe-radii for easy conversion to kilometers; hence the horizontal extent is -pi to pi.
You can get oblique Mercator projections by specifying the origin
or
roll
options; this is implemented via t_rot_sphere.
OPTIONS
$t = t_sin_lat(<options>);
(Cartography) Cyl. equal-area projection (cyl.; authalic)
This projection is commonly used in solar Carrington plots; but not much for terrestrial mapping.
OPTIONS
$t = t_sinusoidal(<options>);
(Cartography) Sinusoidal projection (authalic)
Sinusoidal projection preserves the latitude scale but scales longitude according to sin(lat); in this respect it is the companion to t_sin_lat, which is also authalic but preserves the longitude scale instead.
OPTIONS
$t = t_conic(<options>)
(Cartography) Simple conic projection (conic; equidistant)
This is the simplest conic projection, with parallels mapped to equidistant concentric circles. It is neither authalic nor conformal. This transformation is also referred to as the ``Modified Transverse Mercator'' projection in several maps of Alaska published by the USGS; and the American State of New Mexico re-invented the projection in 1936 for an official map of that State.
OPTIONS
parallel(s)
(where the cone intersects
the surface of the sphere). If you specify only one then the other is
taken to be the nearest pole. If you specify both of them to be one
pole then you get an equidistant azimuthal map. If you specify both
of them to be opposite and equidistant from the equator you get a
Plate Caree projection.
$t = t_albers(<options>)
(Cartography) Albers conic projection (conic; authalic)
This is the standard projection used by the US Geological Survey for sectionals of the 50 contiguous United States of America.
The projection reduces to the Lambert equal-area conic (infrequently used and not to be confused with the Lambert conformal conic, t_lambert!) if the pole is used as one of the two standard parallels.
Notionally, this is a conic projection onto a cone that intersects the sphere at the two standard parallels; it works best when the two parallels straddle the region of interest.
OPTIONS
The default values for the standard parallels are those chosen by Adams for maps of the lower 48 US states: (29.5,45.5). The USGS recommends (55,65) for maps of Alaska and (8,18) for maps of Hawaii -- these latter are chosen to also include the Canal Zone and Philippine Islands farther south, which is why both of those parallels are south of the Hawaiian islands.
The transformation reduces to the cylindrical equal-area (sin-lat) transformation in the case where the standard parallels are opposite and equidistant from the equator, and in fact this is implemented by a call to t_sin_lat.
$t = t_lambert(<options>);
(Cartography) Lambert conic projection (conic; conformal)
Lambert conformal conic projection is widely used in aeronautical charts and state base maps published by the USA's FAA and USGS. It's especially useful for mid-latitude charts. In particular, straight lines approximate (but are not exactly) great circle routes of up to ~2 radians.
The default standard parallels are 33 and 45 to match the USGS state 1:500,000 base maps of the United States. At scales of 1:500,000 and larger, discrepancies between the spherical and ellipsoidal projections become important; use care with this projection on spheres.
OPTIONS
parallel(s)
for the conic projection.
The transform reduces to the Mercator projection in the case where the
standard parallels are opposite and equidistant from the equator, and
in fact this is implemented by a call to t_mercator.
$t = t_stereographic(<options>);
(Cartography) Stereographic projection (az.; conf.; persp.)
The stereographic projection is a true perspective (planar) projection from a point on the spherical surface opposite the origin of the map.
OPTIONS
$t = t_gnomonic(<options>);
(Cartography) Gnomonic (focal-plane) projection (az.; persp.)
The gnomonic projection projects a hemisphere onto a tangent plane. It is useful in cartography for the property that straight lines are great circles; and it is useful in scientific imaging because it is the projection generated by a simple optical system with a flat focal plane.
OPTIONS
$t = t_az_eqd(<options>);
(Cartography) Azimuthal equidistant projection (az.; equi.)
Basic azimuthal projection preserving length along radial lines from the origin (meridians, in the original polar aspect). Hence, both azimuth and distance are correct for journeys beginning at the origin.
Applied to the celestial sphere, this is the projection made by
fisheye lenses; it is also the projection into which t_vertical
puts perspective views.
The projected plane scale is normally taken to be planetary radii; this is useful for cartographers but not so useful for scientific observers. Setting the 't=>1' option causes the output scale to shift to camera angular coordinates (the angular unit is determined by the standard 'Units' option; default is degrees).
OPTIONS
$t = t_az_eqa(<options>);
(Cartography) Azimuthal equal-area projection (az.; auth.)
OPTIONS
(Cartography) Hammer/Aitoff elliptical projection (az.; auth.)
The Hammer/Aitoff projection is often used to display the Celestial
sphere. It is mathematically related to the Lambert Azimuthal Equal-Area
projection (t_az_eqa), and maps the sphere to an ellipse of unit
eccentricity, with vertical radius sqrt(2)
and horizontal radius of
2 sqrt(2).
OPTIONS
$t = t_vertical(<options>);
(Cartography) Vertical perspective projection (az.; persp.)
Vertical perspective projection is a generalization of gnomonic and stereographic projection, and a special case of perspective projection. It is a projection from the sphere onto a focal plane at the camera location.
OPTIONS
$t = t_perspective(<options>);
(Cartography) Arbitrary perspective projection
Perspective projection onto a focal plane from an arbitrary location within or without the sphere, with an arbitary central look direction, and with correction for magnification within the optical system.
In the forward direction, t_perspective generates perspective views of a sphere given (lon/lat) mapping or vector information. In the reverse direction, t_perspective produces (lon/lat) maps from aerial or distant photographs of spherical objects.
Viewpoints outside the sphere treat the sphere as opaque by default, though you can use the 'm' option to specify either the near or far surface (relative to the origin). Viewpoints below the surface treat the sphere as transparent and undergo a mirror reversal for consistency with projections that are special cases of the perspective projection (e.g. t_gnomonic for r0=0 or t_stereographic for r0=-1).
Magnification correction handles the extra edge distortion due to higher angles between the focal plane and focused rays within the optical system of your camera. If you do not happen to know the magnification of your camera, a simple rule of thumb is that the magnification of a reflective telescope is roughly its focal length (plate scale) divided by its physical length; and the magnification of a compound refractive telescope is roughly twice its physical length divided by its focal length. Simple optical sytems with a single optic have magnification = 1. Fisheye lenses have magnification < 1.
This transformation was derived by direct geometrical calculation rather than being translated from Voxland & Snyder.
OPTIONS
The 'roll' option is the roll angle about the sub-camera point, for consistency with the other projectons.
Be careful not to confuse 'p' (pointing) with 'P' (P angle, a standard synonym for roll).
tan(theta)
rather
than theta itself); this is appropriate because wide-field optics more
often conform to the equidistant azimuthal approximation than to the
tangent plane approximation. If you need more detailed control of
the relationship between incident angle and focal-plane position,
use mag=1.0 and compose the transform with something else to tweak the
angles.
EXAMPLES
Model a camera looking at the Sun through a 10x telescope from Earth (~230 solar radii from the Sun), with an 0.5 degree field of view and a solar P (roll) angle of 30 degrees, in February (sub-Earth solar latitude is 7 degrees south). Convert a solar FITS image taken with that camera to a FITS lon/lat map of the Sun with 20 pixels/degree latitude:
# Define map output header (no need if you don't want a FITS output map) $maphdr = {NAXIS1=>7200,NAXIS2=>3600, # Size of image CTYPE1=>longitude,CTYPE2=>latitude, # Type of axes CUNIT1=>deg,CUNIT2=>deg, # Unit of axes CDELT1=>0.05,CDELT2=>0.05, # Scale of axes CRPIX1=>3601,CRPIX2=>1801, # Center of map CRVAL1=>0,CRVAL2=>0 # (lon,lat) of center };
# Set up the perspective transformation, and apply it. $t = t_perspective(r0=>229,fov=>0.5,mag=>10,P=>30,B=>-7); $map = $im->map( $t , $maphdr );
Draw an aerial-view map of the Chesapeake Bay, as seen from a sounding rocket at an altitude of 100km, looking NNE from ~200km south of Washington (the radius of Earth is 6378 km; Washington D.C. is at roughly 77W,38N). Superimpose a linear coastline map on a photographic map.
$a = graticule(1,0.1)->glue(1,earth_coast()); $t = t_perspective(r0=>6478/6378.0,fov=>60,cam=>[22.5,-20],o=>[-77,36]) $w = pgwin(size=>[10,6],J=>1); $w->fits_imag(earth_image()->map($t,[800,500],{m=>linear})); $w->hold; $w->lines($a->apply($t),{xt=>'Degrees',yt=>'Degrees'}); $w->release;
Model a 5x telescope looking at Betelgeuse with a 10 degree field of view (since the telescope is looking at the Celestial sphere, r is 0 and this is just an expensive modified-gnomonic projection).
$t = t_perspective(r0=>0,fov=>10,mag=>5,o=>[88.79,7.41])
PDL::Transform::Cartography - Useful cartographic projections |