summaryrefslogtreecommitdiff
path: root/ipl/procs/cartog.icn
diff options
context:
space:
mode:
Diffstat (limited to 'ipl/procs/cartog.icn')
-rw-r--r--ipl/procs/cartog.icn533
1 files changed, 533 insertions, 0 deletions
diff --git a/ipl/procs/cartog.icn b/ipl/procs/cartog.icn
new file mode 100644
index 0000000..010ebc9
--- /dev/null
+++ b/ipl/procs/cartog.icn
@@ -0,0 +1,533 @@
+############################################################################
+#
+# File: cartog.icn
+#
+# Subject: Procedures for cartographic projection
+#
+# Authors: Gregg M. Townsend and William S. Evans
+#
+# Date: February 19, 2003
+#
+############################################################################
+#
+# 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
+#
+############################################################################
+#
+# 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
+
+
+
+# sbsize(p, x, y, u, maxl) -- calculate scalebar size
+
+procedure sbsize(p, x, y, u, maxl) #: calculate scalebar size
+ local d, i, m, r
+
+ m := 1
+ repeat {
+ r := project(p, [x, y, x + m * u, y])
+ d := r[3] - r[1]
+ if d > maxl then
+ m := m / 10.0
+ else if d * 10 >= maxl
+ then break
+ else
+ m := m * 10
+ }
+
+ if maxl >= d * (i := 5 | 4 | 3 | 2) then
+ m *:= i
+
+ return m
+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
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-15 13:13:29 +0000