[bknr-cvs] r2392 - branches/bos/projects/bos/m2
ksprotte at common-lisp.net
ksprotte at common-lisp.net
Wed Jan 23 07:13:32 UTC 2008
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)))))
+
+
More information about the Bknr-cvs
mailing list