############################################################################ # # 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