Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/orthodrome.py: 75%
427 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-03-07 11:54 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7Some basic geodetic functions.
8'''
10import math
11import numpy as num
13from .moment_tensor import euler_to_matrix
14from .config import config
15from .plot.beachball import spoly_cut
17from matplotlib.path import Path
19d2r = math.pi/180.
20r2d = 1./d2r
21earth_oblateness = 1./298.257223563
22earthradius_equator = 6378.14 * 1000.
23earthradius = config().earthradius
24d2m = earthradius_equator*math.pi/180.
25m2d = 1./d2m
27_testpath = Path([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)], closed=True)
29raise_if_slow_path_contains_points = False
32class Slow(Exception):
33 pass
36if hasattr(_testpath, 'contains_points') and num.all(
37 _testpath.contains_points([(0.5, 0.5), (1.5, 0.5)]) == [True, False]):
39 def path_contains_points(verts, points):
40 p = Path(verts, closed=True)
41 return p.contains_points(points).astype(bool)
43else:
44 # work around missing contains_points and bug in matplotlib ~ v1.2.0
46 def path_contains_points(verts, points):
47 if raise_if_slow_path_contains_points:
48 # used by unit test to skip slow gshhg_example.py
49 raise Slow()
51 p = Path(verts, closed=True)
52 result = num.zeros(points.shape[0], dtype=bool)
53 for i in range(result.size):
54 result[i] = p.contains_point(points[i, :])
56 return result
59try:
60 cbrt = num.cbrt
61except AttributeError:
62 def cbrt(x):
63 return x**(1./3.)
66def float_array_broadcast(*args):
67 return num.broadcast_arrays(*[
68 num.asarray(x, dtype=float) for x in args])
71class Loc(object):
72 '''
73 Simple location representation.
75 :attrib lat: Latitude in [deg].
76 :attrib lon: Longitude in [deg].
77 '''
78 def __init__(self, lat, lon):
79 self.lat = lat
80 self.lon = lon
83def clip(x, mi, ma):
84 '''
85 Clip values of an array.
87 :param x: Continunous data to be clipped.
88 :param mi: Clip minimum.
89 :param ma: Clip maximum.
90 :type x: :py:class:`numpy.ndarray`
91 :type mi: float
92 :type ma: float
94 :return: Clipped data.
95 :rtype: :py:class:`numpy.ndarray`
96 '''
97 return num.minimum(num.maximum(mi, x), ma)
100def wrap(x, mi, ma):
101 '''
102 Wrapping continuous data to fundamental phase values.
104 .. math::
105 x_{\\mathrm{wrapped}} = x_{\\mathrm{cont},i} -
106 \\frac{ x_{\\mathrm{cont},i} - r_{\\mathrm{min}} }
107 { r_{\\mathrm{max}} - r_{\\mathrm{min}}}
108 \\cdot ( r_{\\mathrm{max}} - r_{\\mathrm{min}}),\\quad
109 x_{\\mathrm{wrapped}}\\; \\in
110 \\;[ r_{\\mathrm{min}},\\, r_{\\mathrm{max}}].
112 :param x: Continunous data to be wrapped.
113 :param mi: Minimum value of wrapped data.
114 :param ma: Maximum value of wrapped data.
115 :type x: :py:class:`numpy.ndarray`
116 :type mi: float
117 :type ma: float
119 :return: Wrapped data.
120 :rtype: :py:class:`numpy.ndarray`
121 '''
122 return x - num.floor((x-mi)/(ma-mi)) * (ma-mi)
125def _latlon_pair(args):
126 if len(args) == 2:
127 a, b = args
128 return a.lat, a.lon, b.lat, b.lon
130 elif len(args) == 4:
131 return args
134def cosdelta(*args):
135 '''
136 Cosine of the angular distance between two points ``a`` and ``b`` on a
137 sphere.
139 This function (find implementation below) returns the cosine of the
140 distance angle 'delta' between two points ``a`` and ``b``, coordinates of
141 which are expected to be given in geographical coordinates and in degrees.
142 For numerical stability a maximum of 1.0 is enforced.
144 .. math::
146 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
147 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
148 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
149 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\[0.5cm]
151 \\cos(\\Delta) = \\min( 1.0, \\quad \\sin( A_{\\mathrm{lat'}})
152 \\sin( B_{\\mathrm{lat'}} ) +
153 \\cos(A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
154 \\cos( B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )
156 :param a: Location point A.
157 :type a: :py:class:`pyrocko.orthodrome.Loc`
158 :param b: Location point B.
159 :type b: :py:class:`pyrocko.orthodrome.Loc`
161 :return: Cosdelta.
162 :rtype: float
163 '''
165 alat, alon, blat, blon = _latlon_pair(args)
167 return min(
168 1.0,
169 math.sin(alat*d2r) * math.sin(blat*d2r) +
170 math.cos(alat*d2r) * math.cos(blat*d2r) *
171 math.cos(d2r*(blon-alon)))
174def cosdelta_numpy(a_lats, a_lons, b_lats, b_lons):
175 '''
176 Cosine of the angular distance between two points ``a`` and ``b`` on a
177 sphere.
179 This function returns the cosines of the distance
180 angles *delta* between two points ``a`` and ``b`` given as
181 :py:class:`numpy.ndarray`.
182 The coordinates are expected to be given in geographical coordinates
183 and in degrees. For numerical stability a maximum of ``1.0`` is enforced.
185 Please find the details of the implementation in the documentation of
186 the function :py:func:`pyrocko.orthodrome.cosdelta` above.
188 :param a_lats: Latitudes in [deg] point A.
189 :param a_lons: Longitudes in [deg] point A.
190 :param b_lats: Latitudes in [deg] point B.
191 :param b_lons: Longitudes in [deg] point B.
192 :type a_lats: :py:class:`numpy.ndarray`
193 :type a_lons: :py:class:`numpy.ndarray`
194 :type b_lats: :py:class:`numpy.ndarray`
195 :type b_lons: :py:class:`numpy.ndarray`
197 :return: Cosdelta.
198 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
199 '''
200 return num.minimum(
201 1.0,
202 num.sin(a_lats*d2r) * num.sin(b_lats*d2r) +
203 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
204 num.cos(d2r*(b_lons-a_lons)))
207def azimuth(*args):
208 '''
209 Azimuth calculation.
211 This function (find implementation below) returns azimuth ...
212 between points ``a`` and ``b``, coordinates of
213 which are expected to be given in geographical coordinates and in degrees.
215 .. math::
217 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
218 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
219 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
220 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\
222 \\varphi_{\\mathrm{azi},AB} = \\frac{180}{\\pi} \\arctan \\left[
223 \\frac{
224 \\cos( A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
225 \\sin(B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )}
226 {\\sin ( B_{\\mathrm{lat'}} ) - \\sin( A_{\\mathrm{lat'}}
227 cosdelta) } \\right]
229 :param a: Location point A.
230 :type a: :py:class:`pyrocko.orthodrome.Loc`
231 :param b: Location point B.
232 :type b: :py:class:`pyrocko.orthodrome.Loc`
234 :return: Azimuth in degree
235 '''
237 alat, alon, blat, blon = _latlon_pair(args)
239 return r2d*math.atan2(
240 math.cos(alat*d2r) * math.cos(blat*d2r) *
241 math.sin(d2r*(blon-alon)),
242 math.sin(d2r*blat) - math.sin(d2r*alat) * cosdelta(
243 alat, alon, blat, blon))
246def azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta=None):
247 '''
248 Calculation of the azimuth (*track angle*) from a location A towards B.
250 This function returns azimuths (*track angles*) from locations A towards B
251 given in :py:class:`numpy.ndarray`. Coordinates are expected to be given in
252 geographical coordinates and in degrees.
254 Please find the details of the implementation in the documentation of the
255 function :py:func:`pyrocko.orthodrome.azimuth`.
258 :param a_lats: Latitudes in [deg] point A.
259 :param a_lons: Longitudes in [deg] point A.
260 :param b_lats: Latitudes in [deg] point B.
261 :param b_lons: Longitudes in [deg] point B.
262 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
263 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
264 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
265 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
267 :return: Azimuths in [deg].
268 :rtype: :py:class:`numpy.ndarray`, ``(N)``
269 '''
270 if _cosdelta is None:
271 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
273 return r2d*num.arctan2(
274 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
275 num.sin(d2r*(b_lons-a_lons)),
276 num.sin(d2r*b_lats) - num.sin(d2r*a_lats) * _cosdelta)
279def azibazi(*args, **kwargs):
280 '''
281 Azimuth and backazimuth from location A towards B and back.
283 :returns: Azimuth in [deg] from A to B, back azimuth in [deg] from B to A.
284 :rtype: tuple[float, float]
285 '''
287 alat, alon, blat, blon = _latlon_pair(args)
288 if alat == blat and alon == blon:
289 return 0., 180.
291 implementation = kwargs.get('implementation', 'c')
292 assert implementation in ('c', 'python')
293 if implementation == 'c':
294 from pyrocko import orthodrome_ext
295 return orthodrome_ext.azibazi(alat, alon, blat, blon)
297 cd = cosdelta(alat, alon, blat, blon)
298 azi = r2d*math.atan2(
299 math.cos(alat*d2r) * math.cos(blat*d2r) *
300 math.sin(d2r*(blon-alon)),
301 math.sin(d2r*blat) - math.sin(d2r*alat) * cd)
302 bazi = r2d*math.atan2(
303 math.cos(blat*d2r) * math.cos(alat*d2r) *
304 math.sin(d2r*(alon-blon)),
305 math.sin(d2r*alat) - math.sin(d2r*blat) * cd)
307 return azi, bazi
310def azibazi_numpy(a_lats, a_lons, b_lats, b_lons, implementation='c'):
311 '''
312 Azimuth and backazimuth from location A towards B and back.
314 Arguments are given as :py:class:`numpy.ndarray`.
316 :param a_lats: Latitude(s) in [deg] of point A.
317 :type a_lats: :py:class:`numpy.ndarray`
318 :param a_lons: Longitude(s) in [deg] of point A.
319 :type a_lons: :py:class:`numpy.ndarray`
320 :param b_lats: Latitude(s) in [deg] of point B.
321 :type b_lats: :py:class:`numpy.ndarray`
322 :param b_lons: Longitude(s) in [deg] of point B.
323 :type b_lons: :py:class:`numpy.ndarray`
325 :returns: Azimuth(s) in [deg] from A to B,
326 back azimuth(s) in [deg] from B to A.
327 :rtype: :py:class:`numpy.ndarray`, :py:class:`numpy.ndarray`
328 '''
330 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
331 a_lats, a_lons, b_lats, b_lons)
333 assert implementation in ('c', 'python')
334 if implementation == 'c':
335 from pyrocko import orthodrome_ext
336 return orthodrome_ext.azibazi_numpy(a_lats, a_lons, b_lats, b_lons)
338 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
339 azis = azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta)
340 bazis = azimuth_numpy(b_lats, b_lons, a_lats, a_lons, _cosdelta)
342 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
343 ii_eq = num.where(eq)[0]
344 azis[ii_eq] = 0.0
345 bazis[ii_eq] = 180.0
346 return azis, bazis
349def azidist_numpy(*args):
350 '''
351 Calculation of the azimuth (*track angle*) and the distance from locations
352 A towards B on a sphere.
354 The assisting functions used are :py:func:`pyrocko.orthodrome.cosdelta` and
355 :py:func:`pyrocko.orthodrome.azimuth`
357 :param a_lats: Latitudes in [deg] point A.
358 :param a_lons: Longitudes in [deg] point A.
359 :param b_lats: Latitudes in [deg] point B.
360 :param b_lons: Longitudes in [deg] point B.
361 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
362 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
363 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
364 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
366 :return: Azimuths in [deg], distances in [deg].
367 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
368 '''
369 _cosdelta = cosdelta_numpy(*args)
370 _azimuths = azimuth_numpy(_cosdelta=_cosdelta, *args)
371 return _azimuths, r2d*num.arccos(_cosdelta)
374def distance_accurate50m(*args, **kwargs):
375 '''
376 Accurate distance calculation based on a spheroid of rotation.
378 Function returns distance in meter between points A and B, coordinates of
379 which must be given in geographical coordinates and in degrees.
380 The returned distance should be accurate to 50 m using WGS84.
381 Values for the Earth's equator radius and the Earth's oblateness
382 (``f_oblate``) are defined in the pyrocko configuration file
383 :py:class:`pyrocko.config`.
385 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
387 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
388 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
390 .. math::
392 F = \\frac{\\pi}{180}
393 \\frac{(A_{lat} + B_{lat})}{2}, \\quad
394 G = \\frac{\\pi}{180}
395 \\frac{(A_{lat} - B_{lat})}{2}, \\quad
396 l = \\frac{\\pi}{180}
397 \\frac{(A_{lon} - B_{lon})}{2} \\quad
398 \\\\[0.5cm]
399 S = \\sin^2(G) \\cdot \\cos^2(l) +
400 \\cos^2(F) \\cdot \\sin^2(l), \\quad \\quad
401 C = \\cos^2(G) \\cdot \\cos^2(l) +
402 \\sin^2(F) \\cdot \\sin^2(l)
404 .. math::
406 w = \\arctan \\left( \\sqrt{ \\frac{S}{C}} \\right) , \\quad
407 r = \\sqrt{\\frac{S}{C} }
409 The spherical-earth distance D between A and B, can be given with:
411 .. math::
413 D_{sphere} = 2w \\cdot R_{equator}
415 The oblateness of the Earth requires some correction with
416 correction factors h1 and h2:
418 .. math::
420 h_1 = \\frac{3r - 1}{2C}, \\quad
421 h_2 = \\frac{3r +1 }{2S}\\\\[0.5cm]
423 D = D_{\\mathrm{sphere}} \\cdot [ 1 + h_1 \\,f_{\\mathrm{oblate}}
424 \\cdot \\sin^2(F)
425 \\cos^2(G) - h_2\\, f_{\\mathrm{oblate}}
426 \\cdot \\cos^2(F) \\sin^2(G)]
429 :param a: Location point A.
430 :type a: :py:class:`pyrocko.orthodrome.Loc`
431 :param b: Location point B.
432 :type b: :py:class:`pyrocko.orthodrome.Loc`
434 :return: Distance in [m].
435 :rtype: float
436 '''
438 alat, alon, blat, blon = _latlon_pair(args)
440 implementation = kwargs.get('implementation', 'c')
441 assert implementation in ('c', 'python')
442 if implementation == 'c':
443 from pyrocko import orthodrome_ext
444 return orthodrome_ext.distance_accurate50m(alat, alon, blat, blon)
446 f = (alat + blat)*d2r / 2.
447 g = (alat - blat)*d2r / 2.
448 h = (alon - blon)*d2r / 2.
450 s = math.sin(g)**2 * math.cos(h)**2 + math.cos(f)**2 * math.sin(h)**2
451 c = math.cos(g)**2 * math.cos(h)**2 + math.sin(f)**2 * math.sin(h)**2
453 w = math.atan(math.sqrt(s/c))
455 if w == 0.0:
456 return 0.0
458 r = math.sqrt(s*c)/w
459 d = 2.*w*earthradius_equator
460 h1 = (3.*r-1.)/(2.*c)
461 h2 = (3.*r+1.)/(2.*s)
463 return d * (1. +
464 earth_oblateness * h1 * math.sin(f)**2 * math.cos(g)**2 -
465 earth_oblateness * h2 * math.cos(f)**2 * math.sin(g)**2)
468def distance_accurate50m_numpy(
469 a_lats, a_lons, b_lats, b_lons, implementation='c'):
471 '''
472 Accurate distance calculation based on a spheroid of rotation.
474 Function returns distance in meter between points ``a`` and ``b``,
475 coordinates of which must be given in geographical coordinates and in
476 degrees.
477 The returned distance should be accurate to 50 m using WGS84.
478 Values for the Earth's equator radius and the Earth's oblateness
479 (``f_oblate``) are defined in the pyrocko configuration file
480 :py:class:`pyrocko.config`.
482 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
484 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
485 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
487 .. math::
489 F_i = \\frac{\\pi}{180}
490 \\frac{(a_{lat,i} + a_{lat,i})}{2}, \\quad
491 G_i = \\frac{\\pi}{180}
492 \\frac{(a_{lat,i} - b_{lat,i})}{2}, \\quad
493 l_i= \\frac{\\pi}{180}
494 \\frac{(a_{lon,i} - b_{lon,i})}{2} \\\\[0.5cm]
495 S_i = \\sin^2(G_i) \\cdot \\cos^2(l_i) +
496 \\cos^2(F_i) \\cdot \\sin^2(l_i), \\quad \\quad
497 C_i = \\cos^2(G_i) \\cdot \\cos^2(l_i) +
498 \\sin^2(F_i) \\cdot \\sin^2(l_i)
500 .. math::
502 w_i = \\arctan \\left( \\sqrt{\\frac{S_i}{C_i}} \\right), \\quad
503 r_i = \\sqrt{\\frac{S_i}{C_i} }
505 The spherical-earth distance ``D`` between ``a`` and ``b``,
506 can be given with:
508 .. math::
510 D_{\\mathrm{sphere},i} = 2w_i \\cdot R_{\\mathrm{equator}}
512 The oblateness of the Earth requires some correction with
513 correction factors ``h1`` and ``h2``:
515 .. math::
517 h_{1.i} = \\frac{3r - 1}{2C_i}, \\quad
518 h_{2,i} = \\frac{3r +1 }{2S_i}\\\\[0.5cm]
520 D_{AB,i} = D_{\\mathrm{sphere},i} \\cdot [1 + h_{1,i}
521 \\,f_{\\mathrm{oblate}}
522 \\cdot \\sin^2(F_i)
523 \\cos^2(G_i) - h_{2,i}\\, f_{\\mathrm{oblate}}
524 \\cdot \\cos^2(F_i) \\sin^2(G_i)]
526 :param a_lats: Latitudes in [deg] point A.
527 :param a_lons: Longitudes in [deg] point A.
528 :param b_lats: Latitudes in [deg] point B.
529 :param b_lons: Longitudes in [deg] point B.
530 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
531 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
532 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
533 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
535 :return: Distances in [m].
536 :rtype: :py:class:`numpy.ndarray`, ``(N)``
537 '''
539 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
540 a_lats, a_lons, b_lats, b_lons)
542 assert implementation in ('c', 'python')
543 if implementation == 'c':
544 from pyrocko import orthodrome_ext
545 return orthodrome_ext.distance_accurate50m_numpy(
546 a_lats, a_lons, b_lats, b_lons)
548 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
549 ii_neq = num.where(num.logical_not(eq))[0]
551 if num.all(eq):
552 return num.zeros_like(eq, dtype=float)
554 def extr(x):
555 if isinstance(x, num.ndarray) and x.size > 1:
556 return x[ii_neq]
557 else:
558 return x
560 a_lats = extr(a_lats)
561 a_lons = extr(a_lons)
562 b_lats = extr(b_lats)
563 b_lons = extr(b_lons)
565 f = (a_lats + b_lats)*d2r / 2.
566 g = (a_lats - b_lats)*d2r / 2.
567 h = (a_lons - b_lons)*d2r / 2.
569 s = num.sin(g)**2 * num.cos(h)**2 + num.cos(f)**2 * num.sin(h)**2
570 c = num.cos(g)**2 * num.cos(h)**2 + num.sin(f)**2 * num.sin(h)**2
572 w = num.arctan(num.sqrt(s/c))
574 r = num.sqrt(s*c)/w
576 d = 2.*w*earthradius_equator
577 h1 = (3.*r-1.)/(2.*c)
578 h2 = (3.*r+1.)/(2.*s)
580 dists = num.zeros(eq.size, dtype=float)
581 dists[ii_neq] = d * (
582 1. +
583 earth_oblateness * h1 * num.sin(f)**2 * num.cos(g)**2 -
584 earth_oblateness * h2 * num.cos(f)**2 * num.sin(g)**2)
586 return dists
589def ne_to_latlon(lat0, lon0, north_m, east_m):
590 '''
591 Transform local cartesian coordinates to latitude and longitude.
593 From east and north coordinates (``x`` and ``y`` coordinate
594 :py:class:`numpy.ndarray`) relative to a reference differences in
595 longitude and latitude are calculated, which are effectively changes in
596 azimuth and distance, respectively:
598 .. math::
600 \\text{distance change:}\\; \\Delta {\\bf{a}} &= \\sqrt{{\\bf{y}}^2 +
601 {\\bf{x}}^2 }/ \\mathrm{R_E},
603 \\text{azimuth change:}\\; \\Delta \\bf{\\gamma} &= \\arctan( \\bf{x}
604 / \\bf{y}).
606 The projection used preserves the azimuths of the input points.
608 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
609 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
610 :param north_m: Northing distances from origin in [m].
611 :param east_m: Easting distances from origin in [m].
612 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
613 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
614 :type lat0: float
615 :type lon0: float
617 :return: Array with latitudes and longitudes in [deg].
618 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
620 '''
622 a = num.sqrt(north_m**2+east_m**2)/earthradius
623 gamma = num.arctan2(east_m, north_m)
625 return azidist_to_latlon_rad(lat0, lon0, gamma, a)
628def azidist_to_latlon(lat0, lon0, azimuth_deg, distance_deg):
629 '''
630 Absolute latitudes and longitudes are calculated from relative changes.
632 Convenience wrapper to :py:func:`azidist_to_latlon_rad` with azimuth and
633 distance given in degrees.
635 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
636 :type lat0: float
637 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
638 :type lon0: float
639 :param azimuth_deg: Azimuth from origin in [deg].
640 :type azimuth_deg: :py:class:`numpy.ndarray`, ``(N)``
641 :param distance_deg: Distances from origin in [deg].
642 :type distance_deg: :py:class:`numpy.ndarray`, ``(N)``
644 :return: Array with latitudes and longitudes in [deg].
645 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
646 '''
648 return azidist_to_latlon_rad(
649 lat0, lon0, azimuth_deg/180.*num.pi, distance_deg/180.*num.pi)
652def azidist_to_latlon_rad(lat0, lon0, azimuth_rad, distance_rad):
653 '''
654 Absolute latitudes and longitudes are calculated from relative changes.
656 For numerical stability a range between of ``-1.0`` and ``1.0`` is
657 enforced for ``c`` and ``alpha``.
659 .. math::
661 \\Delta {\\bf a}_i \\; \\text{and} \\; \\Delta \\gamma_i \\;
662 \\text{are relative distances and azimuths from lat0 and lon0 for
663 \\textit{i} source points of a finite source.}
665 .. math::
667 \\mathrm{b} &= \\frac{\\pi}{2} -\\frac{\\pi}{180}\\;\\mathrm{lat_0}\\\\
668 {\\bf c}_i &=\\arccos[\\; \\cos(\\Delta {\\bf{a}}_i)
669 \\cos(\\mathrm{b}) + |\\Delta \\gamma_i| \\,
670 \\sin(\\Delta {\\bf a}_i)
671 \\sin(\\mathrm{b})\\; ] \\\\
672 \\mathrm{lat}_i &= \\frac{180}{\\pi}
673 \\left(\\frac{\\pi}{2} - {\\bf c}_i \\right)
675 .. math::
677 \\alpha_i &= \\arcsin \\left[ \\; \\frac{ \\sin(\\Delta {\\bf a}_i )
678 \\sin(|\\Delta \\gamma_i|)}{\\sin({\\bf c}_i)}\\;
679 \\right] \\\\
680 \\alpha_i &= \\begin{cases}
681 \\alpha_i, &\\text{if} \\; \\cos(\\Delta {\\bf a}_i) -
682 \\cos(\\mathrm{b}) \\cos({\\bf{c}}_i) > 0, \\;
683 \\text{else} \\\\
684 \\pi - \\alpha_i, & \\text{if} \\; \\alpha_i > 0,\\;
685 \\text{else}\\\\
686 -\\pi - \\alpha_i, & \\text{if} \\; \\alpha_i < 0.
687 \\end{cases} \\\\
688 \\mathrm{lon}_i &= \\mathrm{lon_0} +
689 \\frac{180}{\\pi} \\,
690 \\frac{\\Delta \\gamma_i }{|\\Delta \\gamma_i|}
691 \\cdot \\alpha_i
692 \\text{, with $\\alpha_i \\in [-\\pi,\\pi]$}
694 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
695 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
696 :param distance_rad: Distances from origin in [rad].
697 :param azimuth_rad: Azimuth from origin in [rad].
698 :type distance_rad: :py:class:`numpy.ndarray`, ``(N)``
699 :type azimuth_rad: :py:class:`numpy.ndarray`, ``(N)``
700 :type lat0: float
701 :type lon0: float
703 :return: Array with latitudes and longitudes in [deg].
704 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
705 '''
707 a = distance_rad
708 gamma = azimuth_rad
710 b = math.pi/2.-lat0*d2r
712 alphasign = 1.
713 alphasign = num.where(gamma < 0, -1., 1.)
714 gamma = num.abs(gamma)
716 c = num.arccos(clip(
717 num.cos(a)*num.cos(b)+num.sin(a)*num.sin(b)*num.cos(gamma), -1., 1.))
719 alpha = num.arcsin(clip(
720 num.sin(a)*num.sin(gamma)/num.sin(c), -1., 1.))
722 alpha = num.where(
723 num.cos(a)-num.cos(b)*num.cos(c) < 0,
724 num.where(alpha > 0, math.pi-alpha, -math.pi-alpha),
725 alpha)
727 lat = r2d * (math.pi/2. - c)
728 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
730 return lat, lon
733def ne_to_latlon_alternative_method(lat0, lon0, north_m, east_m):
734 '''
735 Transform local cartesian coordinates to latitude and longitude.
737 Like :py:func:`pyrocko.orthodrome.ne_to_latlon`,
738 but this method (implementation below), although it should be numerically
739 more stable, suffers problems at points which are *across the pole*
740 as seen from the cartesian origin.
742 .. math::
744 \\text{distance change:}\\; \\Delta {{\\bf a}_i} &=
745 \\sqrt{{\\bf{y}}^2_i + {\\bf{x}}^2_i }/ \\mathrm{R_E},\\\\
746 \\text{azimuth change:}\\; \\Delta {\\bf \\gamma}_i &=
747 \\arctan( {\\bf x}_i {\\bf y}_i). \\\\
748 \\mathrm{b} &=
749 \\frac{\\pi}{2} -\\frac{\\pi}{180} \\;\\mathrm{lat_0}\\\\
751 .. math::
753 {{\\bf z}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i -
754 \\mathrm{b}}{2} \\right)}
755 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
756 {{\\bf n}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i +
757 \\mathrm{b}}{2} \\right)}
758 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
759 {{\\bf z}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i -
760 \\mathrm{b}}{2} \\right)}
761 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
762 {{\\bf n}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i +
763 \\mathrm{b}}{2} \\right)}
764 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
765 {{\\bf t}_1}_i &= \\arctan{\\left( \\frac{{{\\bf z}_1}_i}
766 {{{\\bf n}_1}_i} \\right) }\\\\
767 {{\\bf t}_2}_i &= \\arctan{\\left( \\frac{{{\\bf z}_2}_i}
768 {{{\\bf n}_2}_i} \\right) } \\\\[0.5cm]
769 c &= \\begin{cases}
770 2 \\cdot \\arccos \\left( {{\\bf z}_1}_i / \\sin({{\\bf t}_1}_i)
771 \\right),\\; \\text{if }
772 |\\sin({{\\bf t}_1}_i)| >
773 |\\sin({{\\bf t}_2}_i)|,\\; \\text{else} \\\\
774 2 \\cdot \\arcsin{\\left( {{\\bf z}_2}_i /
775 \\sin({{\\bf t}_2}_i) \\right)}.
776 \\end{cases}\\\\
778 .. math::
780 {\\bf {lat}}_i &= \\frac{180}{ \\pi } \\left( \\frac{\\pi}{2}
781 - {\\bf {c}}_i \\right) \\\\
782 {\\bf {lon}}_i &= {\\bf {lon}}_0 + \\frac{180}{ \\pi }
783 \\frac{\\gamma_i}{|\\gamma_i|},
784 \\text{ with}\\; \\gamma \\in [-\\pi,\\pi]
786 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
787 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
788 :param north_m: Northing distances from origin in [m].
789 :param east_m: Easting distances from origin in [m].
790 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
791 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
792 :type lat0: float
793 :type lon0: float
795 :return: Array with latitudes and longitudes in [deg].
796 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
797 '''
799 b = math.pi/2.-lat0*d2r
800 a = num.sqrt(north_m**2+east_m**2)/earthradius
802 gamma = num.arctan2(east_m, north_m)
803 alphasign = 1.
804 alphasign = num.where(gamma < 0., -1., 1.)
805 gamma = num.abs(gamma)
807 z1 = num.cos((a-b)/2.)*num.cos(gamma/2.)
808 n1 = num.cos((a+b)/2.)*num.sin(gamma/2.)
809 z2 = num.sin((a-b)/2.)*num.cos(gamma/2.)
810 n2 = num.sin((a+b)/2.)*num.sin(gamma/2.)
811 t1 = num.arctan2(z1, n1)
812 t2 = num.arctan2(z2, n2)
814 alpha = t1 + t2
816 sin_t1 = num.sin(t1)
817 sin_t2 = num.sin(t2)
818 c = num.where(
819 num.abs(sin_t1) > num.abs(sin_t2),
820 num.arccos(z1/sin_t1)*2.,
821 num.arcsin(z2/sin_t2)*2.)
823 lat = r2d * (math.pi/2. - c)
824 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
825 return lat, lon
828def latlon_to_ne(*args):
829 '''
830 Relative cartesian coordinates with respect to a reference location.
832 For two locations, a reference location A and another location B, given in
833 geographical coordinates in degrees, the corresponding cartesian
834 coordinates are calculated.
835 Assisting functions are :py:func:`pyrocko.orthodrome.azimuth` and
836 :py:func:`pyrocko.orthodrome.distance_accurate50m`.
838 .. math::
840 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
841 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
842 \\mathrm{)}\\\\[0.3cm]
844 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180}
845 \\varphi_{\\mathrm{azi},AB} )\\\\
846 e &= D_{AB} \\cdot
847 \\sin( \\frac{\\pi }{180} \\varphi_{\\mathrm{azi},AB})
849 :param refloc: Location reference point.
850 :type refloc: :py:class:`pyrocko.orthodrome.Loc`
851 :param loc: Location of interest.
852 :type loc: :py:class:`pyrocko.orthodrome.Loc`
854 :return: Northing and easting from refloc to location in [m].
855 :rtype: tuple[float, float]
857 '''
859 azi = azimuth(*args)
860 dist = distance_accurate50m(*args)
861 n, e = math.cos(azi*d2r)*dist, math.sin(azi*d2r)*dist
862 return n, e
865def latlon_to_ne_numpy(lat0, lon0, lat, lon):
866 '''
867 Relative cartesian coordinates with respect to a reference location.
869 For two locations, a reference location (``lat0``, ``lon0``) and another
870 location B, given in geographical coordinates in degrees,
871 the corresponding cartesian coordinates are calculated.
872 Assisting functions are :py:func:`azimuth`
873 and :py:func:`distance_accurate50m`.
875 :param lat0: Latitude of the reference location in [deg].
876 :param lon0: Longitude of the reference location in [deg].
877 :param lat: Latitude of the absolute location in [deg].
878 :param lon: Longitude of the absolute location in [deg].
880 :return: ``(n, e)``: relative north and east positions in [m].
881 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
883 Implemented formulations:
885 .. math::
887 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
888 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
889 \\mathrm{)}\\\\[0.3cm]
891 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180} \\varphi_{
892 \\mathrm{azi},AB} )\\\\
893 e &= D_{AB} \\cdot \\sin( \\frac{\\pi }{180} \\varphi_{
894 \\mathrm{azi},AB} )
895 '''
897 azi = azimuth_numpy(lat0, lon0, lat, lon)
898 dist = distance_accurate50m_numpy(lat0, lon0, lat, lon)
899 n = num.cos(azi*d2r)*dist
900 e = num.sin(azi*d2r)*dist
901 return n, e
904_wgs84 = None
907def get_wgs84():
908 global _wgs84
909 if _wgs84 is None:
910 from geographiclib.geodesic import Geodesic
911 _wgs84 = Geodesic.WGS84
913 return _wgs84
916def amap(n):
917 def wrap(f):
918 if n == 1:
919 def func(*args):
920 it = num.nditer(args + (None,))
921 for ops in it:
922 ops[-1][...] = f(*ops[:-1])
924 return it.operands[-1]
925 elif n == 2:
926 def func(*args):
927 it = num.nditer(args + (None, None))
928 for ops in it:
929 ops[-2][...], ops[-1][...] = f(*ops[:-2])
931 return it.operands[-2], it.operands[-1]
933 return func
934 return wrap
937@amap(2)
938def ne_to_latlon2(lat0, lon0, north_m, east_m):
939 wgs84 = get_wgs84()
940 az = num.arctan2(east_m, north_m)*r2d
941 dist = num.sqrt(east_m**2 + north_m**2)
942 x = wgs84.Direct(lat0, lon0, az, dist)
943 return x['lat2'], x['lon2']
946@amap(2)
947def latlon_to_ne2(lat0, lon0, lat1, lon1):
948 wgs84 = get_wgs84()
949 x = wgs84.Inverse(lat0, lon0, lat1, lon1)
950 dist = x['s12']
951 az = x['azi1']
952 n = num.cos(az*d2r)*dist
953 e = num.sin(az*d2r)*dist
954 return n, e
957@amap(1)
958def distance_accurate15nm(lat1, lon1, lat2, lon2):
959 wgs84 = get_wgs84()
960 return wgs84.Inverse(lat1, lon1, lat2, lon2)['s12']
963def positive_region(region):
964 '''
965 Normalize parameterization of a rectangular geographical region.
967 :param region: ``(west, east, south, north)`` in [deg].
968 :type region: :py:class:`tuple` of :py:class:`float`
970 :returns: ``(west, east, south, north)``, where ``west <= east`` and
971 where ``west`` and ``east`` are in the range ``[-180., 180.+360.]``.
972 :rtype: :py:class:`tuple` of :py:class:`float`
973 '''
974 west, east, south, north = [float(x) for x in region]
976 assert -180. - 360. <= west < 180.
977 assert -180. < east <= 180. + 360.
978 assert -90. <= south < 90.
979 assert -90. < north <= 90.
981 if east < west:
982 east += 360.
984 if west < -180.:
985 west += 360.
986 east += 360.
988 return (west, east, south, north)
991def points_in_region(p, region):
992 '''
993 Check what points are contained in a rectangular geographical region.
995 :param p: ``(lat, lon)`` pairs in [deg].
996 :type p: :py:class:`numpy.ndarray` ``(N, 2)``
997 :param region: ``(west, east, south, north)`` region boundaries in [deg].
998 :type region: :py:class:`tuple` of :py:class:`float`
1000 :returns: Mask, returning ``True`` for each point within the region.
1001 :rtype: :py:class:`numpy.ndarray` of bool, shape ``(N)``
1002 '''
1004 w, e, s, n = positive_region(region)
1005 return num.logical_and(
1006 num.logical_and(s <= p[:, 0], p[:, 0] <= n),
1007 num.logical_or(
1008 num.logical_and(w <= p[:, 1], p[:, 1] <= e),
1009 num.logical_and(w-360. <= p[:, 1], p[:, 1] <= e-360.)))
1012def point_in_region(p, region):
1013 '''
1014 Check if a point is contained in a rectangular geographical region.
1016 :param p: ``(lat, lon)`` in [deg].
1017 :type p: :py:class:`tuple` of :py:class:`float`
1018 :param region: ``(west, east, south, north)`` region boundaries in [deg].
1019 :type region: :py:class:`tuple` of :py:class:`float`
1021 :returns: ``True``, if point is in region, else ``False``.
1022 :rtype: bool
1023 '''
1025 w, e, s, n = positive_region(region)
1026 return num.logical_and(
1027 num.logical_and(s <= p[0], p[0] <= n),
1028 num.logical_or(
1029 num.logical_and(w <= p[1], p[1] <= e),
1030 num.logical_and(w-360. <= p[1], p[1] <= e-360.)))
1033def radius_to_region(lat, lon, radius):
1034 '''
1035 Get a rectangular region which fully contains a given circular region.
1037 :param lat: Latitude of the center point of circular region in [deg].
1038 :type lat: float
1039 :param lon: Longitude of the center point of circular region in [deg].
1040 :type lon: float
1041 :param radius: Radius of circular region in [m].
1042 :type radius: float
1044 :returns: Rectangular region as ``(east, west, south, north)`` in [deg] or
1045 ``None``.
1046 :rtype: tuple[float, float, float, float]
1047 '''
1048 radius_deg = radius * m2d
1049 if radius_deg < 45.:
1050 lat_min = max(-90., lat - radius_deg)
1051 lat_max = min(90., lat + radius_deg)
1052 absmaxlat = max(abs(lat_min), abs(lat_max))
1053 if absmaxlat > 89:
1054 lon_min = -180.
1055 lon_max = 180.
1056 else:
1057 lon_min = max(
1058 -180. - 360.,
1059 lon - radius_deg / math.cos(absmaxlat*d2r))
1060 lon_max = min(
1061 180. + 360.,
1062 lon + radius_deg / math.cos(absmaxlat*d2r))
1064 lon_min, lon_max, lat_min, lat_max = positive_region(
1065 (lon_min, lon_max, lat_min, lat_max))
1067 return lon_min, lon_max, lat_min, lat_max
1069 else:
1070 return None
1073def geographic_midpoint(
1074 lats, lons, weights=None, depths=None, earthradius=earthradius):
1076 '''
1077 Calculate geographic midpoints by finding the center of gravity.
1079 This method suffers from instabilities if points are centered around the
1080 poles.
1082 :param lats: Latitudes in [deg].
1083 :param lons: Longitudes in [deg].
1084 :param weights: Weighting factors.
1085 :param depths: Depths in [m].
1086 :type lats: :py:class:`numpy.ndarray`, ``(N)``
1087 :type lons: :py:class:`numpy.ndarray`, ``(N)``
1088 :type weights: optional, :py:class:`numpy.ndarray`, ``(N)``
1089 :type depths: optional, :py:class:`numpy.ndarray`, ``(N)``
1091 :return: Latitudes and longitudes of the midpoint in [deg] (and depth [m]
1092 if depths are given).
1093 :rtype: ``(lat, lon)`` or ``(lat, lon, depth)``
1094 '''
1095 if not weights:
1096 weights = num.ones(len(lats))
1098 total_weigth = num.sum(weights)
1099 weights /= total_weigth
1100 lats = lats * d2r
1101 lons = lons * d2r
1102 if depths is not None:
1103 radii = (earthradius - depths) / earthradius
1104 else:
1105 radii = 1.0
1107 x = num.sum(num.cos(lats) * num.cos(lons) * weights * radii)
1108 y = num.sum(num.cos(lats) * num.sin(lons) * weights * radii)
1109 z = num.sum(num.sin(lats) * weights * radii)
1111 lon = num.arctan2(y, x)
1112 hyp = num.sqrt(x**2 + y**2)
1113 lat = num.arctan2(z, hyp)
1114 depth = earthradius - num.sqrt(x**2 + y**2 + z**2) * earthradius
1116 if depths is None:
1117 return lat/d2r, lon/d2r
1118 else:
1119 return lat/d2r, lon/d2r, depth
1122def geographic_midpoint_locations(
1123 locations, weights=None, include_depth=False):
1125 if not include_depth:
1126 coords = num.array([loc.effective_latlon
1127 for loc in locations])
1128 return geographic_midpoint(
1129 coords[:, 0], coords[:, 1], weights)
1130 else:
1131 coords = num.array([loc.effective_latlon + (loc.depth,)
1132 for loc in locations])
1133 return geographic_midpoint(
1134 coords[:, 0], coords[:, 1], weights=weights, depths=coords[:, 2])
1137def geodetic_to_ecef(lat, lon, alt):
1138 '''
1139 Convert geodetic coordinates to Earth-Centered, Earth-Fixed (ECEF)
1140 Cartesian coordinates. [#1]_ [#2]_
1142 :param lat: Geodetic latitude in [deg].
1143 :param lon: Geodetic longitude in [deg].
1144 :param alt: Geodetic altitude (height) in [m] (positive for points outside
1145 the geoid).
1146 :type lat: float
1147 :type lon: float
1148 :type alt: float
1150 :return: ECEF Cartesian coordinates (X, Y, Z) in [m].
1151 :rtype: tuple[float, float, float]
1153 .. [#1] https://en.wikipedia.org/wiki/ECEF
1154 .. [#2] https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1155 #From_geodetic_to_ECEF_coordinates
1156 '''
1158 f = earth_oblateness
1159 a = earthradius_equator
1160 e2 = 2*f - f**2
1162 lat, lon = num.radians(lat), num.radians(lon)
1163 # Normal (plumb line)
1164 N = a / num.sqrt(1.0 - (e2 * num.sin(lat)**2))
1166 X = (N+alt) * num.cos(lat) * num.cos(lon)
1167 Y = (N+alt) * num.cos(lat) * num.sin(lon)
1168 Z = (N*(1.0-e2) + alt) * num.sin(lat)
1170 return (X, Y, Z)
1173def ecef_to_geodetic(X, Y, Z):
1174 '''
1175 Convert Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates to
1176 geodetic coordinates (Ferrari's solution).
1178 :param X: ECEF X coordinate [m].
1179 :type X: float
1180 :param Y: ECEF Y coordinate [m].
1181 :type Y: float
1182 :param Z: ECEF Z coordinate [m].
1183 :type Z: float
1185 :return: Geodetic coordinates (lat, lon, alt). Latitude and longitude are
1186 in [deg] and altitude is in [m]
1187 (positive for points outside the geoid).
1188 :rtype: :py:class:`tuple` of :py:class:`float`
1190 .. seealso ::
1191 https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1192 #The_application_of_Ferrari.27s_solution
1193 '''
1194 f = earth_oblateness
1195 a = earthradius_equator
1196 b = a * (1. - f)
1197 e2 = 2.*f - f**2
1199 # usefull
1200 a2 = a**2
1201 b2 = b**2
1202 e4 = e2**2
1203 X2 = X**2
1204 Y2 = Y**2
1205 Z2 = Z**2
1207 r = num.sqrt(X2 + Y2)
1208 r2 = r**2
1210 e_prime2 = (a2 - b2)/b2
1211 E2 = a2 - b2
1212 F = 54. * b2 * Z2
1213 G = r2 + (1.-e2)*Z2 - (e2*E2)
1214 C = (e4 * F * r2) / (G**3)
1215 S = cbrt(1. + C + num.sqrt(C**2 + 2.*C))
1216 P = F / (3. * (S + 1./S + 1.)**2 * G**2)
1217 Q = num.sqrt(1. + (2.*e4*P))
1219 dum1 = -(P*e2*r) / (1.+Q)
1220 dum2 = 0.5 * a2 * (1. + 1./Q)
1221 dum3 = (P * (1.-e2) * Z2) / (Q * (1.+Q))
1222 dum4 = 0.5 * P * r2
1223 r0 = dum1 + num.sqrt(dum2 - dum3 - dum4)
1225 U = num.sqrt((r - e2*r0)**2 + Z2)
1226 V = num.sqrt((r - e2*r0)**2 + (1.-e2)*Z2)
1227 Z0 = (b2*Z) / (a*V)
1229 alt = U * (1. - (b2 / (a*V)))
1230 lat = num.arctan((Z + e_prime2 * Z0)/r)
1231 lon = num.arctan2(Y, X)
1233 return (lat*r2d, lon*r2d, alt)
1236class Farside(Exception):
1237 pass
1240def latlon_to_xyz(latlons):
1241 if latlons.ndim == 1:
1242 return latlon_to_xyz(latlons[num.newaxis, :])[0]
1244 points = num.zeros((latlons.shape[0], 3))
1245 lats = latlons[:, 0]
1246 lons = latlons[:, 1]
1247 points[:, 0] = num.cos(lats*d2r) * num.cos(lons*d2r)
1248 points[:, 1] = num.cos(lats*d2r) * num.sin(lons*d2r)
1249 points[:, 2] = num.sin(lats*d2r)
1250 return points
1253def xyz_to_latlon(xyz):
1254 if xyz.ndim == 1:
1255 return xyz_to_latlon(xyz[num.newaxis, :])[0]
1257 latlons = num.zeros((xyz.shape[0], 2))
1258 latlons[:, 0] = num.arctan2(
1259 xyz[:, 2], num.sqrt(xyz[:, 0]**2 + xyz[:, 1]**2)) * r2d
1260 latlons[:, 1] = num.arctan2(
1261 xyz[:, 1], xyz[:, 0]) * r2d
1262 return latlons
1265def rot_to_00(lat, lon):
1266 rot0 = euler_to_matrix(0., -90.*d2r, 0.)
1267 rot1 = euler_to_matrix(-d2r*lat, 0., -d2r*lon)
1268 return num.dot(rot0.T, num.dot(rot1, rot0)).T
1271def distances3d(a, b):
1272 return num.sqrt(num.sum((a-b)**2, axis=a.ndim-1))
1275def circulation(points2):
1276 return num.sum(
1277 (points2[1:, 0] - points2[:-1, 0])
1278 * (points2[1:, 1] + points2[:-1, 1]))
1281def stereographic(points):
1282 dists = distances3d(points[1:, :], points[:-1, :])
1283 if dists.size > 0:
1284 maxdist = num.max(dists)
1285 cutoff = maxdist**2 / 2.
1286 else:
1287 cutoff = 1.0e-5
1289 points = points.copy()
1290 if num.any(points[:, 0] < -1. + cutoff):
1291 raise Farside()
1293 points_out = points[:, 1:].copy()
1294 factor = 1.0 / (1.0 + points[:, 0])
1295 points_out *= factor[:, num.newaxis]
1297 return points_out
1300def stereographic_poly(points):
1301 dists = distances3d(points[1:, :], points[:-1, :])
1302 if dists.size > 0:
1303 maxdist = num.max(dists)
1304 cutoff = maxdist**2 / 2.
1305 else:
1306 cutoff = 1.0e-5
1308 points = points.copy()
1309 if num.any(points[:, 0] < -1. + cutoff):
1310 raise Farside()
1312 points_out = points[:, 1:].copy()
1313 factor = 1.0 / (1.0 + points[:, 0])
1314 points_out *= factor[:, num.newaxis]
1316 if circulation(points_out) >= 0:
1317 raise Farside()
1319 return points_out
1322def gnomonic_x(points, cutoff=0.01):
1323 points_out = points[:, 1:].copy()
1324 if num.any(points[:, 0] < cutoff):
1325 raise Farside()
1327 factor = 1.0 / points[:, 0]
1328 points_out *= factor[:, num.newaxis]
1329 return points_out
1332def cneg(i, x):
1333 if i == 1:
1334 return x
1335 else:
1336 return num.logical_not(x)
1339def contains_points(polygon, points):
1340 '''
1341 Test which points are inside polygon on a sphere.
1343 The inside of the polygon is defined as the area which is to the left hand
1344 side of an observer walking the polygon line, points in order, on the
1345 sphere. Lines between the polygon points are treated as great circle paths.
1346 The polygon may be arbitrarily complex, as long as it does not have any
1347 crossings or thin parts with zero width. The polygon may contain the poles
1348 and is allowed to wrap around the sphere multiple times.
1350 The algorithm works by consecutive cutting of the polygon into (almost)
1351 hemispheres and subsequent Gnomonic projections to perform the
1352 point-in-polygon tests on a 2D plane.
1354 :param polygon: Point coordinates defining the polygon [deg].
1355 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1356 0=lat, 1=lon
1357 :param points: Coordinates of points to test [deg].
1358 :type points: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1359 0=lat, 1=lon
1361 :returns: Boolean mask array.
1362 :rtype: :py:class:`numpy.ndarray` of shape ``(N,)``.
1363 '''
1365 and_ = num.logical_and
1367 points_xyz = latlon_to_xyz(points)
1368 mask_x = 0. <= points_xyz[:, 0]
1369 mask_y = 0. <= points_xyz[:, 1]
1370 mask_z = 0. <= points_xyz[:, 2]
1372 result = num.zeros(points.shape[0], dtype=int)
1374 for ix in [-1, 1]:
1375 for iy in [-1, 1]:
1376 for iz in [-1, 1]:
1377 mask = and_(
1378 and_(cneg(ix, mask_x), cneg(iy, mask_y)),
1379 cneg(iz, mask_z))
1381 center_xyz = num.array([ix, iy, iz], dtype=float)
1383 lat, lon = xyz_to_latlon(center_xyz)
1384 rot = rot_to_00(lat, lon)
1386 points_rot_xyz = num.dot(rot, points_xyz[mask, :].T).T
1387 points_rot_pro = gnomonic_x(points_rot_xyz)
1389 offset = 0.01
1391 poly_xyz = latlon_to_xyz(polygon)
1392 poly_rot_xyz = num.dot(rot, poly_xyz.T).T
1393 poly_rot_xyz[:, 0] -= offset
1394 groups = spoly_cut([poly_rot_xyz], axis=0)
1396 for poly_rot_group_xyz in groups[1]:
1397 poly_rot_group_xyz[:, 0] += offset
1399 poly_rot_group_pro = gnomonic_x(
1400 poly_rot_group_xyz)
1402 if circulation(poly_rot_group_pro) > 0:
1403 result[mask] += path_contains_points(
1404 poly_rot_group_pro, points_rot_pro)
1405 else:
1406 result[mask] -= path_contains_points(
1407 poly_rot_group_pro, points_rot_pro)
1409 return result.astype(bool)
1412def contains_point(polygon, point):
1413 '''
1414 Test if point is inside polygon on a sphere.
1416 Convenience wrapper to :py:func:`contains_points` to test a single point.
1418 :param polygon: Point coordinates defining the polygon [deg].
1419 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1420 0=lat, 1=lon
1421 :param point: Coordinates ``(lat, lon)`` of point to test [deg].
1422 :type point: :py:class:`tuple` of :py:class:`float`
1424 :returns: ``True``, if point is located within polygon, else ``False``.
1425 :rtype: bool
1426 '''
1428 return bool(
1429 contains_points(polygon, num.asarray(point)[num.newaxis, :])[0])