Summary

###########################################################################

	File:     cartog.icn

	Subject:  Procedures for cartographic projection 

	Authors:  Gregg M. Townsend and William S. Evans

	Date:     May 24, 2000

###########################################################################

   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 specifed, 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.

###########################################################################

	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.

###########################################################################

	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

###########################################################################

   Links: geodat, io, lu, numbers, strings

###########################################################################

This file is part of the (main) package.

compose(a)

: define composite projection

ctg_comp_inv(p)

ctg_comp_proj(p, L)

ctg_iutm_inv(p)

ctg_iutm_proj(p, L)

ctg_ppt_inv(p)

ctg_ppt_proj(p, L)

ctg_rect_inv(p)

ctg_rect_proj(p, L)

ctg_utm_inv(p)

ctg_utm_proj(p, L)

invp(p)

: return inversion of projection

pptrans(L1, L2)

: define planar projective transformation

project(p, L)

: project a list of coordinates

rectp(x1, y1, x2, y2, xm, ym)

: define rectangular projection

utm(a, f)

: define UTM projection

ctg_comp(proj, inv, projList)

 composition of two projections

ctg_ppt(proj, inv, org, tgt, h11, h12, h13, h21, h22, h23, h31, h32, h33)

 planar projective transformation record

ctg_rect(proj, inv, xmul, ymul, xadd, yadd)

 rectangular projection record

ctg_utm(proj, inv, a, f, e, esq, epsq, c0, c2, c4, c6, c8)

 UTM projection record