Skip to content

Commit

Permalink
Merge branch 'master' of https://github.com/jl2/utm into type-fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@jl2
jl2 committed Feb 10, 2020
2 parents 0ee04f9 + de0f0f1 commit 82d2684
Show file tree
Hide file tree
Showing 2 changed files with 37 additions and 32 deletions.
20 changes: 10 additions & 10 deletions t/package.lisp
View file
Expand Up @@ -110,14 +110,14 @@

(test decimal-deg-min-sec
;; TODO: Add some more test cases...
(is-true (utm-near (deg-min-sec-to-decimal 20 30 40) 20.51111))
(multiple-value-bind (deg min sec) (utm:decimal-to-deg-min-sec 20.5111119)
(is-true (utm-near deg 20.0))
(is-true (utm-near min 30.0))
(is-true (utm-near sec 40.003967)))
(is-true (utm-near (deg-min-sec-to-decimal 78 39 10.8) 78.653))
(multiple-value-bind (deg min sec) (utm:decimal-to-deg-min-sec 78.653)
(is-true (utm-near deg 78.0))
(is-true (utm-near min 39.0))
(is-true (utm-near sec 10.8)))
(is-true (utm-near (deg-min-sec-to-decimal 20 30 40.0d0) 20.51111d0))
(multiple-value-bind (deg min sec) (utm:decimal-to-deg-min-sec 20.5111119d0)
(is-true (utm-near deg 20.0d0))
(is-true (utm-near min 30.0d0))
(is-true (utm-near sec 40.003967d0)))
(is-true (utm-near (deg-min-sec-to-decimal 78 39 10.8d0) 78.653d0))
(multiple-value-bind (deg min sec) (utm:decimal-to-deg-min-sec 78.653d0)
(is-true (utm-near deg 78.0d0))
(is-true (utm-near min 39.0d0))
(is-true (utm-near sec 10.8d0)))
)
49 changes: 27 additions & 22 deletions utm.lisp
View file
Expand Up @@ -64,7 +64,7 @@
;; * http://www.uwgb.edu/dutchs/FieldMethods/UTMSystem.htm
(defun lat-lon-to-utm (lat lon &key (ellipsoid "WGS84") (zone nil))
"Convert a point given as latitude and longitude into UTM using the specified ellipsoid. The default ellipsoid is WGS84."
(declare (optimize (speed 3) (space 0) (safety 0) (debug 3))
(declare (optimize (speed 3))
(type double-float lat lon)
(type (or null fixnum) zone))
(let* ((lat-rad (deg2rad lat))
Expand All @@ -75,15 +75,15 @@
(e-squared (* f (- 2.0d0 f)))
(e-prime-squared (/ e-squared (- 1.0d0 e-squared)))
(nzone (if zone zone (ceiling (/ (+ lon 180.0d0) 6.0d0))))
(long0 (deg2rad (- (* 6.0d0 (- nzone 1)) 177) ))
(long0 (deg2rad (- (* 6.0d0 (1- nzone)) 177) ))
;;(rho (/ (* a (- 1 e-squared)) (expt (- 1.0 (* e-squared (expt (sin lat-rad) 2.0))) (/ 3.0 2.0))))
(nu (realpart (/ a (expt (- 1.0d0 (* e-squared (expt (sin lat-rad) 2.0d0))) 0.5d0))))
(p (- lon-rad long0))

(M (* a (+ (* (- 1.0d0 (/ e-squared 4.0d0) (* (/ 3.0d0 64.0d0) (expt e-squared 2.0d0)) (* (/ 5.0d0 256.0d0) (expt e-squared 3.0d0))) lat-rad)
(M (realpart (* a (+ (* (- 1.0d0 (/ e-squared 4.0d0) (* (/ 3.0d0 64.0d0) (expt e-squared 2.0d0)) (* (/ 5.0d0 256.0d0) (expt e-squared 3.0d0))) lat-rad)
(- (* (sin (* 2.0d0 lat-rad)) (+ (* (/ 3.0d0 8.0d0) e-squared) (* (/ 3.0d0 32.0d0) (expt e-squared 2.0d0)) (* (/ 45.0d0 1024.0d0) (expt e-squared 3.0d0)))))
(* (sin (* 4.0d0 lat-rad)) (+ (* (/ 15.0d0 256.0d0) (expt e-squared 2.0d0)) (* (/ 45.0d0 1024.0d0) (expt e-squared 3.0d0))))
(- (* (sin (* 6.0d0 lat-rad)) (* (/ 35.0d0 3072.0d0) (expt e-squared 3.0d0)))))))
(- (* (sin (* 6.0d0 lat-rad)) (* (/ 35.0d0 3072.0d0) (expt e-squared 3.0d0))))))))

(K1 (* M k0))
(k2 (* k0 nu (sin (* 2.0d0 lat-rad)) 0.250d0))
Expand Down Expand Up @@ -116,49 +116,54 @@
(f (/ (- a b) a))
(e-squared (* f (- 2.0d0 f)))
(e-prime-squared (/ e-squared (- 1.0d0 e-squared)))
(long0 (deg2rad (- (* 6.0d0 (- zone 1)) 177) ))
(long0 (deg2rad (- (* 6.0d0 (1- zone)) 177) ))
(M (/ northing k0))
(mu (/ M (* a (- 1.0d0 (/ e-squared 4.0d0) (* (/ 3.0d0 64.0d0) (expt e-squared 2.0d0)) (* (/ 5.0d0 256.0d0) (expt e-squared 3.0d0))))))
(mu (realpart (/ M (* a (- 1.0d0 (/ e-squared 4.0d0) (* (/ 3.0d0 64.0d0) (expt e-squared 2.0d0)) (* (/ 5.0d0 256.0d0) (expt e-squared 3.0d0)))))))
(e1 (/ (- 1.0d0 (sqrt (- 1.0d0 e-squared))) (+ 1.0d0 (sqrt (- 1.0d0 e-squared)))))
(j1 (- (* (/ 3.0d0 2.0d0) e1) (* (/ 27.0d0 32.0d0) (expt e1 3.0d0))))
(j2 (- (* (/ 21.0d0 16.0d0) (expt e1 2.0d0)) (* (/ 55.0d0 32.0d0) (expt e1 4.0d0))))
(j3 (* (/ 151 96) (expt e1 3.0d0)))
(j4 (* (/ 1097.0d0 512.0d0) (expt e1 4.0d0)))
(j1 (realpart (- (* (/ 3.0d0 2.0d0) e1) (* (/ 27.0d0 32.0d0) (expt e1 3.0d0)))))
(j2 (realpart (- (* (/ 21.0d0 16.0d0) (expt e1 2.0d0)) (* (/ 55.0d0 32.0d0) (expt e1 4.0d0)))))
(j3 (realpart (* (/ 151 96) (expt e1 3.0d0))))
(j4 (realpart (* (/ 1097.0d0 512.0d0) (expt e1 4.0d0))))
(fp (+ mu
(* j1 (sin (* 2.0d0 mu )))
(* j2 (sin (* 4.0d0 mu )))
(* j3 (sin (* 6.0d0 mu )))
(* j4 (sin (* 8.0d0 mu )))))
(c1 (* e-prime-squared (expt (cos fp) 2.0d0)))
(t1 (expt (tan fp) 2.0d0))
(r1 (* a (/ (- 1.0d0 e-squared) (expt (- 1.0d0 (* e-squared (expt (sin fp) 2.0d0))) (/ 3.0d0 2.0d0)))))
(n1 (/ a (expt (- 1.0d0 (* e-squared (expt (sin fp) 2.0d0))) 0.5d0)))
(c1 (realpart (* e-prime-squared (expt (cos fp) 2.0d0))))
(t1 (realpart (expt (tan fp) 2.0d0)))
(r1 (realpart (* a (/ (- 1.0d0 e-squared) (expt (- 1.0d0 (* e-squared (expt (sin fp) 2.0d0))) (/ 3.0d0 2.0d0))))))
(n1 (realpart (/ a (expt (- 1.0d0 (* e-squared (expt (sin fp) 2.0d0))) 0.5))))
(D (/ reasting (* n1 k0)))
(q1 (* n1 (/ (tan fp) r1)))
(q2 (/ (expt d 2.0d0) 2.0d0))
(q3 (/ (* (+ 5 (* 3.0d0 t1) (* 10.0d0 c1) (* -4 (expt c1 2.0d0)) (* -9.0d0 e-prime-squared)) (expt D 4.0d0)) 24.0d0))
(q4 (/ (* (+ 61 (* 90.0d0 t1) (* 298.0d0 c1) (* 45.0d0 (expt t1 2.0d0)) (* -3.0d0 (expt c1 2.0d0)) (* -252 e-prime-squared)) (expt D 6.0d0)) 720.0d0))
(q2 (realpart (/ (expt d 2.0d0) 2.0d0)))
(q3 (realpart (/ (* (+ 5 (* 3.0d0 t1) (* 10.0d0 c1) (* -4 (expt c1 2.0d0)) (* -9.0d0 e-prime-squared)) (expt D 4.0d0)) 24.0d0)))
(q4 (realpart (/ (* (+ 61 (* 90.0d0 t1) (* 298.0d0 c1) (* 45.0d0 (expt t1 2.0d0)) (* -3.0d0 (expt c1 2.0d0)) (* -252 e-prime-squared)) (expt D 6.0d0)) 720.0d0)))
(q5 d)
(q6 (/ (* (+ 1.0d0 (* 2.0d0 t1) c1) (expt d 3.0d0)) 6.0d0))
(q7 (/ (* (+ 5.0d0 (* -2.0d0 c1) (* 28.0d0 t1) (* -3 (expt c1 2.0d0)) (* 8.0d0 e-prime-squared) (* 24.0d0 (expt t1 2.0d0))) (expt d 5.0d0)) 120.0d0))
(q6 (realpart (/ (* (+ 1.0d0 (* 2.0d0 t1) c1) (expt d 3.0d0)) 6.0d0)))
(q7 (realpart (/ (* (+ 5.0d0 (* -2.0d0 c1) (* 28.0d0 t1) (* -3 (expt c1 2.0d0)) (* 8.0d0 e-prime-squared) (* 24.0d0 (expt t1 2.0d0))) (expt d 5.0d0)) 120.0d0)))
(lat (- fp (* q1 (+ q2 (- q3) q4))))
(lon (+ long0 (/ (+ q5 (- q6) q7) (cos fp)))))
(declare (type double-float reasting a b f e-squared e-prime-squared long0 M mu
e1 j1 j2 j3 j4 fp c1 t1 r1 n1 D q1 q2 q3 q4 q5 q6 q7 lat lon))
(list (rad2deg (realpart lat)) (rad2deg (realpart lon)))))

(defun deg-min-sec-to-decimal (degree minute second)
"Convert degree, minute, second format to decimal."
(declare (optimize (speed 3)))
"Convert degree, minute, second format to decimal.nn"
(declare (optimize (speed 3))
(type number degree minute second))
(+ degree (/ minute 60.0d0) (/ second (* 60.0d0 60.0d0))))

(defun decimal-to-deg-min-sec (decimal)
"Convert degree, minute, second format to decimal."
(declare (optimize (speed 3))
(type double-float decimal))
(type number decimal))
(multiple-value-bind (degrees min-secs) (truncate decimal)
(declare (type fixnum degrees)
(type double-float min-secs))
(when (< degrees 0)
(setf min-secs (- min-secs))
(decf degrees))
(multiple-value-bind (minutes secs) (truncate (* 60 min-secs))
(declare (type fixnum minutes)
(type double-float secs))
(values degrees minutes (* secs 60.0d0)))))

0 comments on commit 82d2684

Please sign in to comment.