############################################################################
#
# 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
#
############################################################################
link geodat
link io
link lu
link numbers
link strings
# Procedures and globals named with a "ctg_" prefix are
# not intended for access outside this file.
global ctg_eps_ptab # table of [axis, flatng], keyed by eps name
#################### General Projection Support ####################
# project(p, L) projects a list of coordinates, returning a new list.
procedure project(p, L) #: project a list of coordinates
return p.proj(p, L)
end
# invp(p) returns the inverse of projection p.
procedure invp(p) #: return inversion of projection
return (\p.inv)(p)
end
#################### Rectangular Projection ####################
record ctg_rect( # rectangular projection record
proj, # projection procedure
inv, # inversion procedure
xmul, # x multiplier
ymul, # y multiplier
xadd, # x additive factor
yadd # y additive factor
)
# rectp(x1, y1, x2, y2, xm, ym) -- define rectangular projection
procedure rectp(x1, y1, x2, y2, xm, ym) #: define rectangular projection
local p
/xm := 1.0
/ym := xm
p := ctg_rect()
p.proj := ctg_rect_proj
p.inv := ctg_rect_inv
p.xmul := real(xm)
p.ymul := real(ym)
p.xadd := x2 - x1 * xm
p.yadd := y2 - y1 * ym
return p
end
# ctg_rect_proj(p, L) -- project using rectangular projection
procedure ctg_rect_proj(p, L)
local i, a, xmul, ymul, xadd, yadd
a := list()
xmul := p.xmul
ymul := p.ymul
xadd := p.xadd
yadd := p.yadd
every i := 1 to *L by 2 do {
put(a, xmul * L[i] + xadd)
put(a, ymul * L[i+1] + yadd)
}
return a
end
# ctg_rect_inv(p) -- invert rectangular projection
procedure ctg_rect_inv(p)
local q
q := copy(p)
q.xmul := 1.0 / p.xmul
q.ymul := 1.0 / p.ymul
q.xadd := -p.xadd / p.xmul
q.yadd := -p.yadd / p.ymul
return q
end
################ Planar Projective Transformation ###############
record ctg_ppt( # planar projective transformation record
proj, # projection procedure
inv, # inversion procedure
org, # origin points
tgt, # target points
h11, h12, h13, # transformation matrix: (x' y' 1) = H (x y 1)
h21, h22, h23,
h31, h32, h33
)
# pptrans(L1, L2) -- define planar projective transformation
procedure pptrans(L1, L2) #: define planar projective transformation
local p, M, I, B
local x1, x2, x3, x4, y1, y2, y3, y4
local x1p, x2p, x3p, x4p, y1p, y2p, y3p, y4p
*L1 = 8 | runerr(205, L1)
*L2 = 8 | runerr(205, L2)
p := ctg_ppt()
p.proj := ctg_ppt_proj
p.inv := ctg_ppt_inv
p.org := copy(L1)
p.tgt := copy(L2)
B := copy(L1)
every (x1 | y1 | x2 | y2 | x3 | y3 | x4 | y4) := get(B)
B := copy(L2)
every (x1p | y1p | x2p | y2p | x3p | y3p | x4p | y4p) := get(B)
M := [
[ x1, y1, 1., 0., 0., 0., -x1p * x1, -x1p * y1],
[ 0., 0., 0., x1, y1, 1., -y1p * x1, -y1p * y1],
[ x2, y2, 1., 0., 0., 0., -x2p * x2, -x2p * y2],
[ 0., 0., 0., x2, y2, 1., -y2p * x2, -y2p * y2],
[ x3, y3, 1., 0., 0., 0., -x3p * x3, -x3p * y3],
[ 0., 0., 0., x3, y3, 1., -y3p * x3, -y3p * y3],
[ x4, y4, 1., 0., 0., 0., -x4p * x4, -x4p * y4],
[ 0., 0., 0., x4, y4, 1., -y4p * x4, -y4p * y4]
]
I := list(8)
B := copy(L2)
lu_decomp(M, I) | fail # if singular, fail
lu_back_sub(M, I, B)
every (p.h11 | p.h12 | p.h13 | p.h21 | p.h22 | p.h23 | p.h31 | p.h32) :=
get(B)
p.h33 := 1.0
return p
end
# ctg_ppt_proj(p, L) -- project using planar projective transformation
procedure ctg_ppt_proj(p, L)
local a, i, x, y, d, h11, h12, h13, h21, h22, h23, h31, h32, h33
h11 := p.h11
h12 := p.h12
h13 := p.h13
h21 := p.h21
h22 := p.h22
h23 := p.h23
h31 := p.h31
h32 := p.h32
h33 := p.h33
a := list()
every i := 1 to *L by 2 do {
x := L[i]
y := L[i+1]
d := h31 * x + h32 * y + h33
put(a, (h11 * x + h12 * y + h13) / d, (h21 * x + h22 * y + h23) / d)
}
return a
end
# ctg_ppt_inv(p, L) -- invert planar projective transformation
procedure ctg_ppt_inv(p)
return pptrans(p.tgt, p.org)
end
############### Universal Transverse Mercator Projection ###############
# UTM conversion parameters
$define k0 0.9996 # central meridian scaling factor for UTM
$define M0 0.0 # M0 = 0 because y origin is at phi=0
record ctg_utm( # UTM projection record
proj, # projection procedure
inv, # inversion procedure
a, # polar radius
f, # flattening
e, # eccentricity
esq, # eccentricity squared
epsq, # e prime squared
c0, c2, c4, c6, c8 # other conversion constants
)
# utm(a, f) -- define UTM projection
procedure utm(a, f) #: define UTM projection
local p, e, af
p := ctg_utm()
p.proj := ctg_utm_proj
p.inv := ctg_utm_inv
if /f then {
af := ellipsoid(a) | fail
a := af[1]
f := af[2]
}
p.a := a # p.a = equatorial radius
p.f := f # p.f = flattening
p.esq := 2 * f - f ^ 2 # p.esq = eccentricity squared
p.epsq := p.esq / (1 - p.esq)
p.e := sqrt(p.esq) # p.e = eccentricity
p.c0 := p.a * (1 - (p.e^2) / 4 - 3 * (p.e^4) / 64 - 5 * (p.e^6) / 256)
p.c2 := p.a * (3 * (p.e^2) / 8 + 3 * (p.e^4) / 32 + 45 * (p.e^6) / 1024)
p.c4 := p.a * (15 * (p.e^4) / 256 + 45 * (p.e^6) / 1024)
p.c6 := p.a * (35 * (p.e^6) / 3072)
return p
end
# ctg_utm_proj(p, L) -- project using UTM projection (Snyder, p61)
procedure ctg_utm_proj(p, L)
local ulist, epsq, lat, lon, zone, phi, lambda, lamzero, cosphi
local i, N, T, C, A, M, x, u, y
ulist := list()
epsq := p.epsq
every i := 1 to *L by 2 do {
lon := numeric(L[i])
lat := numeric(L[i+1])
zone := (185 + integer(lon)) / 6
phi := dtor(lat) # latitude in radians
lambda := dtor(lon) # longitude in radians
lamzero := dtor(-183 + 6 * zone) # central meridian of zone
N := p.a / sqrt(1 - p.esq * sin(phi) ^ 2) # (8-12)
T := tan(phi) ^ 2 # (4-20)
cosphi := cos(phi)
C := epsq * cosphi ^ 2 # (8-13)
A := (lambda - lamzero) * cosphi # (8-15)
M := p.c0*phi - p.c2*sin(2.*phi) + p.c4*sin(4.*phi) - p.c6*sin(6.*phi)
x := k0 * N * (A + (1 - T + C) * A^3 / 6. +
(5. - 18. * T + T^2 + 72. * C - 58. * epsq) * A^5 / 120.)
u := A^2 / 2 + (5 - T + 9 * C + 4 * C^2) * A^4 / 24 +
(61. - 58. * T + T^2 + 600. * C - 330. * epsq) * A^6 / 720.
y := k0 * (M - M0 + N * tan(phi) * u)
put(ulist, zone, x, y)
}
return ulist
end
# ctg_utm_inv(p) -- invert UTM projection
procedure ctg_utm_inv(p)
local q, e, e1
q := copy(p)
q.proj := ctg_iutm_proj
q.inv := ctg_iutm_inv
e := q.e
e1 := (1 - sqrt(1 - e^2)) / (1 + sqrt(1 - e^2))
q.c0 := q.a * (1 - e^2 / 4. - 3. * e^4 / 64. - 5. * e^6 / 256.)
q.c2 := 3. * e1 / 2. - 27. * e1^3 / 32.
q.c4 := 21. * e1^2 / 16. - 55. * e1^4 / 32.
q.c6 := 151. * e1^3 / 96.
q.c8 := 1097. * e1^4 / 512.
return q
end
# ctg_iutm_proj(p, L) -- project using inverse UTM projection (Snyder, p63)
procedure ctg_iutm_proj(p, L)
local a, esq, epsq
local lllist, i, x, y, zone
local lam0, mu, phi1, sin1, cos1, tan1, phi, lam, t1, t2, C1, T1, N1, R1, D
a := p.a
esq := p.esq
epsq := p.epsq
lllist := list()
every i := 1 to *L by 3 do {
zone := L[i]
x := L[i + 1]
y := L[i + 2]
lam0 := dtor(-183 + 6 * zone) # central meridian of zone
mu := y / (k0 * p.c0)
phi1 := mu + p.c2 * sin(2. * mu) + p.c4 * sin(4. * mu) +
p.c6 * sin(6. * mu) + p.c8 * sin(8. * mu)
sin1 := sin(phi1)
cos1 := cos(phi1)
tan1 := tan(phi1)
t1 := 1 - esq * sin1^2
t2 := sqrt(t1)
C1 := epsq * cos1^2
T1 := tan1^2
N1 := a / t2
R1 := a * (1 - esq) / (t1 * t2)
D := x / (N1 * k0)
phi := phi1 - (N1 * tan1 / R1) *
(D^2 / 2. - (5. + 3.*T1 + 10.*C1 - 4.*C1*C1 - 9.*epsq) * D^4 / 24. +
(61. + 90.*T1 + 298.*C1 + 45.*T1*T1 - 252.*epsq - 3. * C1*C1) *
D^6 / 720.)
lam := lam0 + (D - (1 + 2 * T1 + C1) * D^3 / 6. +
(5. - 2. * C1 + 28. * T1 - 3. * C1 * C1 +
8. * epsq + 24. * T1 * T1) * D^5 / 120.) / cos1
put(lllist, rtod(lam), rtod(phi))
}
return lllist
end
# ctg_iutm_inv(p, L) -- invert inverse UTM projection
procedure ctg_iutm_inv(p)
return utm(p.a, p.f)
end
################## Composing projections #############################
record ctg_comp( # composition of two projections
proj, # projection procedure (always ctg_comp_proj)
inv, # inverse (always ctg_comp_inv)
projList # list of projections in composition,
# first is applied first, etc.
)
# compose -- produce a projection that applies the LAST projection
# in a[] first, etc.
procedure compose(a[]) #: define composite projection
local q, r
q := ctg_comp()
q.proj := ctg_comp_proj
q.inv := ctg_comp_inv
q.projList := []
every r := !a do
push(q.projList, r)
return q
end
procedure ctg_comp_proj(p, L)
local r
every r := !(p.projList) do
L := project(r, L)
return L
end
procedure ctg_comp_inv(p)
local q, r
q := ctg_comp()
q.proj := ctg_comp_proj
q.inv := ctg_comp_inv
q.projList := []
every r := !(p.projList) do
push(q.projList, invp(r))
return q
end
This page produced by UniDoc on 2021/04/15 @ 23:59:44.