geo-utm.lisp

Revision 3656, 12.4 kB (checked in by ksprotte, 4 months ago)
whitespace cleanup

Line
1	(in-package :geo-utm)
2
3	;; Converted from Javascript originally almost automatically
4	;; has now been tuned manually
5
6	;; Origin: http:;;home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html
7	;; Copyright 1997-1998 by Charles L. Taylor
8
9	;; Page says:
10
11	;; <P>Programmers: The JavaScript source code in this document may be copied
12	;; and reused without restriction.</P>
13
14	(eval-when (:compile-toplevel :load-toplevel :execute)
15	(defparameter initial-read-default-float-format read-default-float-format)
16	(setq read-default-float-format 'double-float))
17
18	(deftype limited-double-float ()
19	'(double-float
20	(#.(- (expt 2d0 64)))
21	(#.(expt 2d0 64))))
22
23	(deftype float-pair ()
24	'(simple-array double-float (2)))
25
26	(defun make-float-pair ()
27	(make-array 2 :element-type 'double-float :initial-element 0d0))
28
29	(defconstant sm-a 6378137.0)
30	(defconstant sm-b 6356752.314)
31	(defconstant sm-eccsquared 6.69437999013e-03)
32
33	(defconstant utmscale-factor 0.9996)
34
35	(define-modify-macro multiplyf (x) *)
36	(define-modify-macro dividef (x) /)
37
38	(declaim (ftype (function (double-float double-float float-pair)
39	(values float-pair t t))
40	lon-lat-to-utm-x-y*))
41
42	(locally
43	(declare (optimize (speed 3) (safety 1) (debug 1)))
44	(macrolet ((let-double-floats (symbols &body body)
45	`(let (,@(mapcar #'(lambda (symbol) `(,symbol 0d0)) symbols))
46	(declare (double-float ,@symbols))
47	,@body)))
48	(labels ((deg-to-rad (deg)
49	(declare (double-float deg))
50	(* (/ deg 180.0) pi))
51
52	(rad-to-deg (rad)
53	(declare (double-float rad))
54	(* (/ rad pi) 180.0))
55
56	(arc-length-of-meridian (phi)
57	(declare (double-float phi))
58	(let* ((n (/ (- sm-a sm-b) (+ sm-a sm-b)))
59	alpha beta gamma delta epsilon)
60	(setq alpha
61	(* (/ (+ sm-a sm-b) 2.0)
62	(+ (+ 1.0 (/ (expt n 2) 4.0))
63	(/ (expt n 4) 64.0))))
64	(setq beta
65	(+
66	(+ (/ (* (- 3.0) n) 2.0) (/ (* 9.0 (expt n 3)) 16.0))
67	(/ (* (- 3.0) (expt n 5)) 32.0)))
68	(setq gamma
69	(+ (/ (* 15.0 (expt n 2)) 16.0)
70	(/ (* (- 15.0) (expt n 4)) 32.0)))
71	(setq delta
72	(+ (/ (* (- 35.0) (expt n 3)) 48.0)
73	(/ (* 105.0 (expt n 5)) 256.0)))
74	(setq epsilon (/ (* 315.0 (expt n 4)) 512.0))
75	(* alpha
76	(+
77	(+
78	(+ (+ phi (* beta (sin (the limited-double-float (* 2.0 phi)))))
79	(* gamma (sin (the limited-double-float (* 4.0 phi)))))
80	(* delta (sin (the limited-double-float (* 6.0 phi)))))
81	(* epsilon (sin (the limited-double-float (* 8.0 phi))))))))
82
83	(utmcentral-meridian (zone)
84	(declare (fixnum zone))
85	(deg-to-rad (+ (- 183.0) (* zone 6.0))))
86
87	(footpoint-latitude (y)
88	(declare (double-float y))
89	(let* ((n (/ (- sm-a sm-b) (+ sm-a sm-b)))
90	(alpha_ (* (/ (+ sm-a sm-b) 2.0)
91	(+ (+ 1 (/ (expt n 2) 4)) (/ (expt n 4) 64))))
92	(y_ (/ y alpha_))
93	beta_ gamma_ delta_ epsilon_)
94	(setq beta_
95	(+
96	(+ (/ (* 3.0 n) 2.0) (/ (* (- 27.0) (expt n 3)) 32.0))
97	(/ (* 269.0 (expt n 5)) 512.0)))
98	(setq gamma_
99	(+ (/ (* 21.0 (expt n 2)) 16.0)
100	(/ (* (- 55.0) (expt n 4)) 32.0)))
101	(setq delta_
102	(+ (/ (* 151.0 (expt n 3)) 96.0)
103	(/ (* (- 417.0) (expt n 5)) 128.0)))
104	(setq epsilon_ (/ (* 1097.0 (expt n 4)) 512.0))
105	(+
106	(+
107	(+ (+ y_ (* beta_ (sin (the limited-double-float (* 2.0 y_)))))
108	(* gamma_ (sin (the limited-double-float (* 4.0 y_)))))
109	(* delta_ (sin (the limited-double-float (* 6.0 y_)))))
110	(* epsilon_ (sin (the limited-double-float (* 8.0 y_)))))))
111
112	(map-lat-lon-to-xy (phi lambda lambda0)
113	(declare (double-float phi lambda lambda0))
114	(let* ((ep2 (/ (- (expt sm-a 2) (expt sm-b 2))
115	(expt sm-b 2)))
116	(nu2 (* ep2 (expt (cos (the limited-double-float phi)) 2)))
117	(n (/ (expt sm-a 2) (* sm-b (sqrt (+ 1 nu2)))))
118	(%t (tan (the limited-double-float phi)))
119	(t2 (* %t %t))
120	;; (tmp (- (* (* t2 t2) t2) (expt %t 6)))
121	(l (- lambda lambda0))
122	(l3coef (+ (- 1.0 t2) nu2))
123	(l4coef (+ (+ (- 5.0 t2) (* 9 nu2)) (* 4.0 (* nu2 nu2))))
124	(l5coef (- (+ (+ (- 5.0 (* 18.0 t2)) (* t2 t2)) (* 14.0 nu2))
125	(* (* 58.0 t2) nu2)))
126	(l6coef (- (+ (+ (- 61.0 (* 58.0 t2)) (* t2 t2)) (* 270.0 nu2))
127	(* (* 330.0 t2) nu2)))
128	(l7coef (- (+ (- 61.0 (* 479.0 t2)) (* 179.0 (* t2 t2)))
129	(* (* t2 t2) t2)))
130	(l8coef (- (+ (- 1385.0 (* 3111.0 t2)) (* 543.0 (* t2 t2)))
131	(* (* t2 t2) t2)))
132	(l4 (expt (abs l) 4))
133	(l5 (* l l4))
134	(l6 (* l l5))
135	(l7 (* l l6))
136	(l8 (* l l7))
137	(cos-phi (cos (the limited-double-float phi)))
138	(expt-cos-phi-2 (expt cos-phi 2))
139	(expt-cos-phi-3 (* expt-cos-phi-2 cos-phi))
140	(expt-cos-phi-4 (* expt-cos-phi-3 cos-phi))
141	(expt-cos-phi-5 (* expt-cos-phi-4 cos-phi))
142	(expt-cos-phi-6 (* expt-cos-phi-5 cos-phi))
143	(expt-cos-phi-7 (* expt-cos-phi-6 cos-phi))
144	(expt-cos-phi-8 (* expt-cos-phi-7 cos-phi)))
145	(values
146	(+
147	(+
148	(+ (* (* n (cos (the limited-double-float phi))) l)
149	(* (* (* (/ n 6.0) expt-cos-phi-3) l3coef)
150	(expt l 3)))
151	(* (* (* (/ n 120.0) expt-cos-phi-5) l5coef)
152	l5))
153	(* (* (* (/ n 5040.0) expt-cos-phi-7) l7coef)
154	l7))
155	(+
156	(+
157	(+
158	(+ (arc-length-of-meridian phi)
159	(* (* (* (/ %t 2.0) n) expt-cos-phi-2)
160	(expt l 2)))
161	(* (* (* (* (/ %t 24.0) n) expt-cos-phi-4) l4coef)
162	l4))
163	(* (* (* (* (/ %t 720.0) n) expt-cos-phi-6) l6coef)
164	l6))
165	(* (* (* (* (/ %t 40320.0) n) expt-cos-phi-8) l8coef)
166	l8)))))
167
168	(map-xyto-lat-lon (x y lambda0)
169	(declare (double-float x y lambda0))
170	(let-double-floats
171	(phif nf nfpow nuf2 ep2 tf tf2 tf4 cf
172	x1frac x2frac x3frac x4frac x5frac x6frac x7frac x8frac
173	x2poly x3poly x4poly x5poly x6poly x7poly x8poly)
174	(setq phif (footpoint-latitude y))
175	(setq ep2
176	(/ (- (expt sm-a 2) (expt sm-b 2))
177	(expt sm-b 2)))
178	(setq cf (cos (the limited-double-float phif)))
179	(setq nuf2 (* ep2 (expt cf 2)))
180	(setq nf (/ (expt sm-a 2) (* sm-b (sqrt (+ 1 nuf2)))))
181	(setq nfpow nf)
182	(setq tf (tan (the limited-double-float phif)))
183	(setq tf2 (* tf tf))
184	(setq tf4 (* tf2 tf2))
185	(setq x1frac (/ 1.0 (* nfpow cf)))
186	(multiplyf nfpow nf)
187	(setq x2frac (/ tf (* 2.0 nfpow)))
188	(multiplyf nfpow nf)
189	(setq x3frac (/ 1.0 (* (* 6.0 nfpow) cf)))
190	(multiplyf nfpow nf)
191	(setq x4frac (/ tf (* 24.0 nfpow)))
192	(multiplyf nfpow nf)
193	(setq x5frac (/ 1.0 (* (* 120.0 nfpow) cf)))
194	(multiplyf nfpow nf)
195	(setq x6frac (/ tf (* 720.0 nfpow)))
196	(multiplyf nfpow nf)
197	(setq x7frac (/ 1.0 (* (* 5040.0 nfpow) cf)))
198	(multiplyf nfpow nf)
199	(setq x8frac (/ tf (* 40320.0 nfpow)))
200	(setq x2poly (- (- 1.0) nuf2))
201	(setq x3poly (- (- (- 1.0) (* 2 tf2)) nuf2))
202	(setq x4poly
203	(-
204	(-
205	(- (+ (+ 5.0 (* 3.0 tf2)) (* 6.0 nuf2))
206	(* (* 6.0 tf2) nuf2))
207	(* 3.0 (* nuf2 nuf2)))
208	(* (* 9.0 tf2) (* nuf2 nuf2))))
209	(setq x5poly
210	(+
211	(+ (+ (+ 5.0 (* 28.0 tf2)) (* 24.0 tf4)) (* 6.0 nuf2))
212	(* (* 8.0 tf2) nuf2)))
213	(setq x6poly
214	(+
215	(- (- (- (- 61.0) (* 90.0 tf2)) (* 45.0 tf4))
216	(* 107.0 nuf2))
217	(* (* 162.0 tf2) nuf2)))
218	(setq x7poly
219	(- (- (- (- 61.0) (* 662.0 tf2)) (* 1320.0 tf4))
220	(* 720.0 (* tf4 tf2))))
221	(setq x8poly
222	(+ (+ (+ 1385.0 (* 3633.0 tf2)) (* 4095.0 tf4))
223	(* 1575 (* tf4 tf2))))
224	(let* ((x4 (expt (abs x) 4))
225	(x5 (* x4 x))
226	(x6 (* x5 x))
227	(x7 (* x6 x))
228	(x8 (* x7 x)))
229	(values
230	(+
231	(+
232	(+ (+ phif (* (* x2frac x2poly) (* x x)))
233	(* (* x4frac x4poly) x4))
234	(* (* x6frac x6poly) x6))
235	(* (* x8frac x8poly) x8))
236	(+
237	(+
238	(+ (+ lambda0 (* x1frac x))
239	(* (* x3frac x3poly) (expt x 3)))
240	(* (* x5frac x5poly) x5))
241	(* (* x7frac x7poly) x7)))))))
242
243	(defun lon-lat-to-utm-x-y* (lon lat float-pair)
244	"Returns list of four values X, Y, ZONE and SOUTHHEMI-P."
245	(declare ((double-float -180d0 180d0) lon)
246	((double-float -90d0 90d0) lat)
247	(float-pair float-pair))
248	(let ((zone (+ (floor (+ lon 180.0) 6.0) 1)))
249	(multiple-value-bind (x y)
250	(map-lat-lon-to-xy (deg-to-rad lat) (deg-to-rad lon) (utmcentral-meridian zone))
251	(setq x (+ (* x utmscale-factor) 500000.0))
252	(setq y (* y utmscale-factor))
253	(when (< y 0.0) (setq y (+ y 1.e7)))
254	(setf (aref float-pair 0) x
255	(aref float-pair 1) y)
256	(values float-pair zone (minusp lat)))))
257
258	(defun utm-x-y-to-lon-lat* (x y zone southhemi-p float-pair)
259	"Returns list (LON LAT)."
260	(declare (double-float x y) (float-pair float-pair))
261	(let ((cmeridian (utmcentral-meridian zone)))
262	(decf x 500000.0)
263	(dividef x utmscale-factor)
264	(if southhemi-p (decf y 1.e7) nil)
265	(dividef y utmscale-factor)
266	(multiple-value-bind (lat lon)
267	(map-xyto-lat-lon x y cmeridian)
268	(setf (aref float-pair 0) (rad-to-deg lon)
269	(aref float-pair 1) (rad-to-deg lat))
270	float-pair))))))
271
272	(defun lon-lat-to-utm-x-y (lon lat)
273	(multiple-value-bind (float-pair zone southhemi-p)
274	(lon-lat-to-utm-x-y* (float lon 0d0) (float lat 0d0) (make-float-pair))
275	(list (aref float-pair 0) (aref float-pair 1)
276	zone southhemi-p)))
277
278	(defun utm-x-y-to-lon-lat (x y zone southhemi-p)
279	"Returns list (LON LAT)."
280	(let ((lon-lat (utm-x-y-to-lon-lat* (float x 0d0) (float y 0d0)
281	zone southhemi-p (make-float-pair))))
282	(list (aref lon-lat 0) (aref lon-lat 1))))
283
284	(eval-when (:compile-toplevel :load-toplevel :execute)
285	(setq read-default-float-format initial-read-default-float-format))

Note: See TracBrowser for help on using the browser.

Download in other formats:

Original Format