procedure project: project a list of coordinates procedure invp: return inversion of projection procedure sbsize: calculate scalebar size procedure rectp: define rectangular projection procedure pptrans: define planar projective transformation procedure utm: define UTM projection procedure compose: define composite projection
link cartog
February 19, 2003; Gregg M. Townsend and William S. Evans
This file is in the public domain.
These procedures project geographic coordinates.
rectp(x1, y1, x2, y2, xm, ym) defines a rectangular projection.
pptrans(L1, L2) defines a planar projective transformation.
utm(a, f) defines a latitude/longitude to UTM projection.
project(p, L) projects a list of coordinates.
invp(p) returns the inverse of projection p.
compose(p1, p2, ...) creates a composite projection.
____________________________________________________________
rectp(x1, y1, x2, y2, xm, ym) returns a rectangular projection
in which the point (x1, y1) maps to (x2, y2). If xm is specified,
distances in the projected coordinate system are scaled by xm. If
ym is also specified, xm scales x values while ym scales y values.
____________________________________________________________
pptrans(L1, L2) returns a planar projective transform that maps
the four points in L1 to the four points in L2. Each of the two
lists contains 8 coordinates: [x1, y1, x2, y2, x3, y3, x4, y4].
____________________________________________________________
utm(a, f) returns a projection from latitude and longitude to
Universal Transverse Mercator (UTM) representation. The reference
ellipsoid is specified by a, the equatorial radius in metres, and f,
the flattening. Alternatively, f can be omitted with a specifying
a string, such as "Clarke66"; if a is also omitted, "WGS84" is used.
See ellipsoid() in geodat.icn for the list of possible strings.
The input list contains signed numeric values: longitude and
latitude, in degrees, in that order (x before y). The output list
contains triples: an integer zone number followed by real-valued
UTM x and y distances in metres. No "false easting" is applied.
UTM conversions are valid between latitudes 72S and 84N, except
for those portions of Norway where the UTM grid is irregular.
____________________________________________________________
project(p, L) applies a projection, reading a list of coordinates
and returning a new list of transformed coordinates.
____________________________________________________________
invp(p) returns the inverse of projection p, or fails if no
inverse projection is available.
____________________________________________________________
compose(p1, p2, ..., pn) returns the projection that is the
composition of the projections p1, p2, ..., pn. The composition
applies pn first.
____________________________________________________________
sbsize(p, x, y, u, maxl) calculates a scale bar size for use with
projection p at input coordinates (x, y). Given u, the size of
an unprojected convenient unit (meter, foot, mile, etc.) at (x, y),
sbsize() returns the maximum "round number" N such that
-- N is of the form i * 10 ^ k for i in {1,2,3,4,5}
-- the projected length of the segment (x, y, x + N * u, y)
does not exceed maxl
____________________________________________________________
UTM conversion algorithms are based on:
Map Projections: A Working Manual
John P. Snyder
U.S. Geological Survey Professional Paper 1395
Washington: Superintendent of Documents, 1987
Planar projective transformation calculations come from:
Computing Plane Projective Transformations (Method 1)
Andrew Zisserman, Robotics Research Group, Oxford
in CVOnline (R. Fisher, ed.), found 22 February 2000 at:
http://www.dai.ed.ac.uk/CVonline/LOCAL_COPIES/EPSRC_SSAZ/node11.html