-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathsph_subs.py
419 lines (350 loc) · 11.8 KB
/
sph_subs.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
"""
sphere_subs.py
included in Package sphere/
/Dropbox/Dropbox/gprieto/python/Modules/sphere
SPHERE_SUBS is set of Python function definitions to compute distances
and angles on a spherical Earth. All angles are in degrees.
Latitude used on standard maps is geographic latitude; this may be
converted to the geocentric latitude used in these routines
by using the SPH_GEOCENTRIC function. Longitude input to
these routines may be from either -180 to 180 or from 0 to 360.
Longitude returned from these routines will be from 0 to 360.
Based on Fortran subroutines from Peter M. Shearer
Last Modified
German A. Prieto
May 2017
Contains definitions of
sph_loc - finds location of second point on sphere, given range
and azimuth at first point.
sph_dist - computes angular separation of two points on sphere
sph_azi - computes distance and azimuth between two points
on the sphere
sph_mid - finds midpoint between two surface points on sphere
and azimuth at midpoint to second point
"""
#--------------------------------------------------------------------
# sph_loc
def sph_loc(flat1,flon1,delta,azi):
"""
SPH_LOC finds location of second point on sphere, given range
and azimuth at first point.
Inputs: flat1 = latitude of first point (degrees)
flon1 = longitude of first point (degrees)
del = angular separation between points (degrees)
azi = azimuth at 1st point to 2nd point, from N (deg.)
Returns: flat2 = latitude of second point (degrees)
flon2 = longitude of second point (degrees)
"""
# import appropriate functions
import math as mt
import numpy as np
if (delta == 0):
flat2 = flat1
flon2 = flon1
return flat2, flon2
raddeg = mt.pi/180.
delr = delta*raddeg
azr = azi*raddeg
theta1=(90.-flat1)*raddeg
phi1=flon1*raddeg
ctheta2 = mt.sin(delr)*mt.sin(theta1)*mt.cos(azr) + mt.cos(theta1)*mt.cos(delr)
theta2 = mt.acos(ctheta2)
if (theta1 == 0.):
phi2 = azr
elif (theta2 == 0.):
phi2=0.0
else:
sphi2 = mt.sin(delr)*mt.sin(azr)/mt.sin(theta2)
cphi2 = (mt.cos(delr)-mt.cos(theta1)*ctheta2)/(mt.sin(theta1)*mt.sin(theta2))
phi2 = phi1+mt.atan2(sphi2,cphi2)
flat2=90.-theta2/raddeg
flon2=phi2/raddeg
if (flon2 > 360.):
flon2=flon2-360.
if (flon2 < 0.):
flon2=flon2+360.
return flat2, flon2
# end sph_loc
#--------------------------------------------------------------------
#--------------------------------------------------------------------
# sph_dist
def sph_dist(flat1,flon1,flat2,flon2):
"""
SPH_DIST computes angular separation of two points on sphere
Inputs: flat1 = latitude of first point (degrees)
flon1 = longitude of first point (degrees)
flat2 = latitude of second point (degrees)
flon2 = longitude of second point (degrees)
Returns: delta = angular separation between points (degrees)
Note: This routine is accurate depending on the precision of the
real numbers used. Python should be accurate to real(8) precision
For greater accuracy, perform a separate calculation for close
ranges using Cartesian geometry.
"""
# import appropriate functions
import math as mt
import numpy as np
if ( (flat1 == flat2 and flon1 == flon2)
or (flat1 == 90. and flat2 == 90.)
or (flat1 == -90. and flat2 == -90.) ):
delta = 0.
return delta
raddeg=mt.pi/180.
theta1=(90.-flat1)*raddeg
theta2=(90.-flat2)*raddeg
phi1=flon1*raddeg
phi2=flon2*raddeg
stheta1=mt.sin(theta1)
stheta2=mt.sin(theta2)
ctheta1=mt.cos(theta1)
ctheta2=mt.cos(theta2)
cang=stheta1*stheta2*mt.cos(phi2-phi1)+ctheta1*ctheta2
ang=mt.acos(cang)
delta=ang/raddeg
return delta
# end sph_dist
#--------------------------------------------------------------------
#--------------------------------------------------------------------
# sph_azi
def sph_azi(flat1,flon1,flat2,flon2):
"""
SPH_AZI computes distance and azimuth between two points
on the sphere
Inputs: flat1 = latitude of first point (degrees)
flon2 = longitude of first point (degrees)
flat2 = latitude of second point (degrees)
flon2 = longitude of second point (degrees)
Returns: del = angular separation between points (degrees)
azi = azimuth at 1st point to 2nd point, from N (deg.)
Notes:
(1) applies to geocentric not geographic lat,lon on Earth
(2) This routine is accurate depending on the precision of the
real numbers used. Python should be accurate to real(8) precision
For greater accuracy, perform a separate calculation for close
ranges using Cartesian geometry.
"""
# import appropriate functions
import math as mt
import numpy as np
if ( (flat1 == flat2 and flon1 == flon2)
or (flat1 == 90. and flat2 == 90.)
or (flat1 == -90. and flat2 == -90.) ):
delta = 0.
azi = 0.
return [delta,azi]
# Perform calculation
raddeg=mt.pi/180.
theta1=(90.-flat1)*raddeg
theta2=(90.-flat2)*raddeg
phi1=flon1*raddeg
phi2=flon2*raddeg
stheta1=mt.sin(theta1)
stheta2=mt.sin(theta2)
ctheta1=mt.cos(theta1)
ctheta2=mt.cos(theta2)
cang=stheta1*stheta2*mt.cos(phi2-phi1)+ctheta1*ctheta2
ang=mt.acos(cang)
delta=ang/raddeg
sang=mt.sqrt(1.-cang*cang)
caz=(ctheta2-ctheta1*cang)/(sang*stheta1)
saz=-stheta2*mt.sin(phi1-phi2)/sang
az=mt.atan2(saz,caz)
azi=az/raddeg
if (azi < 0.):
azi=azi+360.
return [delta, azi]
# end sph_azi
#--------------------------------------------------------------------
#--------------------------------------------------------------------
# sph_mid
def sph_mid(flat1,flon1,flat2,flon2):
"""
SPH_MID finds midpoint between two surface points on sphere
and azimuth at midpoint to second point
Requires: SPH_AZI, SPH_LOC
Inputs: flat1 = latitude of first point (degrees)
flon1 = longitude of first point (degrees)
flat2 = latitude of second point (degrees)
flon2 = longitude of second point (degrees)
Returns: delta = angular separation between points (degrees)
flat3 = latitude of midpoint (degrees)
flon3 = longitude of midpoint (degrees)
azi = azimuth at midpoint to first point
calls sph_azi, sph_loc
"""
delta,azi0 = sph_azi(flat1,flon1,flat2,flon2)
del0 = delta/2.
flat3,flon3 = sph_loc(flat1,flon1,del0,azi0)
del2,azi = sph_azi(flat3,flon3,flat2,flon2)
return delta,flat3,flon3,azi
# end sph_mid
#--------------------------------------------------------------------
#***************************
# Codes to work on
#***************************
#c SPH_GEOCENTRIC converts Earth's geographic latitude to geocentric latitude
#c
#c Inputs: flat1 = geographic latitude (degrees)
#c (relative to local vertical)
#c Returns: flat2 = geocentric latitude (degrees)
#c (relative to Earth's center)
#c
# subroutine SPH_GEOCENTRIC(flat1,flat2)
# if (flat1.eq.90.) then
# flat2=90.
# return
# end if
# pi=3.141592654
# degrad=180./pi
# factor=0.9933056 !=(1-1/298.256)**2
# phi=flat1/degrad
# theta=atan(factor*tan(phi))
# flat2=theta*degrad
# return
# end
#3
#
#c SPH_GEOGRAPHIC converts Earth's geocentric latitude to geographic latitude
#c
#c Inputs: flat1 = geocentric latitude (degrees)
#c (relative to Earth's center)
#c Returns: flat2 = geographic latitude (degrees)
#c (relative to local vertical)
#c
# subroutine SPH_GEOGRAPHIC(flat1,flat2)
# if (flat1.eq.90.) then
# flat2=90.
# return
# end if
# pi=3.141592654
# degrad=180./pi
# factor=0.9933056 !=(1-1/298.256)**2
# theta=flat1/degrad
# phi=atan(tan(theta)/factor)
# flat2=phi*degrad
# return
# end
#
#
#c SPH_CART converts from spherical (lat,lon,r) to cartesian (x,y,z)
#c
#c Inputs: flat = latitude (degrees)
#c flon = longitude (degrees)
#c r = radius
#c Returns: x,y,z = x,y,z coordinates (z=up)
#c
# subroutine SPH_CART(flat,flon,r,x,y,z)
# degrad=180./3.1415927
# theta=(90.-flat)/degrad
# phi=flon/degrad
# stheta=sin(theta)
# x=stheta*cos(phi)*r
# y=stheta*sin(phi)*r
# z=cos(theta)*r
# return
# end
#
#
#c SPH_SPH converts from cartesian (x,y,z) to spherical (lat,lon,r)
#c
#c Inputs: x,y,z = x,y,z coordinates (z=up)
#c Returns: flat = latitude (degrees)
#c flon = longitude (degrees)
#c r = radius
#c
# subroutine SPH_SPH(x,y,z,flat,flon,r)
# degrad=180./3.1415927
# r=sqrt(x*x+y*y+z*z)
# if (r.eq.0.) then
# flat=0.
# flon=0.
# return
# end if
# d=sqrt(x*x+y*y)
# theta=atan2(d,z)
# phi=atan2(y,x)
# flat=90.-theta*degrad
# flon=phi*degrad
# if (flon.lt.0.) flon=flon+360.
# return
# end
#
#
#c SPH_POLE finds pole of great circle joining two points on sphere
#c
#c Requires: SPH_CART, SPH_SPH
#c
#c Inputs: flat1 = latitude of first point (degrees)
#c flon2 = longitude of first point (degrees)
#c flat2 = latitude of second point (degrees)
#c flon2 = longitude of second point (degrees)
#c Returns: flat3 = latitude of great circle pole (degrees)
#c flon3 = longitude of great circle pole (degrees)
#c
# subroutine SPH_POLE(flat1,flon1,flat2,flon2,flat3,flon3)
# call SPH_CART(flat1,flon1,1.,x1,y1,z1)
# call SPH_CART(flat2,flon2,1.,x2,y2,z2)
# x3=y1*z2-y2*z1
# y3=x2*z1-x1*z2
# z3=x1*y2-x2*y1
# call SPH_SPH(x3,y3,z3,flat3,flon3,r)
# return
# end
#
#
#c SPH_TO_GCIRC rotates (lat,lon) to new coordinates using
#c specified great circle as new equator
#c
#c Requires: SPH_CART, SPH_SPH
#c
#c Inputs: flat = latitude of point (degrees)
#c flon = longitude of point (degrees)
#c plat = latitude of great circle pole (degrees)
#c plon = longitude of great circle pole (degrees)
#c Returns: glat = latitude of point relative to great circle (degrees)
#c glon = longitude of point relative to great circle (degrees)
#c
# subroutine SPH_TO_GCIRC(flat,flon,plat,plon,glat,glon)
# degrad=180./3.1415927
# flon1=flon-plon
# if (flon1.lt.0.) flon1=flon1+360.
# call SPH_CART(flat,flon1,1.,x,y,z)
# pcolat=90.-plat
# the=-pcolat/degrad
# costhe=cos(the)
# sinthe=sin(the)
# yr=y
# xr=costhe*x+sinthe*z
# zr=costhe*z-sinthe*x
# call SPH_SPH(xr,yr,zr,glat,glon,r)
# return
# end
#
#
#c SPH_FROM_GCIRC rotates (lat,lon) from great circle coordinates
#c back to original coordinates
#c
#c Requires: SPH_CART, SPH_SPH
#c
#c Inputs: glat = latitude of point relative to great circle (degrees)
#c glon = longitude of point relative to great circle (degrees)
#c plat = latitude of great circle pole (degrees)
#c plon = longitude of great circle pole (degrees)
#c Returns: flat = latitude of point (degrees)
#c flon = longitude of point (degrees)
#c
# subroutine SPH_FROM_GCIRC(glat,glon,plat,plon,flat,flon)
# degrad=180./3.1415927
# call SPH_CART(glat,glon,1.,xr,yr,zr)
# pcolat=90.-plat
# the=-pcolat/degrad
# costhe=cos(the)
# sinthe=sin(the)
# y=yr
# x=costhe*xr-sinthe*zr
# z=costhe*zr+sinthe*xr
# call SPH_SPH(x,y,z,flat,flon1,r)
# flon=flon1+plon
# if (flon.gt.360.) flon=flon-360.
# return
# end