1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import division, absolute_import
7from functools import wraps, lru_cache
8import math
9import numpy as num
11from .moment_tensor import euler_to_matrix
12from .config import config
13from .plot.beachball import spoly_cut
15from matplotlib.path import Path
17d2r = math.pi/180.
18r2d = 1./d2r
19earth_oblateness = 1./298.257223563
20earthradius_equator = 6378.14 * 1000.
21earthradius = config().earthradius
22d2m = earthradius_equator*math.pi/180.
23m2d = 1./d2m
25_testpath = Path([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)], closed=True)
27raise_if_slow_path_contains_points = False
30class Slow(Exception):
31 pass
34if hasattr(_testpath, 'contains_points') and num.all(
35 _testpath.contains_points([(0.5, 0.5), (1.5, 0.5)]) == [True, False]):
37 def path_contains_points(verts, points):
38 p = Path(verts, closed=True)
39 return p.contains_points(points).astype(num.bool)
41else:
42 # work around missing contains_points and bug in matplotlib ~ v1.2.0
44 def path_contains_points(verts, points):
45 if raise_if_slow_path_contains_points:
46 # used by unit test to skip slow gshhg_example.py
47 raise Slow()
49 p = Path(verts, closed=True)
50 result = num.zeros(points.shape[0], dtype=num.bool)
51 for i in range(result.size):
52 result[i] = p.contains_point(points[i, :])
54 return result
57try:
58 cbrt = num.cbrt
59except AttributeError:
60 def cbrt(x):
61 return x**(1./3.)
64def float_array_broadcast(*args):
65 return num.broadcast_arrays(*[
66 num.asarray(x, dtype=float) for x in args])
69class Loc(object):
70 '''
71 Simple location representation.
73 :attrib lat: Latitude in [deg].
74 :attrib lon: Longitude in [deg].
75 '''
76 __slots__ = ['lat', 'lon']
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
134@lru_cache
135def cosdelta(*args):
136 '''
137 Cosine of the angular distance between two points ``a`` and ``b`` on a
138 sphere.
140 This function (find implementation below) returns the cosine of the
141 distance angle 'delta' between two points ``a`` and ``b``, coordinates of
142 which are expected to be given in geographical coordinates and in degrees.
143 For numerical stability a maximum of 1.0 is enforced.
145 .. math::
147 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
148 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
149 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
150 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\[0.5cm]
152 \\cos(\\Delta) = \\min( 1.0, \\quad \\sin( A_{\\mathrm{lat'}})
153 \\sin( B_{\\mathrm{lat'}} ) +
154 \\cos(A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
155 \\cos( B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )
157 :param a: Location point A.
158 :type a: :py:class:`pyrocko.orthodrome.Loc`
159 :param b: Location point B.
160 :type b: :py:class:`pyrocko.orthodrome.Loc`
162 :return: Cosdelta.
163 :rtype: float
164 '''
166 alat, alon, blat, blon = _latlon_pair(args)
168 return min(
169 1.0,
170 math.sin(alat*d2r) * math.sin(blat*d2r) +
171 math.cos(alat*d2r) * math.cos(blat*d2r) *
172 math.cos(d2r*(blon-alon)))
175def cosdelta_numpy(a_lats, a_lons, b_lats, b_lons):
176 '''
177 Cosine of the angular distance between two points ``a`` and ``b`` on a
178 sphere.
180 This function returns the cosines of the distance
181 angles *delta* between two points ``a`` and ``b`` given as
182 :py:class:`numpy.ndarray`.
183 The coordinates are expected to be given in geographical coordinates
184 and in degrees. For numerical stability a maximum of ``1.0`` is enforced.
186 Please find the details of the implementation in the documentation of
187 the function :py:func:`pyrocko.orthodrome.cosdelta` above.
189 :param a_lats: Latitudes in [deg] point A.
190 :param a_lons: Longitudes in [deg] point A.
191 :param b_lats: Latitudes in [deg] point B.
192 :param b_lons: Longitudes in [deg] point B.
193 :type a_lats: :py:class:`numpy.ndarray`
194 :type a_lons: :py:class:`numpy.ndarray`
195 :type b_lats: :py:class:`numpy.ndarray`
196 :type b_lons: :py:class:`numpy.ndarray`
198 :return: Cosdelta.
199 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
200 '''
201 return num.minimum(
202 1.0,
203 num.sin(a_lats*d2r) * num.sin(b_lats*d2r) +
204 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
205 num.cos(d2r*(b_lons-a_lons)))
208@lru_cache
209def azimuth(*args):
210 '''
211 Azimuth calculation.
213 This function (find implementation below) returns azimuth ...
214 between points ``a`` and ``b``, coordinates of
215 which are expected to be given in geographical coordinates and in degrees.
217 .. math::
219 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
220 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
221 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
222 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\
224 \\varphi_{\\mathrm{azi},AB} = \\frac{180}{\\pi} \\arctan \\left[
225 \\frac{
226 \\cos( A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
227 \\sin(B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )}
228 {\\sin ( B_{\\mathrm{lat'}} ) - \\sin( A_{\\mathrm{lat'}}
229 cosdelta) } \\right]
231 :param a: Location point A.
232 :type a: :py:class:`pyrocko.orthodrome.Loc`
233 :param b: Location point B.
234 :type b: :py:class:`pyrocko.orthodrome.Loc`
236 :return: Azimuth in degree
237 '''
239 alat, alon, blat, blon = _latlon_pair(args)
241 return r2d*math.atan2(
242 math.cos(alat*d2r) * math.cos(blat*d2r) *
243 math.sin(d2r*(blon-alon)),
244 math.sin(d2r*blat) - math.sin(d2r*alat) * cosdelta(
245 alat, alon, blat, blon))
248def azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta=None):
249 '''
250 Calculation of the azimuth (*track angle*) from a location A towards B.
252 This function returns azimuths (*track angles*) from locations A towards B
253 given in :py:class:`numpy.ndarray`. Coordinates are expected to be given in
254 geographical coordinates and in degrees.
256 Please find the details of the implementation in the documentation of the
257 function :py:func:`pyrocko.orthodrome.azimuth`.
260 :param a_lats: Latitudes in [deg] point A.
261 :param a_lons: Longitudes in [deg] point A.
262 :param b_lats: Latitudes in [deg] point B.
263 :param b_lons: Longitudes in [deg] point B.
264 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
265 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
266 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
267 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
269 :return: Azimuths in [deg].
270 :rtype: :py:class:`numpy.ndarray`, ``(N)``
271 '''
272 if _cosdelta is None:
273 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
275 return r2d*num.arctan2(
276 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
277 num.sin(d2r*(b_lons-a_lons)),
278 num.sin(d2r*b_lats) - num.sin(d2r*a_lats) * _cosdelta)
281@lru_cache
282def azibazi(*args, **kwargs):
283 '''
284 Azimuth and backazimuth from location A towards B and back.
286 :returns: Azimuth in [deg] from A to B, back azimuth in [deg] from B to A.
287 :rtype: tuple[float, float]
288 '''
290 alat, alon, blat, blon = _latlon_pair(args)
291 if alat == blat and alon == blon:
292 return 0., 180.
294 implementation = kwargs.get('implementation', 'c')
295 assert implementation in ('c', 'python')
296 if implementation == 'c':
297 from pyrocko import orthodrome_ext
298 return orthodrome_ext.azibazi(alat, alon, blat, blon)
300 cd = cosdelta(alat, alon, blat, blon)
301 azi = r2d*math.atan2(
302 math.cos(alat*d2r) * math.cos(blat*d2r) *
303 math.sin(d2r*(blon-alon)),
304 math.sin(d2r*blat) - math.sin(d2r*alat) * cd)
305 bazi = r2d*math.atan2(
306 math.cos(blat*d2r) * math.cos(alat*d2r) *
307 math.sin(d2r*(alon-blon)),
308 math.sin(d2r*alat) - math.sin(d2r*blat) * cd)
310 return azi, bazi
313def azibazi_numpy(a_lats, a_lons, b_lats, b_lons, implementation='c'):
314 '''
315 Azimuth and backazimuth from location A towards B and back.
317 Arguments are given as :py:class:`numpy.ndarray`.
319 :param a_lats: Latitude(s) in [deg] of point A.
320 :type a_lats: :py:class:`numpy.ndarray`
321 :param a_lons: Longitude(s) in [deg] of point A.
322 :type a_lons: :py:class:`numpy.ndarray`
323 :param b_lats: Latitude(s) in [deg] of point B.
324 :type b_lats: :py:class:`numpy.ndarray`
325 :param b_lons: Longitude(s) in [deg] of point B.
326 :type b_lons: :py:class:`numpy.ndarray`
328 :returns: Azimuth(s) in [deg] from A to B,
329 back azimuth(s) in [deg] from B to A.
330 :rtype: :py:class:`numpy.ndarray`, :py:class:`numpy.ndarray`
331 '''
333 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
334 a_lats, a_lons, b_lats, b_lons)
336 assert implementation in ('c', 'python')
337 if implementation == 'c':
338 from pyrocko import orthodrome_ext
339 return orthodrome_ext.azibazi_numpy(a_lats, a_lons, b_lats, b_lons)
341 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
342 azis = azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta)
343 bazis = azimuth_numpy(b_lats, b_lons, a_lats, a_lons, _cosdelta)
345 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
346 ii_eq = num.where(eq)[0]
347 azis[ii_eq] = 0.0
348 bazis[ii_eq] = 180.0
349 return azis, bazis
352def azidist_numpy(*args):
353 '''
354 Calculation of the azimuth (*track angle*) and the distance from locations
355 A towards B on a sphere.
357 The assisting functions used are :py:func:`pyrocko.orthodrome.cosdelta` and
358 :py:func:`pyrocko.orthodrome.azimuth`
360 :param a_lats: Latitudes in [deg] point A.
361 :param a_lons: Longitudes in [deg] point A.
362 :param b_lats: Latitudes in [deg] point B.
363 :param b_lons: Longitudes in [deg] point B.
364 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
365 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
366 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
367 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
369 :return: Azimuths in [deg], distances in [deg].
370 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
371 '''
372 _cosdelta = cosdelta_numpy(*args)
373 _azimuths = azimuth_numpy(_cosdelta=_cosdelta, *args)
374 return _azimuths, r2d*num.arccos(_cosdelta)
377def distance_accurate50m(*args, **kwargs):
378 '''
379 Accurate distance calculation based on a spheroid of rotation.
381 Function returns distance in meter between points A and B, coordinates of
382 which must be given in geographical coordinates and in degrees.
383 The returned distance should be accurate to 50 m using WGS84.
384 Values for the Earth's equator radius and the Earth's oblateness
385 (``f_oblate``) are defined in the pyrocko configuration file
386 :py:class:`pyrocko.config`.
388 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
390 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
391 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
393 .. math::
395 F = \\frac{\\pi}{180}
396 \\frac{(A_{lat} + B_{lat})}{2}, \\quad
397 G = \\frac{\\pi}{180}
398 \\frac{(A_{lat} - B_{lat})}{2}, \\quad
399 l = \\frac{\\pi}{180}
400 \\frac{(A_{lon} - B_{lon})}{2} \\quad
401 \\\\[0.5cm]
402 S = \\sin^2(G) \\cdot \\cos^2(l) +
403 \\cos^2(F) \\cdot \\sin^2(l), \\quad \\quad
404 C = \\cos^2(G) \\cdot \\cos^2(l) +
405 \\sin^2(F) \\cdot \\sin^2(l)
407 .. math::
409 w = \\arctan \\left( \\sqrt{ \\frac{S}{C}} \\right) , \\quad
410 r = \\sqrt{\\frac{S}{C} }
412 The spherical-earth distance D between A and B, can be given with:
414 .. math::
416 D_{sphere} = 2w \\cdot R_{equator}
418 The oblateness of the Earth requires some correction with
419 correction factors h1 and h2:
421 .. math::
423 h_1 = \\frac{3r - 1}{2C}, \\quad
424 h_2 = \\frac{3r +1 }{2S}\\\\[0.5cm]
426 D = D_{\\mathrm{sphere}} \\cdot [ 1 + h_1 \\,f_{\\mathrm{oblate}}
427 \\cdot \\sin^2(F)
428 \\cos^2(G) - h_2\\, f_{\\mathrm{oblate}}
429 \\cdot \\cos^2(F) \\sin^2(G)]
432 :param a: Location point A.
433 :type a: :py:class:`pyrocko.orthodrome.Loc`
434 :param b: Location point B.
435 :type b: :py:class:`pyrocko.orthodrome.Loc`
437 :return: Distance in [m].
438 :rtype: float
439 '''
441 alat, alon, blat, blon = _latlon_pair(args)
443 implementation = kwargs.get('implementation', 'c')
444 assert implementation in ('c', 'python')
445 if implementation == 'c':
446 from pyrocko import orthodrome_ext
447 return orthodrome_ext.distance_accurate50m(alat, alon, blat, blon)
449 f = (alat + blat)*d2r / 2.
450 g = (alat - blat)*d2r / 2.
451 h = (alon - blon)*d2r / 2.
453 s = math.sin(g)**2 * math.cos(h)**2 + math.cos(f)**2 * math.sin(h)**2
454 c = math.cos(g)**2 * math.cos(h)**2 + math.sin(f)**2 * math.sin(h)**2
456 w = math.atan(math.sqrt(s/c))
458 if w == 0.0:
459 return 0.0
461 r = math.sqrt(s*c)/w
462 d = 2.*w*earthradius_equator
463 h1 = (3.*r-1.)/(2.*c)
464 h2 = (3.*r+1.)/(2.*s)
466 return d * (1. +
467 earth_oblateness * h1 * math.sin(f)**2 * math.cos(g)**2 -
468 earth_oblateness * h2 * math.cos(f)**2 * math.sin(g)**2)
471def distance_accurate50m_numpy(
472 a_lats, a_lons, b_lats, b_lons, implementation='c'):
474 '''
475 Accurate distance calculation based on a spheroid of rotation.
477 Function returns distance in meter between points ``a`` and ``b``,
478 coordinates of which must be given in geographical coordinates and in
479 degrees.
480 The returned distance should be accurate to 50 m using WGS84.
481 Values for the Earth's equator radius and the Earth's oblateness
482 (``f_oblate``) are defined in the pyrocko configuration file
483 :py:class:`pyrocko.config`.
485 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
487 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
488 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
490 .. math::
492 F_i = \\frac{\\pi}{180}
493 \\frac{(a_{lat,i} + a_{lat,i})}{2}, \\quad
494 G_i = \\frac{\\pi}{180}
495 \\frac{(a_{lat,i} - b_{lat,i})}{2}, \\quad
496 l_i= \\frac{\\pi}{180}
497 \\frac{(a_{lon,i} - b_{lon,i})}{2} \\\\[0.5cm]
498 S_i = \\sin^2(G_i) \\cdot \\cos^2(l_i) +
499 \\cos^2(F_i) \\cdot \\sin^2(l_i), \\quad \\quad
500 C_i = \\cos^2(G_i) \\cdot \\cos^2(l_i) +
501 \\sin^2(F_i) \\cdot \\sin^2(l_i)
503 .. math::
505 w_i = \\arctan \\left( \\sqrt{\\frac{S_i}{C_i}} \\right), \\quad
506 r_i = \\sqrt{\\frac{S_i}{C_i} }
508 The spherical-earth distance ``D`` between ``a`` and ``b``,
509 can be given with:
511 .. math::
513 D_{\\mathrm{sphere},i} = 2w_i \\cdot R_{\\mathrm{equator}}
515 The oblateness of the Earth requires some correction with
516 correction factors ``h1`` and ``h2``:
518 .. math::
520 h_{1.i} = \\frac{3r - 1}{2C_i}, \\quad
521 h_{2,i} = \\frac{3r +1 }{2S_i}\\\\[0.5cm]
523 D_{AB,i} = D_{\\mathrm{sphere},i} \\cdot [1 + h_{1,i}
524 \\,f_{\\mathrm{oblate}}
525 \\cdot \\sin^2(F_i)
526 \\cos^2(G_i) - h_{2,i}\\, f_{\\mathrm{oblate}}
527 \\cdot \\cos^2(F_i) \\sin^2(G_i)]
529 :param a_lats: Latitudes in [deg] point A.
530 :param a_lons: Longitudes in [deg] point A.
531 :param b_lats: Latitudes in [deg] point B.
532 :param b_lons: Longitudes in [deg] point B.
533 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
534 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
535 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
536 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
538 :return: Distances in [m].
539 :rtype: :py:class:`numpy.ndarray`, ``(N)``
540 '''
542 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
543 a_lats, a_lons, b_lats, b_lons)
545 assert implementation in ('c', 'python')
546 if implementation == 'c':
547 from pyrocko import orthodrome_ext
548 return orthodrome_ext.distance_accurate50m_numpy(
549 a_lats, a_lons, b_lats, b_lons)
551 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
552 ii_neq = num.where(num.logical_not(eq))[0]
554 if num.all(eq):
555 return num.zeros_like(eq, dtype=float)
557 def extr(x):
558 if isinstance(x, num.ndarray) and x.size > 1:
559 return x[ii_neq]
560 else:
561 return x
563 a_lats = extr(a_lats)
564 a_lons = extr(a_lons)
565 b_lats = extr(b_lats)
566 b_lons = extr(b_lons)
568 f = (a_lats + b_lats)*d2r / 2.
569 g = (a_lats - b_lats)*d2r / 2.
570 h = (a_lons - b_lons)*d2r / 2.
572 s = num.sin(g)**2 * num.cos(h)**2 + num.cos(f)**2 * num.sin(h)**2
573 c = num.cos(g)**2 * num.cos(h)**2 + num.sin(f)**2 * num.sin(h)**2
575 w = num.arctan(num.sqrt(s/c))
577 r = num.sqrt(s*c)/w
579 d = 2.*w*earthradius_equator
580 h1 = (3.*r-1.)/(2.*c)
581 h2 = (3.*r+1.)/(2.*s)
583 dists = num.zeros(eq.size, dtype=float)
584 dists[ii_neq] = d * (
585 1. +
586 earth_oblateness * h1 * num.sin(f)**2 * num.cos(g)**2 -
587 earth_oblateness * h2 * num.cos(f)**2 * num.sin(g)**2)
589 return dists
592def ne_to_latlon(lat0, lon0, north_m, east_m):
593 '''
594 Transform local cartesian coordinates to latitude and longitude.
596 From east and north coordinates (``x`` and ``y`` coordinate
597 :py:class:`numpy.ndarray`) relative to a reference differences in
598 longitude and latitude are calculated, which are effectively changes in
599 azimuth and distance, respectively:
601 .. math::
603 \\text{distance change:}\\; \\Delta {\\bf{a}} &= \\sqrt{{\\bf{y}}^2 +
604 {\\bf{x}}^2 }/ \\mathrm{R_E},
606 \\text{azimuth change:}\\; \\Delta \\bf{\\gamma} &= \\arctan( \\bf{x}
607 / \\bf{y}).
609 The projection used preserves the azimuths of the input points.
611 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
612 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
613 :param north_m: Northing distances from origin in [m].
614 :param east_m: Easting distances from origin in [m].
615 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
616 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
617 :type lat0: float
618 :type lon0: float
620 :return: Array with latitudes and longitudes in [deg].
621 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
623 '''
625 a = num.sqrt(north_m**2+east_m**2)/earthradius
626 gamma = num.arctan2(east_m, north_m)
628 return azidist_to_latlon_rad(lat0, lon0, gamma, a)
631def azidist_to_latlon(lat0, lon0, azimuth_deg, distance_deg):
632 '''
633 Absolute latitudes and longitudes are calculated from relative changes.
635 Convenience wrapper to :py:func:`azidist_to_latlon_rad` with azimuth and
636 distance given in degrees.
638 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
639 :type lat0: float
640 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
641 :type lon0: float
642 :param azimuth_deg: Azimuth from origin in [deg].
643 :type azimuth_deg: :py:class:`numpy.ndarray`, ``(N)``
644 :param distance_deg: Distances from origin in [deg].
645 :type distance_deg: :py:class:`numpy.ndarray`, ``(N)``
647 :return: Array with latitudes and longitudes in [deg].
648 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
649 '''
651 return azidist_to_latlon_rad(
652 lat0, lon0, azimuth_deg/180.*num.pi, distance_deg/180.*num.pi)
655def azidist_to_latlon_rad(lat0, lon0, azimuth_rad, distance_rad):
656 '''
657 Absolute latitudes and longitudes are calculated from relative changes.
659 For numerical stability a range between of ``-1.0`` and ``1.0`` is
660 enforced for ``c`` and ``alpha``.
662 .. math::
664 \\Delta {\\bf a}_i \\; \\text{and} \\; \\Delta \\gamma_i \\;
665 \\text{are relative distances and azimuths from lat0 and lon0 for
666 \\textit{i} source points of a finite source.}
668 .. math::
670 \\mathrm{b} &= \\frac{\\pi}{2} -\\frac{\\pi}{180}\\;\\mathrm{lat_0}\\\\
671 {\\bf c}_i &=\\arccos[\\; \\cos(\\Delta {\\bf{a}}_i)
672 \\cos(\\mathrm{b}) + |\\Delta \\gamma_i| \\,
673 \\sin(\\Delta {\\bf a}_i)
674 \\sin(\\mathrm{b})\\; ] \\\\
675 \\mathrm{lat}_i &= \\frac{180}{\\pi}
676 \\left(\\frac{\\pi}{2} - {\\bf c}_i \\right)
678 .. math::
680 \\alpha_i &= \\arcsin \\left[ \\; \\frac{ \\sin(\\Delta {\\bf a}_i )
681 \\sin(|\\Delta \\gamma_i|)}{\\sin({\\bf c}_i)}\\;
682 \\right] \\\\
683 \\alpha_i &= \\begin{cases}
684 \\alpha_i, &\\text{if} \\; \\cos(\\Delta {\\bf a}_i) -
685 \\cos(\\mathrm{b}) \\cos({\\bf{c}}_i) > 0, \\;
686 \\text{else} \\\\
687 \\pi - \\alpha_i, & \\text{if} \\; \\alpha_i > 0,\\;
688 \\text{else}\\\\
689 -\\pi - \\alpha_i, & \\text{if} \\; \\alpha_i < 0.
690 \\end{cases} \\\\
691 \\mathrm{lon}_i &= \\mathrm{lon_0} +
692 \\frac{180}{\\pi} \\,
693 \\frac{\\Delta \\gamma_i }{|\\Delta \\gamma_i|}
694 \\cdot \\alpha_i
695 \\text{, with $\\alpha_i \\in [-\\pi,\\pi]$}
697 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
698 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
699 :param distance_rad: Distances from origin in [rad].
700 :param azimuth_rad: Azimuth from origin in [rad].
701 :type distance_rad: :py:class:`numpy.ndarray`, ``(N)``
702 :type azimuth_rad: :py:class:`numpy.ndarray`, ``(N)``
703 :type lat0: float
704 :type lon0: float
706 :return: Array with latitudes and longitudes in [deg].
707 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
708 '''
710 a = distance_rad
711 gamma = azimuth_rad
713 b = math.pi/2.-lat0*d2r
715 alphasign = 1.
716 alphasign = num.where(gamma < 0, -1., 1.)
717 gamma = num.abs(gamma)
719 c = num.arccos(clip(
720 num.cos(a)*num.cos(b)+num.sin(a)*num.sin(b)*num.cos(gamma), -1., 1.))
722 alpha = num.arcsin(clip(
723 num.sin(a)*num.sin(gamma)/num.sin(c), -1., 1.))
725 alpha = num.where(
726 num.cos(a)-num.cos(b)*num.cos(c) < 0,
727 num.where(alpha > 0, math.pi-alpha, -math.pi-alpha),
728 alpha)
730 lat = r2d * (math.pi/2. - c)
731 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
733 return lat, lon
736def crosstrack_distance(lat_begin, lon_begin, lat_end, lon_end,
737 lat_point, lon_point):
738 '''Calculate distance of a point to a great-circle path.
740 The sign of the results shows side of the path the point is on. Negative
741 distance is right of the path, positive left.
743 .. math ::
745 d_{xt} = \\arcsin \\left \\sin \\left \\Delta_{13} \\right *
746 \\sin \\left \\gamma_{13} - \\gamma_{12} \\right \\right
748 :param lat_begin: Latitude of the great circle start in [deg].
749 :param lon_begin: Longitude of the great circle start in [deg].
750 :param lat_end: Latitude of the great circle end in [deg].
751 :param lon_end: Longitude of the great circle end in [deg].
752 :param lat_point: Latitude of the point in [deg].
753 :param lon_point: Longitude of the point in [deg].
754 :type lat_begin: float
755 :type lon_begin: float
756 :type lat_end: float
757 :type lon_end: float
758 :type lat_point: float
759 :type lon_point: float
761 :return: Distance of the point to the great-circle path in [deg].
762 :rtype: float
763 '''
764 start = Loc(lat_begin, lon_begin)
765 end = Loc(lat_end, lon_end)
766 point = Loc(lat_point, lon_point)
768 dist_ang = math.acos(cosdelta(start, point))
769 azi_point = azimuth(start, point) * d2r
770 azi_end = azimuth(start, end) * d2r
772 return math.asin(math.sin(dist_ang) * math.sin(azi_point - azi_end)) * r2d
775def alongtrack_distance(lat_begin, lon_begin, lat_end, lon_end,
776 lat_point, lon_point):
777 '''Calculate distance of a point along a great-circle path in [deg].
779 Distance is relative to the beginning of the path.
781 .. math ::
783 d_{At} = \\arccos \\left \\frac{\\cos \\left \\Delta_{13} \\right}
784 {\\cos \\left \\Delta_{xt}} \\right \\right}
786 :param lat_begin: Latitude of the great circle start in [deg].
787 :param lon_begin: Longitude of the great circle start in [deg].
788 :param lat_end: Latitude of the great circle end in [deg].
789 :param lon_end: Longitude of the great circle end in [deg].
790 :param lat_point: Latitude of the point in [deg].
791 :param lon_point: Longitude of the point in [deg].
792 :type lat_begin: float
793 :type lon_begin: float
794 :type lat_end: float
795 :type lon_end: float
796 :type lat_point: float
797 :type lon_point: float
799 :return: Distance of the point along the great-circle path in [deg].
800 :rtype: float
801 '''
802 start = Loc(lat_begin, lon_begin)
803 point = Loc(lat_point, lon_point)
804 cos_dist_ang = cosdelta(start, point)
805 dist_rad = crosstrack_distance(
806 lat_begin, lon_begin, lat_end, lon_end, lat_point, lon_point) * d2r
807 return math.acos(cos_dist_ang / math.cos(dist_rad)) * r2d
810def alongtrack_distance_m(lat_begin, lon_begin, lat_end, lon_end,
811 lat_point, lon_point):
812 '''Calculate distance of a point along a great-circle path in [m].
814 Distance is relative to the beginning of the path.
816 .. math ::
818 d_{At} = \\arccos \\left \\frac{\\cos \\left \\Delta_{13} \\right }
819 { \\cos \\left \\Delta_{xt} \\right \\right}
821 :param lat_begin: Latitude of the great circle start in [deg].
822 :param lon_begin: Longitude of the great circle start in [deg].
823 :param lat_end: Latitude of the great circle end in [deg].
824 :param lon_end: Longitude of the great circle end in [deg].
825 :param lat_point: Latitude of the point in [deg].
826 :param lon_point: Longitude of the point in [deg].
827 :type lat_begin: float
828 :type lon_begin: float
829 :type lat_end: float
830 :type lon_end: float
831 :type lat_point: float
832 :type lon_point: float
834 :return: Distance of the point along the great-circle path in [m].
835 :rtype: float
836 '''
837 start = Loc(lat_begin, lon_begin)
838 end = Loc(lat_end, lon_end)
839 azi_end = azimuth(start, end)
840 dist_deg = alongtrack_distance(
841 lat_begin, lon_begin, lat_end, lon_end,
842 lat_point, lon_point)
843 along_point = Loc(
844 *azidist_to_latlon(lat_begin, lon_begin, azi_end, dist_deg))
846 return distance_accurate50m(start, along_point)
849def ne_to_latlon_alternative_method(lat0, lon0, north_m, east_m):
850 '''
851 Transform local cartesian coordinates to latitude and longitude.
853 Like :py:func:`pyrocko.orthodrome.ne_to_latlon`,
854 but this method (implementation below), although it should be numerically
855 more stable, suffers problems at points which are *across the pole*
856 as seen from the cartesian origin.
858 .. math::
860 \\text{distance change:}\\; \\Delta {{\\bf a}_i} &=
861 \\sqrt{{\\bf{y}}^2_i + {\\bf{x}}^2_i }/ \\mathrm{R_E},\\\\
862 \\text{azimuth change:}\\; \\Delta {\\bf \\gamma}_i &=
863 \\arctan( {\\bf x}_i {\\bf y}_i). \\\\
864 \\mathrm{b} &=
865 \\frac{\\pi}{2} -\\frac{\\pi}{180} \\;\\mathrm{lat_0}\\\\
867 .. math::
869 {{\\bf z}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i -
870 \\mathrm{b}}{2} \\right)}
871 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
872 {{\\bf n}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i +
873 \\mathrm{b}}{2} \\right)}
874 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
875 {{\\bf z}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i -
876 \\mathrm{b}}{2} \\right)}
877 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
878 {{\\bf n}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i +
879 \\mathrm{b}}{2} \\right)}
880 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
881 {{\\bf t}_1}_i &= \\arctan{\\left( \\frac{{{\\bf z}_1}_i}
882 {{{\\bf n}_1}_i} \\right) }\\\\
883 {{\\bf t}_2}_i &= \\arctan{\\left( \\frac{{{\\bf z}_2}_i}
884 {{{\\bf n}_2}_i} \\right) } \\\\[0.5cm]
885 c &= \\begin{cases}
886 2 \\cdot \\arccos \\left( {{\\bf z}_1}_i / \\sin({{\\bf t}_1}_i)
887 \\right),\\; \\text{if }
888 |\\sin({{\\bf t}_1}_i)| >
889 |\\sin({{\\bf t}_2}_i)|,\\; \\text{else} \\\\
890 2 \\cdot \\arcsin{\\left( {{\\bf z}_2}_i /
891 \\sin({{\\bf t}_2}_i) \\right)}.
892 \\end{cases}\\\\
894 .. math::
896 {\\bf {lat}}_i &= \\frac{180}{ \\pi } \\left( \\frac{\\pi}{2}
897 - {\\bf {c}}_i \\right) \\\\
898 {\\bf {lon}}_i &= {\\bf {lon}}_0 + \\frac{180}{ \\pi }
899 \\frac{\\gamma_i}{|\\gamma_i|},
900 \\text{ with}\\; \\gamma \\in [-\\pi,\\pi]
902 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
903 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
904 :param north_m: Northing distances from origin in [m].
905 :param east_m: Easting distances from origin in [m].
906 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
907 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
908 :type lat0: float
909 :type lon0: float
911 :return: Array with latitudes and longitudes in [deg].
912 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
913 '''
915 b = math.pi/2.-lat0*d2r
916 a = num.sqrt(north_m**2+east_m**2)/earthradius
918 gamma = num.arctan2(east_m, north_m)
919 alphasign = 1.
920 alphasign = num.where(gamma < 0., -1., 1.)
921 gamma = num.abs(gamma)
923 z1 = num.cos((a-b)/2.)*num.cos(gamma/2.)
924 n1 = num.cos((a+b)/2.)*num.sin(gamma/2.)
925 z2 = num.sin((a-b)/2.)*num.cos(gamma/2.)
926 n2 = num.sin((a+b)/2.)*num.sin(gamma/2.)
927 t1 = num.arctan2(z1, n1)
928 t2 = num.arctan2(z2, n2)
930 alpha = t1 + t2
932 sin_t1 = num.sin(t1)
933 sin_t2 = num.sin(t2)
934 c = num.where(
935 num.abs(sin_t1) > num.abs(sin_t2),
936 num.arccos(z1/sin_t1)*2.,
937 num.arcsin(z2/sin_t2)*2.)
939 lat = r2d * (math.pi/2. - c)
940 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
941 return lat, lon
944def latlon_to_ne(*args):
945 '''
946 Relative cartesian coordinates with respect to a reference location.
948 For two locations, a reference location A and another location B, given in
949 geographical coordinates in degrees, the corresponding cartesian
950 coordinates are calculated.
951 Assisting functions are :py:func:`pyrocko.orthodrome.azimuth` and
952 :py:func:`pyrocko.orthodrome.distance_accurate50m`.
954 .. math::
956 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
957 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
958 \\mathrm{)}\\\\[0.3cm]
960 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180}
961 \\varphi_{\\mathrm{azi},AB} )\\\\
962 e &= D_{AB} \\cdot
963 \\sin( \\frac{\\pi }{180} \\varphi_{\\mathrm{azi},AB})
965 :param refloc: Location reference point.
966 :type refloc: :py:class:`pyrocko.orthodrome.Loc`
967 :param loc: Location of interest.
968 :type loc: :py:class:`pyrocko.orthodrome.Loc`
970 :return: Northing and easting from refloc to location in [m].
971 :rtype: tuple[float, float]
973 '''
975 azi = azimuth(*args)
976 dist = distance_accurate50m(*args)
977 n, e = math.cos(azi*d2r)*dist, math.sin(azi*d2r)*dist
978 return n, e
981def latlon_to_ne_numpy(lat0, lon0, lat, lon):
982 '''
983 Relative cartesian coordinates with respect to a reference location.
985 For two locations, a reference location (``lat0``, ``lon0``) and another
986 location B, given in geographical coordinates in degrees,
987 the corresponding cartesian coordinates are calculated.
988 Assisting functions are :py:func:`azimuth`
989 and :py:func:`distance_accurate50m`.
991 :param lat0: Latitude of the reference location in [deg].
992 :param lon0: Longitude of the reference location in [deg].
993 :param lat: Latitude of the absolute location in [deg].
994 :param lon: Longitude of the absolute location in [deg].
996 :return: ``(n, e)``: relative north and east positions in [m].
997 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
999 Implemented formulations:
1001 .. math::
1003 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
1004 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
1005 \\mathrm{)}\\\\[0.3cm]
1007 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180} \\varphi_{
1008 \\mathrm{azi},AB} )\\\\
1009 e &= D_{AB} \\cdot \\sin( \\frac{\\pi }{180} \\varphi_{
1010 \\mathrm{azi},AB} )
1011 '''
1013 azi = azimuth_numpy(lat0, lon0, lat, lon)
1014 dist = distance_accurate50m_numpy(lat0, lon0, lat, lon)
1015 n = num.cos(azi*d2r)*dist
1016 e = num.sin(azi*d2r)*dist
1017 return n, e
1020_wgs84 = None
1023def get_wgs84():
1024 global _wgs84
1025 if _wgs84 is None:
1026 from geographiclib.geodesic import Geodesic
1027 _wgs84 = Geodesic.WGS84
1029 return _wgs84
1032def amap(n):
1034 def wrap(f):
1035 if n == 1:
1036 @wraps(f)
1037 def func(*args):
1038 it = num.nditer(args + (None,))
1039 for ops in it:
1040 ops[-1][...] = f(*ops[:-1])
1042 return it.operands[-1]
1043 elif n == 2:
1044 @wraps(f)
1045 def func(*args):
1046 it = num.nditer(args + (None, None))
1047 for ops in it:
1048 ops[-2][...], ops[-1][...] = f(*ops[:-2])
1050 return it.operands[-2], it.operands[-1]
1051 else:
1052 raise ValueError("Cannot wrap %s" % f.__qualname__)
1054 return func
1055 return wrap
1058@amap(2)
1059def ne_to_latlon2(lat0, lon0, north_m, east_m):
1060 wgs84 = get_wgs84()
1061 az = num.arctan2(east_m, north_m)*r2d
1062 dist = num.sqrt(east_m**2 + north_m**2)
1063 x = wgs84.Direct(lat0, lon0, az, dist)
1064 return x['lat2'], x['lon2']
1067@amap(2)
1068def latlon_to_ne2(lat0, lon0, lat1, lon1):
1069 wgs84 = get_wgs84()
1070 x = wgs84.Inverse(lat0, lon0, lat1, lon1)
1071 dist = x['s12']
1072 az = x['azi1']
1073 n = num.cos(az*d2r)*dist
1074 e = num.sin(az*d2r)*dist
1075 return n, e
1078@amap(1)
1079def distance_accurate15nm(lat1, lon1, lat2, lon2):
1080 wgs84 = get_wgs84()
1081 return wgs84.Inverse(lat1, lon1, lat2, lon2)['s12']
1084def positive_region(region):
1085 '''
1086 Normalize parameterization of a rectangular geographical region.
1088 :param region: ``(west, east, south, north)`` in [deg].
1089 :type region: tuple of float
1091 :returns: ``(west, east, south, north)``, where ``west <= east`` and
1092 where ``west`` and ``east`` are in the range ``[-180., 180.+360.]``.
1093 :rtype: tuple of float
1094 '''
1095 west, east, south, north = [float(x) for x in region]
1097 assert -180. - 360. <= west < 180.
1098 assert -180. < east <= 180. + 360.
1099 assert -90. <= south < 90.
1100 assert -90. < north <= 90.
1102 if east < west:
1103 east += 360.
1105 if west < -180.:
1106 west += 360.
1107 east += 360.
1109 return (west, east, south, north)
1112def points_in_region(p, region):
1113 '''
1114 Check what points are contained in a rectangular geographical region.
1116 :param p: ``(lat, lon)`` pairs in [deg].
1117 :type p: :py:class:`numpy.ndarray` ``(N, 2)``
1118 :param region: ``(west, east, south, north)`` region boundaries in [deg].
1119 :type region: tuple of float
1121 :returns: Mask, returning ``True`` for each point within the region.
1122 :rtype: :py:class:`numpy.ndarray` of bool, shape ``(N)``
1123 '''
1125 w, e, s, n = positive_region(region)
1126 return num.logical_and(
1127 num.logical_and(s <= p[:, 0], p[:, 0] <= n),
1128 num.logical_or(
1129 num.logical_and(w <= p[:, 1], p[:, 1] <= e),
1130 num.logical_and(w-360. <= p[:, 1], p[:, 1] <= e-360.)))
1133def point_in_region(p, region):
1134 '''
1135 Check if a point is contained in a rectangular geographical region.
1137 :param p: ``(lat, lon)`` in [deg].
1138 :type p: tuple of float
1139 :param region: ``(west, east, south, north)`` region boundaries in [deg].
1140 :type region: tuple of float
1142 :returns: ``True``, if point is in region, else ``False``.
1143 :rtype: bool
1144 '''
1146 w, e, s, n = positive_region(region)
1147 return num.logical_and(
1148 num.logical_and(s <= p[0], p[0] <= n),
1149 num.logical_or(
1150 num.logical_and(w <= p[1], p[1] <= e),
1151 num.logical_and(w-360. <= p[1], p[1] <= e-360.)))
1154def radius_to_region(lat, lon, radius):
1155 '''
1156 Get a rectangular region which fully contains a given circular region.
1158 :param lat: Latitude of the center point of circular region in [deg].
1159 :type lat: float
1160 :param lon: Longitude of the center point of circular region in [deg].
1161 :type lon: float
1162 :param radius: Radius of circular region in [m].
1163 :type radius: float
1165 :returns: Rectangular region as ``(east, west, south, north)`` in [deg] or
1166 ``None``.
1167 :rtype: tuple[float, float, float, float]
1168 '''
1169 radius_deg = radius * m2d
1170 if radius_deg < 45.:
1171 lat_min = max(-90., lat - radius_deg)
1172 lat_max = min(90., lat + radius_deg)
1173 absmaxlat = max(abs(lat_min), abs(lat_max))
1174 if absmaxlat > 89:
1175 lon_min = -180.
1176 lon_max = 180.
1177 else:
1178 lon_min = max(
1179 -180. - 360.,
1180 lon - radius_deg / math.cos(absmaxlat*d2r))
1181 lon_max = min(
1182 180. + 360.,
1183 lon + radius_deg / math.cos(absmaxlat*d2r))
1185 lon_min, lon_max, lat_min, lat_max = positive_region(
1186 (lon_min, lon_max, lat_min, lat_max))
1188 return lon_min, lon_max, lat_min, lat_max
1190 else:
1191 return None
1194def geographic_midpoint(lats, lons, weights=None):
1195 '''
1196 Calculate geographic midpoints by finding the center of gravity.
1198 This method suffers from instabilities if points are centered around the
1199 poles.
1201 :param lats: Latitudes in [deg].
1202 :param lons: Longitudes in [deg].
1203 :param weights: Weighting factors.
1204 :type lats: :py:class:`numpy.ndarray`, ``(N)``
1205 :type lons: :py:class:`numpy.ndarray`, ``(N)``
1206 :type weights: optional, :py:class:`numpy.ndarray`, ``(N)``
1208 :return: Latitudes and longitudes of the midpoints in [deg].
1209 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
1210 '''
1211 if not weights:
1212 weights = num.ones(len(lats))
1214 total_weigth = num.sum(weights)
1215 weights /= total_weigth
1216 lats = lats * d2r
1217 lons = lons * d2r
1218 x = num.sum(num.cos(lats) * num.cos(lons) * weights)
1219 y = num.sum(num.cos(lats) * num.sin(lons) * weights)
1220 z = num.sum(num.sin(lats) * weights)
1222 lon = num.arctan2(y, x)
1223 hyp = num.sqrt(x**2 + y**2)
1224 lat = num.arctan2(z, hyp)
1226 return lat/d2r, lon/d2r
1229def geographic_midpoint_locations(locations, weights=None):
1230 coords = num.array([loc.effective_latlon
1231 for loc in locations])
1232 return geographic_midpoint(coords[:, 0], coords[:, 1], weights)
1235def geodetic_to_ecef(lat, lon, alt):
1236 '''
1237 Convert geodetic coordinates to Earth-Centered, Earth-Fixed (ECEF)
1238 Cartesian coordinates. [#1]_ [#2]_
1240 :param lat: Geodetic latitude in [deg].
1241 :param lon: Geodetic longitude in [deg].
1242 :param alt: Geodetic altitude (height) in [m] (positive for points outside
1243 the geoid).
1244 :type lat: float
1245 :type lon: float
1246 :type alt: float
1248 :return: ECEF Cartesian coordinates (X, Y, Z) in [m].
1249 :rtype: tuple[float, float, float]
1251 .. [#1] https://en.wikipedia.org/wiki/ECEF
1252 .. [#2] https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1253 #From_geodetic_to_ECEF_coordinates
1254 '''
1256 f = earth_oblateness
1257 a = earthradius_equator
1258 e2 = 2*f - f**2
1260 lat, lon = num.radians(lat), num.radians(lon)
1261 # Normal (plumb line)
1262 N = a / num.sqrt(1.0 - (e2 * num.sin(lat)**2))
1264 X = (N+alt) * num.cos(lat) * num.cos(lon)
1265 Y = (N+alt) * num.cos(lat) * num.sin(lon)
1266 Z = (N*(1.0-e2) + alt) * num.sin(lat)
1268 return (X, Y, Z)
1271def ecef_to_geodetic(X, Y, Z):
1272 '''
1273 Convert Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates to
1274 geodetic coordinates (Ferrari's solution).
1276 :param X, Y, Z: Cartesian coordinates in ECEF system in [m].
1277 :type X, Y, Z: float
1279 :return: Geodetic coordinates (lat, lon, alt). Latitude and longitude are
1280 in [deg] and altitude is in [m]
1281 (positive for points outside the geoid).
1282 :rtype: tuple, float
1284 .. seealso ::
1285 https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1286 #The_application_of_Ferrari.27s_solution
1287 '''
1288 f = earth_oblateness
1289 a = earthradius_equator
1290 b = a * (1. - f)
1291 e2 = 2.*f - f**2
1293 # usefull
1294 a2 = a**2
1295 b2 = b**2
1296 e4 = e2**2
1297 X2 = X**2
1298 Y2 = Y**2
1299 Z2 = Z**2
1301 r = num.sqrt(X2 + Y2)
1302 r2 = r**2
1304 e_prime2 = (a2 - b2)/b2
1305 E2 = a2 - b2
1306 F = 54. * b2 * Z2
1307 G = r2 + (1.-e2)*Z2 - (e2*E2)
1308 C = (e4 * F * r2) / (G**3)
1309 S = cbrt(1. + C + num.sqrt(C**2 + 2.*C))
1310 P = F / (3. * (S + 1./S + 1.)**2 * G**2)
1311 Q = num.sqrt(1. + (2.*e4*P))
1313 dum1 = -(P*e2*r) / (1.+Q)
1314 dum2 = 0.5 * a2 * (1. + 1./Q)
1315 dum3 = (P * (1.-e2) * Z2) / (Q * (1.+Q))
1316 dum4 = 0.5 * P * r2
1317 r0 = dum1 + num.sqrt(dum2 - dum3 - dum4)
1319 U = num.sqrt((r - e2*r0)**2 + Z2)
1320 V = num.sqrt((r - e2*r0)**2 + (1.-e2)*Z2)
1321 Z0 = (b2*Z) / (a*V)
1323 alt = U * (1. - (b2 / (a*V)))
1324 lat = num.arctan((Z + e_prime2 * Z0)/r)
1325 lon = num.arctan2(Y, X)
1327 return (lat*r2d, lon*r2d, alt)
1330class Farside(Exception):
1331 pass
1334def latlon_to_xyz(latlons):
1335 if latlons.ndim == 1:
1336 return latlon_to_xyz(latlons[num.newaxis, :])[0]
1338 points = num.zeros((latlons.shape[0], 3))
1339 lats = latlons[:, 0]
1340 lons = latlons[:, 1]
1341 points[:, 0] = num.cos(lats*d2r) * num.cos(lons*d2r)
1342 points[:, 1] = num.cos(lats*d2r) * num.sin(lons*d2r)
1343 points[:, 2] = num.sin(lats*d2r)
1344 return points
1347def xyz_to_latlon(xyz):
1348 if xyz.ndim == 1:
1349 return xyz_to_latlon(xyz[num.newaxis, :])[0]
1351 latlons = num.zeros((xyz.shape[0], 2))
1352 latlons[:, 0] = num.arctan2(
1353 xyz[:, 2], num.sqrt(xyz[:, 0]**2 + xyz[:, 1]**2)) * r2d
1354 latlons[:, 1] = num.arctan2(
1355 xyz[:, 1], xyz[:, 0]) * r2d
1356 return latlons
1359def rot_to_00(lat, lon):
1360 rot0 = euler_to_matrix(0., -90.*d2r, 0.).A
1361 rot1 = euler_to_matrix(-d2r*lat, 0., -d2r*lon).A
1362 return num.dot(rot0.T, num.dot(rot1, rot0)).T
1365def distances3d(a, b):
1366 return num.sqrt(num.sum((a-b)**2, axis=a.ndim-1))
1369def circulation(points2):
1370 return num.sum(
1371 (points2[1:, 0] - points2[:-1, 0])
1372 * (points2[1:, 1] + points2[:-1, 1]))
1375def stereographic(points):
1376 dists = distances3d(points[1:, :], points[:-1, :])
1377 if dists.size > 0:
1378 maxdist = num.max(dists)
1379 cutoff = maxdist**2 / 2.
1380 else:
1381 cutoff = 1.0e-5
1383 points = points.copy()
1384 if num.any(points[:, 0] < -1. + cutoff):
1385 raise Farside()
1387 points_out = points[:, 1:].copy()
1388 factor = 1.0 / (1.0 + points[:, 0])
1389 points_out *= factor[:, num.newaxis]
1391 return points_out
1394def stereographic_poly(points):
1395 dists = distances3d(points[1:, :], points[:-1, :])
1396 if dists.size > 0:
1397 maxdist = num.max(dists)
1398 cutoff = maxdist**2 / 2.
1399 else:
1400 cutoff = 1.0e-5
1402 points = points.copy()
1403 if num.any(points[:, 0] < -1. + cutoff):
1404 raise Farside()
1406 points_out = points[:, 1:].copy()
1407 factor = 1.0 / (1.0 + points[:, 0])
1408 points_out *= factor[:, num.newaxis]
1410 if circulation(points_out) >= 0:
1411 raise Farside()
1413 return points_out
1416def gnomonic_x(points, cutoff=0.01):
1417 points_out = points[:, 1:].copy()
1418 if num.any(points[:, 0] < cutoff):
1419 raise Farside()
1421 factor = 1.0 / points[:, 0]
1422 points_out *= factor[:, num.newaxis]
1423 return points_out
1426def cneg(i, x):
1427 if i == 1:
1428 return x
1429 else:
1430 return num.logical_not(x)
1433def contains_points(polygon, points):
1434 '''
1435 Test which points are inside polygon on a sphere.
1437 The inside of the polygon is defined as the area which is to the left hand
1438 side of an observer walking the polygon line, points in order, on the
1439 sphere. Lines between the polygon points are treated as great circle paths.
1440 The polygon may be arbitrarily complex, as long as it does not have any
1441 crossings or thin parts with zero width. The polygon may contain the poles
1442 and is allowed to wrap around the sphere multiple times.
1444 The algorithm works by consecutive cutting of the polygon into (almost)
1445 hemispheres and subsequent Gnomonic projections to perform the
1446 point-in-polygon tests on a 2D plane.
1448 :param polygon: Point coordinates defining the polygon [deg].
1449 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1450 0=lat, 1=lon
1451 :param points: Coordinates of points to test [deg].
1452 :type points: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1453 0=lat, 1=lon
1455 :returns: Boolean mask array.
1456 :rtype: :py:class:`numpy.ndarray` of shape ``(N,)``.
1457 '''
1459 and_ = num.logical_and
1461 points_xyz = latlon_to_xyz(points)
1462 mask_x = 0. <= points_xyz[:, 0]
1463 mask_y = 0. <= points_xyz[:, 1]
1464 mask_z = 0. <= points_xyz[:, 2]
1466 result = num.zeros(points.shape[0], dtype=int)
1468 for ix in [-1, 1]:
1469 for iy in [-1, 1]:
1470 for iz in [-1, 1]:
1471 mask = and_(
1472 and_(cneg(ix, mask_x), cneg(iy, mask_y)),
1473 cneg(iz, mask_z))
1475 center_xyz = num.array([ix, iy, iz], dtype=float)
1477 lat, lon = xyz_to_latlon(center_xyz)
1478 rot = rot_to_00(lat, lon)
1480 points_rot_xyz = num.dot(rot, points_xyz[mask, :].T).T
1481 points_rot_pro = gnomonic_x(points_rot_xyz)
1483 offset = 0.01
1485 poly_xyz = latlon_to_xyz(polygon)
1486 poly_rot_xyz = num.dot(rot, poly_xyz.T).T
1487 poly_rot_xyz[:, 0] -= offset
1488 groups = spoly_cut([poly_rot_xyz], axis=0)
1490 for poly_rot_group_xyz in groups[1]:
1491 poly_rot_group_xyz[:, 0] += offset
1493 poly_rot_group_pro = gnomonic_x(
1494 poly_rot_group_xyz)
1496 if circulation(poly_rot_group_pro) > 0:
1497 result[mask] += path_contains_points(
1498 poly_rot_group_pro, points_rot_pro)
1499 else:
1500 result[mask] -= path_contains_points(
1501 poly_rot_group_pro, points_rot_pro)
1503 return result.astype(num.bool)
1506def contains_point(polygon, point):
1507 '''
1508 Test if point is inside polygon on a sphere.
1510 Convenience wrapper to :py:func:`contains_points` to test a single point.
1512 :param polygon: Point coordinates defining the polygon [deg].
1513 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1514 0=lat, 1=lon
1515 :param point: Coordinates ``(lat, lon)`` of point to test [deg].
1516 :type point: tuple of float
1518 :returns: ``True``, if point is located within polygon, else ``False``.
1519 :rtype: bool
1520 '''
1522 return bool(
1523 contains_points(polygon, num.asarray(point)[num.newaxis, :])[0])