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{\\sin{\\Delta_{13}} * \\\\
746 \\sin{\\gamma_{13} - \\gamma_{{12}}}}
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{\\cos{\\Delta_{13}} / \\cos{\\Delta_xt{13}}}
785 :param lat_begin: Latitude of the great circle start in [deg].
786 :param lon_begin: Longitude of the great circle start in [deg].
787 :param lat_end: Latitude of the great circle end in [deg].
788 :param lon_end: Longitude of the great circle end in [deg].
789 :param lat_point: Latitude of the point in [deg].
790 :param lon_point: Longitude of the point in [deg].
791 :type lat_begin: float
792 :type lon_begin: float
793 :type lat_end: float
794 :type lon_end: float
795 :type lat_point: float
796 :type lon_point: float
798 :return: Distance of the point along the great-circle path in [deg].
799 :rtype: float
800 '''
801 start = Loc(lat_begin, lon_begin)
802 point = Loc(lat_point, lon_point)
803 cos_dist_ang = cosdelta(start, point)
804 dist_rad = crosstrack_distance(
805 lat_begin, lon_begin, lat_end, lon_end, lat_point, lon_point) * d2r
806 return math.acos(cos_dist_ang / math.cos(dist_rad)) * r2d
809def alongtrack_distance_m(lat_begin, lon_begin, lat_end, lon_end,
810 lat_point, lon_point):
811 '''Calculate distance of a point along a great-circle path in [m].
813 Distance is relative to the beginning of the path.
815 .. math ::
817 d_{At} = \\arccos{\\cos{\\Delta_{13}} / \\cos{\\Delta_xt{13}}}
819 :param lat_begin: Latitude of the great circle start in [deg].
820 :param lon_begin: Longitude of the great circle start in [deg].
821 :param lat_end: Latitude of the great circle end in [deg].
822 :param lon_end: Longitude of the great circle end in [deg].
823 :param lat_point: Latitude of the point in [deg].
824 :param lon_point: Longitude of the point in [deg].
825 :type lat_begin: float
826 :type lon_begin: float
827 :type lat_end: float
828 :type lon_end: float
829 :type lat_point: float
830 :type lon_point: float
832 :return: Distance of the point along the great-circle path in [m].
833 :rtype: float
834 '''
835 start = Loc(lat_begin, lon_begin)
836 end = Loc(lat_end, lon_end)
837 azi_end = azimuth(start, end)
838 dist_deg = alongtrack_distance(
839 lat_begin, lon_begin, lat_end, lon_end,
840 lat_point, lon_point)
841 along_point = Loc(
842 *azidist_to_latlon(lat_begin, lon_begin, azi_end, dist_deg))
844 return distance_accurate50m(start, along_point)
847def ne_to_latlon_alternative_method(lat0, lon0, north_m, east_m):
848 '''
849 Transform local cartesian coordinates to latitude and longitude.
851 Like :py:func:`pyrocko.orthodrome.ne_to_latlon`,
852 but this method (implementation below), although it should be numerically
853 more stable, suffers problems at points which are *across the pole*
854 as seen from the cartesian origin.
856 .. math::
858 \\text{distance change:}\\; \\Delta {{\\bf a}_i} &=
859 \\sqrt{{\\bf{y}}^2_i + {\\bf{x}}^2_i }/ \\mathrm{R_E},\\\\
860 \\text{azimuth change:}\\; \\Delta {\\bf \\gamma}_i &=
861 \\arctan( {\\bf x}_i {\\bf y}_i). \\\\
862 \\mathrm{b} &=
863 \\frac{\\pi}{2} -\\frac{\\pi}{180} \\;\\mathrm{lat_0}\\\\
865 .. math::
867 {{\\bf z}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i -
868 \\mathrm{b}}{2} \\right)}
869 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
870 {{\\bf n}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i +
871 \\mathrm{b}}{2} \\right)}
872 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
873 {{\\bf z}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i -
874 \\mathrm{b}}{2} \\right)}
875 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
876 {{\\bf n}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i +
877 \\mathrm{b}}{2} \\right)}
878 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
879 {{\\bf t}_1}_i &= \\arctan{\\left( \\frac{{{\\bf z}_1}_i}
880 {{{\\bf n}_1}_i} \\right) }\\\\
881 {{\\bf t}_2}_i &= \\arctan{\\left( \\frac{{{\\bf z}_2}_i}
882 {{{\\bf n}_2}_i} \\right) } \\\\[0.5cm]
883 c &= \\begin{cases}
884 2 \\cdot \\arccos \\left( {{\\bf z}_1}_i / \\sin({{\\bf t}_1}_i)
885 \\right),\\; \\text{if }
886 |\\sin({{\\bf t}_1}_i)| >
887 |\\sin({{\\bf t}_2}_i)|,\\; \\text{else} \\\\
888 2 \\cdot \\arcsin{\\left( {{\\bf z}_2}_i /
889 \\sin({{\\bf t}_2}_i) \\right)}.
890 \\end{cases}\\\\
892 .. math::
894 {\\bf {lat}}_i &= \\frac{180}{ \\pi } \\left( \\frac{\\pi}{2}
895 - {\\bf {c}}_i \\right) \\\\
896 {\\bf {lon}}_i &= {\\bf {lon}}_0 + \\frac{180}{ \\pi }
897 \\frac{\\gamma_i}{|\\gamma_i|},
898 \\text{ with}\\; \\gamma \\in [-\\pi,\\pi]
900 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
901 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
902 :param north_m: Northing distances from origin in [m].
903 :param east_m: Easting distances from origin in [m].
904 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
905 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
906 :type lat0: float
907 :type lon0: float
909 :return: Array with latitudes and longitudes in [deg].
910 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
911 '''
913 b = math.pi/2.-lat0*d2r
914 a = num.sqrt(north_m**2+east_m**2)/earthradius
916 gamma = num.arctan2(east_m, north_m)
917 alphasign = 1.
918 alphasign = num.where(gamma < 0., -1., 1.)
919 gamma = num.abs(gamma)
921 z1 = num.cos((a-b)/2.)*num.cos(gamma/2.)
922 n1 = num.cos((a+b)/2.)*num.sin(gamma/2.)
923 z2 = num.sin((a-b)/2.)*num.cos(gamma/2.)
924 n2 = num.sin((a+b)/2.)*num.sin(gamma/2.)
925 t1 = num.arctan2(z1, n1)
926 t2 = num.arctan2(z2, n2)
928 alpha = t1 + t2
930 sin_t1 = num.sin(t1)
931 sin_t2 = num.sin(t2)
932 c = num.where(
933 num.abs(sin_t1) > num.abs(sin_t2),
934 num.arccos(z1/sin_t1)*2.,
935 num.arcsin(z2/sin_t2)*2.)
937 lat = r2d * (math.pi/2. - c)
938 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
939 return lat, lon
942def latlon_to_ne(*args):
943 '''
944 Relative cartesian coordinates with respect to a reference location.
946 For two locations, a reference location A and another location B, given in
947 geographical coordinates in degrees, the corresponding cartesian
948 coordinates are calculated.
949 Assisting functions are :py:func:`pyrocko.orthodrome.azimuth` and
950 :py:func:`pyrocko.orthodrome.distance_accurate50m`.
952 .. math::
954 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
955 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
956 \\mathrm{)}\\\\[0.3cm]
958 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180}
959 \\varphi_{\\mathrm{azi},AB} )\\\\
960 e &= D_{AB} \\cdot
961 \\sin( \\frac{\\pi }{180} \\varphi_{\\mathrm{azi},AB})
963 :param refloc: Location reference point.
964 :type refloc: :py:class:`pyrocko.orthodrome.Loc`
965 :param loc: Location of interest.
966 :type loc: :py:class:`pyrocko.orthodrome.Loc`
968 :return: Northing and easting from refloc to location in [m].
969 :rtype: tuple[float, float]
971 '''
973 azi = azimuth(*args)
974 dist = distance_accurate50m(*args)
975 n, e = math.cos(azi*d2r)*dist, math.sin(azi*d2r)*dist
976 return n, e
979def latlon_to_ne_numpy(lat0, lon0, lat, lon):
980 '''
981 Relative cartesian coordinates with respect to a reference location.
983 For two locations, a reference location (``lat0``, ``lon0``) and another
984 location B, given in geographical coordinates in degrees,
985 the corresponding cartesian coordinates are calculated.
986 Assisting functions are :py:func:`azimuth`
987 and :py:func:`distance_accurate50m`.
989 :param lat0: Latitude of the reference location in [deg].
990 :param lon0: Longitude of the reference location in [deg].
991 :param lat: Latitude of the absolute location in [deg].
992 :param lon: Longitude of the absolute location in [deg].
994 :return: ``(n, e)``: relative north and east positions in [m].
995 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
997 Implemented formulations:
999 .. math::
1001 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
1002 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
1003 \\mathrm{)}\\\\[0.3cm]
1005 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180} \\varphi_{
1006 \\mathrm{azi},AB} )\\\\
1007 e &= D_{AB} \\cdot \\sin( \\frac{\\pi }{180} \\varphi_{
1008 \\mathrm{azi},AB} )
1009 '''
1011 azi = azimuth_numpy(lat0, lon0, lat, lon)
1012 dist = distance_accurate50m_numpy(lat0, lon0, lat, lon)
1013 n = num.cos(azi*d2r)*dist
1014 e = num.sin(azi*d2r)*dist
1015 return n, e
1018_wgs84 = None
1021def get_wgs84():
1022 global _wgs84
1023 if _wgs84 is None:
1024 from geographiclib.geodesic import Geodesic
1025 _wgs84 = Geodesic.WGS84
1027 return _wgs84
1030def amap(n):
1032 def wrap(f):
1033 if n == 1:
1034 @wraps(f)
1035 def func(*args):
1036 it = num.nditer(args + (None,))
1037 for ops in it:
1038 ops[-1][...] = f(*ops[:-1])
1040 return it.operands[-1]
1041 elif n == 2:
1042 @wraps(f)
1043 def func(*args):
1044 it = num.nditer(args + (None, None))
1045 for ops in it:
1046 ops[-2][...], ops[-1][...] = f(*ops[:-2])
1048 return it.operands[-2], it.operands[-1]
1049 else:
1050 raise ValueError("Cannot wrap %s" % f.__qualname__)
1052 return func
1053 return wrap
1056@amap(2)
1057def ne_to_latlon2(lat0, lon0, north_m, east_m):
1058 wgs84 = get_wgs84()
1059 az = num.arctan2(east_m, north_m)*r2d
1060 dist = num.sqrt(east_m**2 + north_m**2)
1061 x = wgs84.Direct(lat0, lon0, az, dist)
1062 return x['lat2'], x['lon2']
1065@amap(2)
1066def latlon_to_ne2(lat0, lon0, lat1, lon1):
1067 wgs84 = get_wgs84()
1068 x = wgs84.Inverse(lat0, lon0, lat1, lon1)
1069 dist = x['s12']
1070 az = x['azi1']
1071 n = num.cos(az*d2r)*dist
1072 e = num.sin(az*d2r)*dist
1073 return n, e
1076@amap(1)
1077def distance_accurate15nm(lat1, lon1, lat2, lon2):
1078 wgs84 = get_wgs84()
1079 return wgs84.Inverse(lat1, lon1, lat2, lon2)['s12']
1082def positive_region(region):
1083 '''
1084 Normalize parameterization of a rectangular geographical region.
1086 :param region: ``(west, east, south, north)`` in [deg].
1087 :type region: tuple of float
1089 :returns: ``(west, east, south, north)``, where ``west <= east`` and
1090 where ``west`` and ``east`` are in the range ``[-180., 180.+360.]``.
1091 :rtype: tuple of float
1092 '''
1093 west, east, south, north = [float(x) for x in region]
1095 assert -180. - 360. <= west < 180.
1096 assert -180. < east <= 180. + 360.
1097 assert -90. <= south < 90.
1098 assert -90. < north <= 90.
1100 if east < west:
1101 east += 360.
1103 if west < -180.:
1104 west += 360.
1105 east += 360.
1107 return (west, east, south, north)
1110def points_in_region(p, region):
1111 '''
1112 Check what points are contained in a rectangular geographical region.
1114 :param p: ``(lat, lon)`` pairs in [deg].
1115 :type p: :py:class:`numpy.ndarray` ``(N, 2)``
1116 :param region: ``(west, east, south, north)`` region boundaries in [deg].
1117 :type region: tuple of float
1119 :returns: Mask, returning ``True`` for each point within the region.
1120 :rtype: :py:class:`numpy.ndarray` of bool, shape ``(N)``
1121 '''
1123 w, e, s, n = positive_region(region)
1124 return num.logical_and(
1125 num.logical_and(s <= p[:, 0], p[:, 0] <= n),
1126 num.logical_or(
1127 num.logical_and(w <= p[:, 1], p[:, 1] <= e),
1128 num.logical_and(w-360. <= p[:, 1], p[:, 1] <= e-360.)))
1131def point_in_region(p, region):
1132 '''
1133 Check if a point is contained in a rectangular geographical region.
1135 :param p: ``(lat, lon)`` in [deg].
1136 :type p: tuple of float
1137 :param region: ``(west, east, south, north)`` region boundaries in [deg].
1138 :type region: tuple of float
1140 :returns: ``True``, if point is in region, else ``False``.
1141 :rtype: bool
1142 '''
1144 w, e, s, n = positive_region(region)
1145 return num.logical_and(
1146 num.logical_and(s <= p[0], p[0] <= n),
1147 num.logical_or(
1148 num.logical_and(w <= p[1], p[1] <= e),
1149 num.logical_and(w-360. <= p[1], p[1] <= e-360.)))
1152def radius_to_region(lat, lon, radius):
1153 '''
1154 Get a rectangular region which fully contains a given circular region.
1156 :param lat: Latitude of the center point of circular region in [deg].
1157 :type lat: float
1158 :param lon: Longitude of the center point of circular region in [deg].
1159 :type lon: float
1160 :param radius: Radius of circular region in [m].
1161 :type radius: float
1163 :returns: Rectangular region as ``(east, west, south, north)`` in [deg] or
1164 ``None``.
1165 :rtype: tuple[float, float, float, float]
1166 '''
1167 radius_deg = radius * m2d
1168 if radius_deg < 45.:
1169 lat_min = max(-90., lat - radius_deg)
1170 lat_max = min(90., lat + radius_deg)
1171 absmaxlat = max(abs(lat_min), abs(lat_max))
1172 if absmaxlat > 89:
1173 lon_min = -180.
1174 lon_max = 180.
1175 else:
1176 lon_min = max(
1177 -180. - 360.,
1178 lon - radius_deg / math.cos(absmaxlat*d2r))
1179 lon_max = min(
1180 180. + 360.,
1181 lon + radius_deg / math.cos(absmaxlat*d2r))
1183 lon_min, lon_max, lat_min, lat_max = positive_region(
1184 (lon_min, lon_max, lat_min, lat_max))
1186 return lon_min, lon_max, lat_min, lat_max
1188 else:
1189 return None
1192def geographic_midpoint(lats, lons, weights=None):
1193 '''
1194 Calculate geographic midpoints by finding the center of gravity.
1196 This method suffers from instabilities if points are centered around the
1197 poles.
1199 :param lats: Latitudes in [deg].
1200 :param lons: Longitudes in [deg].
1201 :param weights: Weighting factors.
1202 :type lats: :py:class:`numpy.ndarray`, ``(N)``
1203 :type lons: :py:class:`numpy.ndarray`, ``(N)``
1204 :type weights: optional, :py:class:`numpy.ndarray`, ``(N)``
1206 :return: Latitudes and longitudes of the midpoints in [deg].
1207 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
1208 '''
1209 if not weights:
1210 weights = num.ones(len(lats))
1212 total_weigth = num.sum(weights)
1213 weights /= total_weigth
1214 lats = lats * d2r
1215 lons = lons * d2r
1216 x = num.sum(num.cos(lats) * num.cos(lons) * weights)
1217 y = num.sum(num.cos(lats) * num.sin(lons) * weights)
1218 z = num.sum(num.sin(lats) * weights)
1220 lon = num.arctan2(y, x)
1221 hyp = num.sqrt(x**2 + y**2)
1222 lat = num.arctan2(z, hyp)
1224 return lat/d2r, lon/d2r
1227def geographic_midpoint_locations(locations, weights=None):
1228 coords = num.array([loc.effective_latlon
1229 for loc in locations])
1230 return geographic_midpoint(coords[:, 0], coords[:, 1], weights)
1233def geodetic_to_ecef(lat, lon, alt):
1234 '''
1235 Convert geodetic coordinates to Earth-Centered, Earth-Fixed (ECEF)
1236 Cartesian coordinates. [#1]_ [#2]_
1238 :param lat: Geodetic latitude in [deg].
1239 :param lon: Geodetic longitude in [deg].
1240 :param alt: Geodetic altitude (height) in [m] (positive for points outside
1241 the geoid).
1242 :type lat: float
1243 :type lon: float
1244 :type alt: float
1246 :return: ECEF Cartesian coordinates (X, Y, Z) in [m].
1247 :rtype: tuple[float, float, float]
1249 .. [#1] https://en.wikipedia.org/wiki/ECEF
1250 .. [#2] https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1251 #From_geodetic_to_ECEF_coordinates
1252 '''
1254 f = earth_oblateness
1255 a = earthradius_equator
1256 e2 = 2*f - f**2
1258 lat, lon = num.radians(lat), num.radians(lon)
1259 # Normal (plumb line)
1260 N = a / num.sqrt(1.0 - (e2 * num.sin(lat)**2))
1262 X = (N+alt) * num.cos(lat) * num.cos(lon)
1263 Y = (N+alt) * num.cos(lat) * num.sin(lon)
1264 Z = (N*(1.0-e2) + alt) * num.sin(lat)
1266 return (X, Y, Z)
1269def ecef_to_geodetic(X, Y, Z):
1270 '''
1271 Convert Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates to
1272 geodetic coordinates (Ferrari's solution).
1274 :param X, Y, Z: Cartesian coordinates in ECEF system in [m].
1275 :type X, Y, Z: float
1277 :return: Geodetic coordinates (lat, lon, alt). Latitude and longitude are
1278 in [deg] and altitude is in [m]
1279 (positive for points outside the geoid).
1280 :rtype: tuple, float
1282 .. seealso ::
1283 https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1284 #The_application_of_Ferrari.27s_solution
1285 '''
1286 f = earth_oblateness
1287 a = earthradius_equator
1288 b = a * (1. - f)
1289 e2 = 2.*f - f**2
1291 # usefull
1292 a2 = a**2
1293 b2 = b**2
1294 e4 = e2**2
1295 X2 = X**2
1296 Y2 = Y**2
1297 Z2 = Z**2
1299 r = num.sqrt(X2 + Y2)
1300 r2 = r**2
1302 e_prime2 = (a2 - b2)/b2
1303 E2 = a2 - b2
1304 F = 54. * b2 * Z2
1305 G = r2 + (1.-e2)*Z2 - (e2*E2)
1306 C = (e4 * F * r2) / (G**3)
1307 S = cbrt(1. + C + num.sqrt(C**2 + 2.*C))
1308 P = F / (3. * (S + 1./S + 1.)**2 * G**2)
1309 Q = num.sqrt(1. + (2.*e4*P))
1311 dum1 = -(P*e2*r) / (1.+Q)
1312 dum2 = 0.5 * a2 * (1. + 1./Q)
1313 dum3 = (P * (1.-e2) * Z2) / (Q * (1.+Q))
1314 dum4 = 0.5 * P * r2
1315 r0 = dum1 + num.sqrt(dum2 - dum3 - dum4)
1317 U = num.sqrt((r - e2*r0)**2 + Z2)
1318 V = num.sqrt((r - e2*r0)**2 + (1.-e2)*Z2)
1319 Z0 = (b2*Z) / (a*V)
1321 alt = U * (1. - (b2 / (a*V)))
1322 lat = num.arctan((Z + e_prime2 * Z0)/r)
1323 lon = num.arctan2(Y, X)
1325 return (lat*r2d, lon*r2d, alt)
1328class Farside(Exception):
1329 pass
1332def latlon_to_xyz(latlons):
1333 if latlons.ndim == 1:
1334 return latlon_to_xyz(latlons[num.newaxis, :])[0]
1336 points = num.zeros((latlons.shape[0], 3))
1337 lats = latlons[:, 0]
1338 lons = latlons[:, 1]
1339 points[:, 0] = num.cos(lats*d2r) * num.cos(lons*d2r)
1340 points[:, 1] = num.cos(lats*d2r) * num.sin(lons*d2r)
1341 points[:, 2] = num.sin(lats*d2r)
1342 return points
1345def xyz_to_latlon(xyz):
1346 if xyz.ndim == 1:
1347 return xyz_to_latlon(xyz[num.newaxis, :])[0]
1349 latlons = num.zeros((xyz.shape[0], 2))
1350 latlons[:, 0] = num.arctan2(
1351 xyz[:, 2], num.sqrt(xyz[:, 0]**2 + xyz[:, 1]**2)) * r2d
1352 latlons[:, 1] = num.arctan2(
1353 xyz[:, 1], xyz[:, 0]) * r2d
1354 return latlons
1357def rot_to_00(lat, lon):
1358 rot0 = euler_to_matrix(0., -90.*d2r, 0.).A
1359 rot1 = euler_to_matrix(-d2r*lat, 0., -d2r*lon).A
1360 return num.dot(rot0.T, num.dot(rot1, rot0)).T
1363def distances3d(a, b):
1364 return num.sqrt(num.sum((a-b)**2, axis=a.ndim-1))
1367def circulation(points2):
1368 return num.sum(
1369 (points2[1:, 0] - points2[:-1, 0])
1370 * (points2[1:, 1] + points2[:-1, 1]))
1373def stereographic(points):
1374 dists = distances3d(points[1:, :], points[:-1, :])
1375 if dists.size > 0:
1376 maxdist = num.max(dists)
1377 cutoff = maxdist**2 / 2.
1378 else:
1379 cutoff = 1.0e-5
1381 points = points.copy()
1382 if num.any(points[:, 0] < -1. + cutoff):
1383 raise Farside()
1385 points_out = points[:, 1:].copy()
1386 factor = 1.0 / (1.0 + points[:, 0])
1387 points_out *= factor[:, num.newaxis]
1389 return points_out
1392def stereographic_poly(points):
1393 dists = distances3d(points[1:, :], points[:-1, :])
1394 if dists.size > 0:
1395 maxdist = num.max(dists)
1396 cutoff = maxdist**2 / 2.
1397 else:
1398 cutoff = 1.0e-5
1400 points = points.copy()
1401 if num.any(points[:, 0] < -1. + cutoff):
1402 raise Farside()
1404 points_out = points[:, 1:].copy()
1405 factor = 1.0 / (1.0 + points[:, 0])
1406 points_out *= factor[:, num.newaxis]
1408 if circulation(points_out) >= 0:
1409 raise Farside()
1411 return points_out
1414def gnomonic_x(points, cutoff=0.01):
1415 points_out = points[:, 1:].copy()
1416 if num.any(points[:, 0] < cutoff):
1417 raise Farside()
1419 factor = 1.0 / points[:, 0]
1420 points_out *= factor[:, num.newaxis]
1421 return points_out
1424def cneg(i, x):
1425 if i == 1:
1426 return x
1427 else:
1428 return num.logical_not(x)
1431def contains_points(polygon, points):
1432 '''
1433 Test which points are inside polygon on a sphere.
1435 The inside of the polygon is defined as the area which is to the left hand
1436 side of an observer walking the polygon line, points in order, on the
1437 sphere. Lines between the polygon points are treated as great circle paths.
1438 The polygon may be arbitrarily complex, as long as it does not have any
1439 crossings or thin parts with zero width. The polygon may contain the poles
1440 and is allowed to wrap around the sphere multiple times.
1442 The algorithm works by consecutive cutting of the polygon into (almost)
1443 hemispheres and subsequent Gnomonic projections to perform the
1444 point-in-polygon tests on a 2D plane.
1446 :param polygon: Point coordinates defining the polygon [deg].
1447 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1448 0=lat, 1=lon
1449 :param points: Coordinates of points to test [deg].
1450 :type points: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1451 0=lat, 1=lon
1453 :returns: Boolean mask array.
1454 :rtype: :py:class:`numpy.ndarray` of shape ``(N,)``.
1455 '''
1457 and_ = num.logical_and
1459 points_xyz = latlon_to_xyz(points)
1460 mask_x = 0. <= points_xyz[:, 0]
1461 mask_y = 0. <= points_xyz[:, 1]
1462 mask_z = 0. <= points_xyz[:, 2]
1464 result = num.zeros(points.shape[0], dtype=int)
1466 for ix in [-1, 1]:
1467 for iy in [-1, 1]:
1468 for iz in [-1, 1]:
1469 mask = and_(
1470 and_(cneg(ix, mask_x), cneg(iy, mask_y)),
1471 cneg(iz, mask_z))
1473 center_xyz = num.array([ix, iy, iz], dtype=float)
1475 lat, lon = xyz_to_latlon(center_xyz)
1476 rot = rot_to_00(lat, lon)
1478 points_rot_xyz = num.dot(rot, points_xyz[mask, :].T).T
1479 points_rot_pro = gnomonic_x(points_rot_xyz)
1481 offset = 0.01
1483 poly_xyz = latlon_to_xyz(polygon)
1484 poly_rot_xyz = num.dot(rot, poly_xyz.T).T
1485 poly_rot_xyz[:, 0] -= offset
1486 groups = spoly_cut([poly_rot_xyz], axis=0)
1488 for poly_rot_group_xyz in groups[1]:
1489 poly_rot_group_xyz[:, 0] += offset
1491 poly_rot_group_pro = gnomonic_x(
1492 poly_rot_group_xyz)
1494 if circulation(poly_rot_group_pro) > 0:
1495 result[mask] += path_contains_points(
1496 poly_rot_group_pro, points_rot_pro)
1497 else:
1498 result[mask] -= path_contains_points(
1499 poly_rot_group_pro, points_rot_pro)
1501 return result.astype(num.bool)
1504def contains_point(polygon, point):
1505 '''
1506 Test if point is inside polygon on a sphere.
1508 Convenience wrapper to :py:func:`contains_points` to test a single point.
1510 :param polygon: Point coordinates defining the polygon [deg].
1511 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1512 0=lat, 1=lon
1513 :param point: Coordinates ``(lat, lon)`` of point to test [deg].
1514 :type point: tuple of float
1516 :returns: ``True``, if point is located within polygon, else ``False``.
1517 :rtype: bool
1518 '''
1520 return bool(
1521 contains_points(polygon, num.asarray(point)[num.newaxis, :])[0])