source: project/release/4/geo-utils/trunk/geo-utils.scm @ 34506

Last change on this file since 34506 was 34506, checked in by kon, 14 months ago

add types, ? for preds

File size: 6.7 KB
Line 
1;;;; geo-utils.scm
2;;;; Kon Lovett, May '17
3;;;; Kon Lovett, Sep '17
4
5(module geo-utils
6
7(;export
8  ;
9  earth-flattening
10  earth-radius-miles
11  earth-radius-kilometers
12  ;
13  pythagorean-distance pythagorean-distance*
14  spherical-surface-distance
15  great-circle-distance great-circle-distance-radians
16  straight-line-distance
17  approximate-ellipsoid-distance
18  ;
19  great-circle-azimuth
20  ;
21  great-circle-position )
22
23(import scheme)
24
25(import chicken)
26
27(use numbers)
28
29(use fp-utils)
30
31(use geopoint)
32
33(include "geo-constants")
34
35;;
36
37(define earth-flattening EARTH-FLATTENING)
38
39(define earth-radius-miles EARTH-RADIUS-MILES)
40(define earth-radius-kilometers EARTH-RADIUS-KILOMETERS)
41
42;;
43
44(: pythagorean-distance (number number number number --> number))
45(define (pythagorean-distance lat1 lon1 lat2 lon2)
46  (sqrt (pythagorean-distance* lat1 lon1 lat2 lon2)) )
47
48(: pythagorean-distance* (number number number number --> number))
49(define (pythagorean-distance* lat1 lon1 lat2 lon2)
50  (let ((a (- lat1 lat2))
51        (b (- lon1 lon2)) )
52    (+ (* a a) (* b b)) ) )
53
54;;
55
56(: spherical-surface-distance (float float float float #!optional float --> float))
57(define (spherical-surface-distance lat1 lon1 lat2 lon2 #!optional (R EARTH-RADIUS-MILES))
58  ;haversine formula : https://en.wikipedia.org/wiki/Haversine_formula
59  (let ((dlat (fpdegree->radian (fpabs (fp- lat2 lat1))))
60        (dlon (fpdegree->radian (fpabs (fp- lon2 lon1))))
61        (lat1 (fpdegree->radian lat1))
62        (lat2 (fpdegree->radian lat2)) )
63    (let* ((a
64            (fp+
65              (fpsqr (fpsin (fp/ dlon 2.0)))
66              (fp*
67                (fp* (fpcos lat1) (fpcos lat2))
68                (fpsqr (fpsin (fp/ dlat 2.0))))) )
69           (c (fp* 2.0 (fpasin (fpmin 1.0 (fpsqrt a)))) ) )
70      (fp* c R) ) ) )
71
72;;
73
74(: great-circle-distance (float float float float #!optional float --> float))
75(define (great-circle-distance lat1 lon1 lat2 lon2 #!optional (R EARTH-RADIUS-MILES))
76  (let ((lat1 (fpdegree->radian lat1))
77        (lon1 (fpdegree->radian lon1))
78        (lat2 (fpdegree->radian lat2))
79        (lon2 (fpdegree->radian lon2)) )
80    (great-circle-distance-radians lat1 lon1 lat2 lon2 R) ) )
81
82; https://en.wikipedia.org/wiki/Great-circle_distance
83(: great-circle-distance-radians (float float float float #!optional float --> float))
84(define (great-circle-distance-radians lat1 lon1 lat2 lon2 #!optional (R EARTH-RADIUS-MILES))
85  (let ((d
86          (fp+
87            (fp* (fpsin lat1) (fpsin lat2))
88            (fp*
89              (fp* (fpcos lat1) (fpcos lat2))
90              (fpcos (fp- lon1 lon2)))) ) )
91    (fp* (fpacos d) R) ) )
92
93#| ;???
94  (atan2
95    (fp/
96      (sqrt ((sqr (cos lat2)) + (sqr ((cos lat1) • ((sin lon2) - (sin lon1)) • (cos (fpabs (lon2 - lon1)))))
97|#
98
99;;
100
101(: straight-line-distance (float float float float float float #!optional float --> float))
102(define (straight-line-distance lat1 lon1 h1 lat2 lon2 h2 #!optional (R EARTH-RADIUS-MILES))
103  ;filler
104  (pythagorean-distance lat1 lon1 lat2 lon2) )
105
106;;
107
108(: approximate-ellipsoid-distance (float float float float #!optional float --> float))
109(define (approximate-ellipsoid-distance lat1 lon1 lat2 lon2 #!optional (R EARTH-RADIUS-MILES))
110  (let ((lat1 (fpdegree->radian lat1))
111        (lon1 (fpdegree->radian lon1))
112        (lat2 (fpdegree->radian lat2))
113        (lon2 (fpdegree->radian lon2)) )
114    (let ((f (fp/ (fp+ lat1 lat2) 2.0))
115          (g (fp/ (fp- lat1 lat2) 2.0))
116          (l (fp/ (fp- lon1 lon2) 2.0)) )
117      (let ((sing (fpsin g))
118            (cosl (fpcos l))
119            (cosf (fpcos f))
120            (sinl (fpsin l))
121            (sinf (fpsin f))
122            (cosg (fpcos g)) )
123        (let ((sing2 (fpsqr sing))
124              (cosf2 (fpsqr cosf))
125              (cosg2 (fpsqr cosg))
126              (sinf2 (fpsqr sinf)) )
127          (let* ((s (fp+ (fp* sing2 (fpsqr cosl))
128                         (fp* cosf2 (fpsqr sinl))))
129                 (c (fp+ (fp* cosg2 (fpsqr cosl))
130                         (fp* sinf2 (fpsqr sinl))))
131                 (w (fpatan2 (fpsqrt s) (fpsqrt c)))
132                 (r (fpsqrt (fp/ (fp* s c) w)))
133                 (h1 (fp/ (fp* 3.0 (fp- r 1.0)) (fp* 2.0 c)))
134                 (h2 (fp/ (fp* 3.0 (fp+ r 1.0)) (fp* 2.0 s)))
135                 (d (fp* 2.0 (fp* w R))) )
136            (fp* d
137                 (fp+ 1.0
138                      (fp- (fp* (fp* h1 EARTH-FLATTENING)
139                                (fp* sinf2 cosg2))
140                           (fp* (fp* h2 EARTH-FLATTENING)
141                                (fp* cosf2 sing2))))) ) ) ) ) ) )
142
143;;
144
145(: great-circle-azimuth (float float float float #!optional (or fixnum float) --> float))
146(define (great-circle-azimuth lat1 lon1 lat2 lon2 #!optional (prec 5))
147  ;
148  (define precfact
149    (fp* 360.0 (fpprecision-factor prec)) )
150  ;
151  (define (clamp n)
152    (inexact->exact (fpround (fp* n precfact))) )
153  ;
154  (let ((ilat1 (clamp lat1))
155        (ilon1 (clamp lon1))
156        (ilat2 (clamp lat2))
157        (ilon2 (clamp lon2)) )
158    (cond
159      ((fx= ilon1 ilon2)
160        ; going up?
161        (if (fx> ilat1 ilat2)
162          180.0
163          ; going down or nowhere
164          0.0 ) )
165      (else
166        (let ((lat1 (fpdegree->radian lat1))
167              (lon1 (fpdegree->radian lon1))
168              (lat2 (fpdegree->radian lat2))
169              (lon2 (fpdegree->radian lon2)) )
170          (let* ((c (great-circle-distance-radians lat1 lon1 lat2 lon2))
171                 (a (fpasin (fp/ (fp* (fpcos lat2) (fpsin (fp- lon2 lon1)))
172                                 (fpsin c))))
173                 (r (fpradian->degree a)) )
174            (cond
175              ((fx> ilat2 ilat1)
176                (cond
177                  #; ;see else
178                  ((fx> ilon2 ilon1)
179                    r )
180                  ((fx< ilon2 ilon1)
181                    (fp+ r 360.0) )
182                  (else
183                    r ) ) )
184              ((fx< ilat2 ilat1)
185                (cond
186                  ((fx< ilon2 ilon1)
187                    (fp- 180.0 r) )
188                  ((fx> ilon2 ilon1)
189                    (fp- 180.0 r) )
190                  (else
191                    r ) ) )
192              (else
193                r ) ) ) ) ) ) ) )
194
195;;
196
197(: great-circle-position (float float float float #!optional float --> float float))
198(define (great-circle-position lat lon dis azi #!optional (R EARTH-RADIUS-MILES))
199  (let ((dlat (fpdegree->radian (fp- 90.0 lat)))
200        (lat (fpdegree->radian lat))
201        (lon (fpdegree->radian lon))
202        (azi (fpdegree->radian azi)) )
203    (let* ((b (fp/ dis R))
204           (sinb (fpsin b))
205           (a (fpacos (fp+ (fp* (fpcos b) (fpcos dlat))
206                           (fp* sinb (fp* (fpsin dlat) (fpcos azi))))))
207           (b (fpasin (fp/ (fp* sinb (fpsin azi)) (fpsin a)))) )
208      (values (fp- 90.0 (fpradian->degree a)) (fpradian->degree (fp+ b lon))) ) ) )
209
210) ;geo-utils
Note: See TracBrowser for help on using the repository browser.