Author: ksprotte Date: Wed Jan 23 02:13:31 2008 New Revision: 2392 Modified: branches/bos/projects/bos/m2/geo-utm.lisp Log: The latest snapshot of the geoutm.js conversion effort. It seems to work so far, but still misses error checking of the arguments to be converted. Modified: branches/bos/projects/bos/m2/geo-utm.lisp ============================================================================== --- branches/bos/projects/bos/m2/geo-utm.lisp (original) +++ branches/bos/projects/bos/m2/geo-utm.lisp Wed Jan 23 02:13:31 2008 @@ -10,143 +10,224 @@ ;; <P>Programmers: The JavaScript source code in this document may be copied ;; and reused without restriction.</P> -(defconstant sm_a 6378137.0) -(defconstant sm_b 6356752.314) -(defconstant sm_EccSquared 6.69437999013e-03) - -(defconstant UTMScaleFactor 0.9996) - -(defun DegToRad (deg) - (* (/ deg 180.0) pi)) - -(defun RadToDeg (rad) - (/ rad (* pi 180.0))) - -(defun ArcLengthOfMeridian (phi) - (let* ((n (/ (- sm_a sm_b) (+ sm_a sm_b))) - (alpha (* (/ (+ sm_a sm_b) 2) - (+ 1 - (/ (expt n 2) 4) - (/ (expt n 4) 64)))) - (beta (+ (/ (* -3 n) 2) - (/ (* 9 (expt n 3)) 16) - (/ (* -3 (expt n 5)) 32))) - (gamma (+ (/ (* 15 (expt n 2)) 16) - (/ (* -15 (expt n 4)) 32))) - (delta (+ (/ (* -35 (expt n 3)) 48) - (/ (* 105 (expt n 5)) 256))) - (epsilon (/ (* 315 (expt n 4)) 512))) - (* alpha - (+ phi - (* beta (sin (* 2 phi))) - (* gamma (sin (* 4 phi))) - (* delta (sin (* 6 phi))) - (* epsilon (sin (* 8 phi))))))) - -(defun UTMCentralMeridian (zone) - (DegToRad (+ -183 (* zone 6)))) - -(defun FootpointLatitude (y) - (let* ((n (/ (- sm_a sm_b) (+ sm_a sm_b))) - (alpha (* (/ (+ sm_a sm_b) 2) - (+ 1 - (/ (expt n 2) 4) - (/ (expt n 4) 64)))) - (y (/ y alpha)) - (beta (+ (/ (* 3 n) 2) - (/ (* -27 (expt n 3)) 32) - (/ (* 269 (expt n 5) 512)))) - (gamma (+ (/ (* 21 (expt n 2)) 16) - (/ (* -55 (expt n 4)) 32))) - (delta (+ (/ (* 151 (expt n 3)) 96) - (/ (* -417 (expt n 5)) 128))) - (epsilon (/ (* 1097 (expt n 4)) 512))) - (+ y - (* beta (sin (* 2 y))) - (* gamma (sin (* 4 y))) - (* delta (sin (* 6 y))) - (* epsilon (sin (* 8 y)))))) - -(defun MapLatLonToXY (phi lambda lambda0) - (let* ((ep2 (/ (- (expt sm_a 2) (expt sm_b 2)) (expt sm_b 2))) - (nu2 (* ep2 (expt (cos phi) 2))) - (N (/ (expt sm_a 2) (* sm_b (sqrt (+ 1 nu2))))) - (t- (tan phi)) - (t2 (* t- t-)) - (l (- lambda lambda0)) - (l3coef (+ 1 (- t2) nu2)) - (l4coef (+ 5 (- t2) (* 9 nu2) (* 4 nu2 nu2))) - (l5coef (+ 5 (- (* 18 t2)) (* t2 t2) (* 14 nu2) (- (* 58 t2 nu2)))) - (l6coef (+ 61 (- (* 58 t2)) (* t2 t2) (* 270 nu2) (- (* 330 t2 nu2)))) - (l7coef (+ 61 (- (* 479 t2)) (* 179 t2 t2) (- (* t2 t2 t2)))) - (l8coef (+ 1385 (- (* 3111 t2)) (* 543 t2 t2) (- (* t2 t2 t2)))) - (easting (+ (* N (cos phi) l) - (/ N (* 6 (expt (cos phi) 3) l3coef (expt l 3))) - (/ N (* 120 (expt (cos phi) 5) l5coef (expt l 5))) - (/ N (* 5040 (expt (cos phi) 7) l7coef (expt l 7))))) - (northing (+ (ArcLengthOfMeridian phi) - (/ t- (* 2 N (expt (cos phi) 2) (expt l 2))) - (/ t- (* 24 N (expt (cos phi) 4) l4coef (expt l 4))) - (/ t- (* 720 N (expt (cos phi) 6) l6coef (expt l 6))) - (/ t- (* 40320 N (expt (cos phi) 8) l8coef (expt l 8)))))) - (list easting northing))) - -(defun MapXYToLatLon (x y lambda0) - (let* ((phif (FootpointLatitude y)) - (ep2 (/ (- (expt sm_a 2) (expt sm_b 2)) (expt sm_b 2))) - (cf (cos phif)) - (nuf2 (* ep2 (expt cf 2))) - (Nf (/ (expt sm_a 2) (* sm_b (sqrt (+ 1 nuf2))))) - (tf (tan phif)) - (tf2 (* tf tf)) - (tf4 (* tf2 tf2)) - (x1frac (/ 1.0 (* Nf cf))) - (x2frac (/ tf (* 2 (expt Nf 2)))) - (x3frac (/ 1.0 (* 6 (expt Nf 3) cf))) - (x4frac (/ tf (* 24 (expt Nf 4)))) - (x5frac (/ 1.0 (* 120 (expt Nf 5) cf))) - (x6frac (/ tf (* 720 (expt Nf 6)))) - (x7frac (/ 1.0 (* 5040 (expt Nf 7) cf))) - (x8frac (/ tf (* 40320 (expt Nf 8)))) - (x2poly (- -1 nuf2)) - (x3poly (- -1 (* 2 tf2) nuf2)) - (x4poly (+ 5 - (* 3 tf2) - (* 6 nuf2) - (- (* 6 tf2 nuf2)) - (- (* 3 nuf2 nuf2)) - (- (* 9 tf2 nuf2 nuf2)))) - (x5poly (+ 5 (* 28 tf2) (* 24 tf4) (* 6 nuf2) (* 8 tf2 nuf2))) - (x6poly (+ -61 (- (* 90 tf2)) (- (* 45 tf4)) (- (* 107 nuf2)) (* 162 tf2 nuf2))) - (x7poly (- 6 (* 662 tf2) (* 1320 tf4) (* 720 tf4 tf2))) - (x8poly (+ 1385 (* 3633 tf2) (* 4095 tf4) (* 1575 tf4 tf2))) - (latitude (+ phif - (* x2frac x2poly x x) - (* x4frac x4poly (expt x 4)) - (* x6frac x6poly (expt x 6)) - (* x8frac x8poly (expt x 8)))) - (longitude (+ lambda0 - (* x1frac x) - (* x3frac x3poly (expt x 3)) - (* x5frac x5poly (expt x 5)) - (* x7frac x7poly (expt x 7))))) - (list latitude longitude))) - -(defun lat-lon-to-utm-x-y (lat lon &optional (zone 50)) - ;; The Javascript version claims that the zone is calculated if not - ;; provided, but I could not find the code that does it. This - ;; should be added. This version of the code requires that the - ;; correct zone is provided. - (destructuring-bind (easting northing) - (MapLatLonToXY (DegToRad lat) (DegToRad lon) (UTMCentralMeridian zone)) - (list (+ (* easting UTMScaleFactor) 500000) - (+ (* northing UTMScaleFactor) - (if (minusp northing) 10000000 0)) - zone))) - -(defun utm-x-y-to-lat-lon (x y zone southhemi) - (destructuring-bind (lat lon) - (MapXYToLatLon (/ (- x 500000) UTMScaleFactor) - (/ (- y (if southhemi 10000000 0)) UTMScaleFactor) - (UTMCentralMeridian zone)) - (list (RadToDeg lat) (RadToDeg lon)))) \ No newline at end of file +(defconstant sm-a 6378137.0) +(defconstant sm-b 6356752.314) +(defconstant sm-eccsquared 6.69437999013e-03) + +(defconstant utmscale-factor 0.9996) + +(define-modify-macro multiplyf (x) *) + +(define-modify-macro dividef (x) /) + +(defun deg-to-rad (deg) (* (/ deg 180.0) pi)) + +(defun rad-to-deg (rad) (* (/ rad pi) 180.0)) + +(defun arc-length-of-meridian (phi) + (let (alpha beta gamma delta epsilon n) + (let (result) + (setq n (/ (- sm-a sm-b) (+ sm-a sm-b))) + (setq alpha + (* (/ (+ sm-a sm-b) 2.0) + (+ (+ 1.0 (/ (expt n 2.0) 4.0)) + (/ (expt n 4.0) 64.0)))) + (setq beta + (+ + (+ (/ (* (- 3.0) n) 2.0) (/ (* 9.0 (expt n 3.0)) 16.0)) + (/ (* (- 3.0) (expt n 5.0)) 32.0))) + (setq gamma + (+ (/ (* 15.0 (expt n 2.0)) 16.0) + (/ (* (- 15.0) (expt n 4.0)) 32.0))) + (setq delta + (+ (/ (* (- 35.0) (expt n 3.0)) 48.0) + (/ (* 105.0 (expt n 5.0)) 256.0))) + (setq epsilon (/ (* 315.0 (expt n 4.0)) 512.0)) + (setq result + (* alpha + (+ + (+ + (+ (+ phi (* beta (sin (* 2.0 phi)))) + (* gamma (sin (* 4.0 phi)))) + (* delta (sin (* 6.0 phi)))) + (* epsilon (sin (* 8.0 phi)))))) + result))) + +(defun utmcentral-meridian (zone) + (let (cmeridian) + (setq cmeridian (deg-to-rad (+ (- 183.0) (* zone 6.0)))) + cmeridian)) + +(defun footpoint-latitude (y) + (let (y_ alpha_ beta_ gamma_ delta_ epsilon_ n) + (let (result) + (setq n (/ (- sm-a sm-b) (+ sm-a sm-b))) + (setq alpha_ + (* (/ (+ sm-a sm-b) 2.0) + (+ (+ 1 (/ (expt n 2.0) 4)) (/ (expt n 4.0) 64)))) + (setq y_ (/ y alpha_)) + (setq beta_ + (+ + (+ (/ (* 3.0 n) 2.0) (/ (* (- 27.0) (expt n 3.0)) 32.0)) + (/ (* 269.0 (expt n 5.0)) 512.0))) + (setq gamma_ + (+ (/ (* 21.0 (expt n 2.0)) 16.0) + (/ (* (- 55.0) (expt n 4.0)) 32.0))) + (setq delta_ + (+ (/ (* 151.0 (expt n 3.0)) 96.0) + (/ (* (- 417.0) (expt n 5.0)) 128.0))) + (setq epsilon_ (/ (* 1097.0 (expt n 4.0)) 512.0)) + (setq result + (+ + (+ + (+ (+ y_ (* beta_ (sin (* 2.0 y_)))) + (* gamma_ (sin (* 4.0 y_)))) + (* delta_ (sin (* 6.0 y_)))) + (* epsilon_ (sin (* 8.0 y_))))) + result))) + +(defun map-lat-lon-to-xy (phi lambda lambda0) + (let (n nu2 ep2 %t t2 l) + (let (l3coef l4coef l5coef l6coef l7coef l8coef) + (let (tmp) + (setq ep2 + (/ (- (expt sm-a 2.0) (expt sm-b 2.0)) + (expt sm-b 2.0))) + (setq nu2 (* ep2 (expt (cos phi) 2.0))) + (setq n (/ (expt sm-a 2.0) (* sm-b (sqrt (+ 1 nu2))))) + (setq %t (tan phi)) + (setq t2 (* %t %t)) + (setq tmp (- (* (* t2 t2) t2) (expt %t 6.0))) + (setq l (- lambda lambda0)) + (setq l3coef (+ (- 1.0 t2) nu2)) + (setq l4coef (+ (+ (- 5.0 t2) (* 9 nu2)) (* 4.0 (* nu2 nu2)))) + (setq l5coef + (- (+ (+ (- 5.0 (* 18.0 t2)) (* t2 t2)) (* 14.0 nu2)) + (* (* 58.0 t2) nu2))) + (setq l6coef + (- (+ (+ (- 61.0 (* 58.0 t2)) (* t2 t2)) (* 270.0 nu2)) + (* (* 330.0 t2) nu2))) + (setq l7coef + (- (+ (- 61.0 (* 479.0 t2)) (* 179.0 (* t2 t2))) + (* (* t2 t2) t2))) + (setq l8coef + (- (+ (- 1385.0 (* 3111.0 t2)) (* 543.0 (* t2 t2))) + (* (* t2 t2) t2))) + (values + (+ + (+ + (+ (* (* n (cos phi)) l) + (* (* (* (/ n 6.0) (expt (cos phi) 3.0)) l3coef) + (expt l 3.0))) + (* (* (* (/ n 120.0) (expt (cos phi) 5.0)) l5coef) + (expt l 5.0))) + (* (* (* (/ n 5040.0) (expt (cos phi) 7.0)) l7coef) + (expt l 7.0))) + (+ + (+ + (+ + (+ (arc-length-of-meridian phi) + (* (* (* (/ %t 2.0) n) (expt (cos phi) 2.0)) + (expt l 2.0))) + (* (* (* (* (/ %t 24.0) n) (expt (cos phi) 4.0)) l4coef) + (expt l 4.0))) + (* (* (* (* (/ %t 720.0) n) (expt (cos phi) 6.0)) l6coef) + (expt l 6.0))) + (* (* (* (* (/ %t 40320.0) n) (expt (cos phi) 8.0)) l8coef) + (expt l 8.0)))))))) + +(defun map-xyto-lat-lon (x y lambda0) + (let (phif nf nfpow nuf2 ep2 tf tf2 tf4 cf) + (let (x1frac x2frac x3frac x4frac x5frac x6frac x7frac x8frac) + (let (x2poly x3poly x4poly x5poly x6poly x7poly x8poly) + (setq phif (footpoint-latitude y)) + (setq ep2 + (/ (- (expt sm-a 2.0) (expt sm-b 2.0)) + (expt sm-b 2.0))) + (setq cf (cos phif)) + (setq nuf2 (* ep2 (expt cf 2.0))) + (setq nf (/ (expt sm-a 2.0) (* sm-b (sqrt (+ 1 nuf2))))) + (setq nfpow nf) + (setq tf (tan phif)) + (setq tf2 (* tf tf)) + (setq tf4 (* tf2 tf2)) + (setq x1frac (/ 1.0 (* nfpow cf))) + (multiplyf nfpow nf) + (setq x2frac (/ tf (* 2.0 nfpow))) + (multiplyf nfpow nf) + (setq x3frac (/ 1.0 (* (* 6.0 nfpow) cf))) + (multiplyf nfpow nf) + (setq x4frac (/ tf (* 24.0 nfpow))) + (multiplyf nfpow nf) + (setq x5frac (/ 1.0 (* (* 120.0 nfpow) cf))) + (multiplyf nfpow nf) + (setq x6frac (/ tf (* 720.0 nfpow))) + (multiplyf nfpow nf) + (setq x7frac (/ 1.0 (* (* 5040.0 nfpow) cf))) + (multiplyf nfpow nf) + (setq x8frac (/ tf (* 40320.0 nfpow))) + (setq x2poly (- (- 1.0) nuf2)) + (setq x3poly (- (- (- 1.0) (* 2 tf2)) nuf2)) + (setq x4poly + (- + (- + (- (+ (+ 5.0 (* 3.0 tf2)) (* 6.0 nuf2)) + (* (* 6.0 tf2) nuf2)) + (* 3.0 (* nuf2 nuf2))) + (* (* 9.0 tf2) (* nuf2 nuf2)))) + (setq x5poly + (+ + (+ (+ (+ 5.0 (* 28.0 tf2)) (* 24.0 tf4)) (* 6.0 nuf2)) + (* (* 8.0 tf2) nuf2))) + (setq x6poly + (+ + (- (- (- (- 61.0) (* 90.0 tf2)) (* 45.0 tf4)) + (* 107.0 nuf2)) + (* (* 162.0 tf2) nuf2))) + (setq x7poly + (- (- (- (- 61.0) (* 662.0 tf2)) (* 1320.0 tf4)) + (* 720.0 (* tf4 tf2)))) + (setq x8poly + (+ (+ (+ 1385.0 (* 3633.0 tf2)) (* 4095.0 tf4)) + (* 1575 (* tf4 tf2)))) + (values + (+ + (+ + (+ (+ phif (* (* x2frac x2poly) (* x x))) + (* (* x4frac x4poly) (expt x 4.0))) + (* (* x6frac x6poly) (expt x 6.0))) + (* (* x8frac x8poly) (expt x 8.0))) + (+ + (+ + (+ (+ lambda0 (* x1frac x)) + (* (* x3frac x3poly) (expt x 3.0))) + (* (* x5frac x5poly) (expt x 5.0))) + (* (* x7frac x7poly) (expt x 7.0)))))))) + +(defun lat-lon-to-utm-x-y (lat lon) + "Returns four values X, Y, ZONE and SOUTHHEMI-P." + (let* ((lat (float lat 0d0)) + (lon (float lon 0d0)) + (zone (+ (floor (/ (+ lon 180.0) 6)) 1))) + (multiple-value-bind (x y) + (map-lat-lon-to-xy (deg-to-rad lat) (deg-to-rad lon) (utmcentral-meridian zone)) + (setq x (+ (* x utmscale-factor) 500000.0)) + (setq y (* y utmscale-factor)) + (if (< y 0.0) (block nil (setq y (+ y 1.e7))) nil) + (values x y zone (minusp lat))))) + +(defun utm-x-y-to-lat-lon (x y zone southhemi-p) + "Returns two values LAT and LON." + (let ((x (float x 0d0)) + (y (float y 0d0)) + cmeridian) + (decf x 500000.0) + (dividef x utmscale-factor) + (if southhemi-p (decf y 1.e7) nil) + (dividef y utmscale-factor) + (setq cmeridian (utmcentral-meridian zone)) + (multiple-value-bind (lat lon) + (map-xyto-lat-lon x y cmeridian) + (values (rad-to-deg lat) (rad-to-deg lon))))) + +