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
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
134def cosdelta(*args):
135 '''
136 Cosine of the angular distance between two points ``a`` and ``b`` on a
137 sphere.
139 This function (find implementation below) returns the cosine of the
140 distance angle 'delta' between two points ``a`` and ``b``, coordinates of
141 which are expected to be given in geographical coordinates and in degrees.
142 For numerical stability a maximum of 1.0 is enforced.
144 .. math::
146 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
147 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
148 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
149 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\[0.5cm]
151 \\cos(\\Delta) = \\min( 1.0, \\quad \\sin( A_{\\mathrm{lat'}})
152 \\sin( B_{\\mathrm{lat'}} ) +
153 \\cos(A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
154 \\cos( B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )
156 :param a: Location point A.
157 :type a: :py:class:`pyrocko.orthodrome.Loc`
158 :param b: Location point B.
159 :type b: :py:class:`pyrocko.orthodrome.Loc`
161 :return: Cosdelta.
162 :rtype: float
163 '''
165 alat, alon, blat, blon = _latlon_pair(args)
167 return min(
168 1.0,
169 math.sin(alat*d2r) * math.sin(blat*d2r) +
170 math.cos(alat*d2r) * math.cos(blat*d2r) *
171 math.cos(d2r*(blon-alon)))
174def cosdelta_numpy(a_lats, a_lons, b_lats, b_lons):
175 '''
176 Cosine of the angular distance between two points ``a`` and ``b`` on a
177 sphere.
179 This function returns the cosines of the distance
180 angles *delta* between two points ``a`` and ``b`` given as
181 :py:class:`numpy.ndarray`.
182 The coordinates are expected to be given in geographical coordinates
183 and in degrees. For numerical stability a maximum of ``1.0`` is enforced.
185 Please find the details of the implementation in the documentation of
186 the function :py:func:`pyrocko.orthodrome.cosdelta` above.
188 :param a_lats: Latitudes in [deg] point A.
189 :param a_lons: Longitudes in [deg] point A.
190 :param b_lats: Latitudes in [deg] point B.
191 :param b_lons: Longitudes in [deg] point B.
192 :type a_lats: :py:class:`numpy.ndarray`
193 :type a_lons: :py:class:`numpy.ndarray`
194 :type b_lats: :py:class:`numpy.ndarray`
195 :type b_lons: :py:class:`numpy.ndarray`
197 :return: Cosdelta.
198 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
199 '''
200 return num.minimum(
201 1.0,
202 num.sin(a_lats*d2r) * num.sin(b_lats*d2r) +
203 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
204 num.cos(d2r*(b_lons-a_lons)))
207def azimuth(*args):
208 '''
209 Azimuth calculation.
211 This function (find implementation below) returns azimuth ...
212 between points ``a`` and ``b``, coordinates of
213 which are expected to be given in geographical coordinates and in degrees.
215 .. math::
217 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
218 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
219 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
220 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\
222 \\varphi_{\\mathrm{azi},AB} = \\frac{180}{\\pi} \\arctan \\left[
223 \\frac{
224 \\cos( A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
225 \\sin(B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )}
226 {\\sin ( B_{\\mathrm{lat'}} ) - \\sin( A_{\\mathrm{lat'}}
227 cosdelta) } \\right]
229 :param a: Location point A.
230 :type a: :py:class:`pyrocko.orthodrome.Loc`
231 :param b: Location point B.
232 :type b: :py:class:`pyrocko.orthodrome.Loc`
234 :return: Azimuth in degree
235 '''
237 alat, alon, blat, blon = _latlon_pair(args)
239 return r2d*math.atan2(
240 math.cos(alat*d2r) * math.cos(blat*d2r) *
241 math.sin(d2r*(blon-alon)),
242 math.sin(d2r*blat) - math.sin(d2r*alat) * cosdelta(
243 alat, alon, blat, blon))
246def azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta=None):
247 '''
248 Calculation of the azimuth (*track angle*) from a location A towards B.
250 This function returns azimuths (*track angles*) from locations A towards B
251 given in :py:class:`numpy.ndarray`. Coordinates are expected to be given in
252 geographical coordinates and in degrees.
254 Please find the details of the implementation in the documentation of the
255 function :py:func:`pyrocko.orthodrome.azimuth`.
257 :param a_lats: Latitudes in [deg] point A.
258 :param a_lons: Longitudes in [deg] point A.
259 :param b_lats: Latitudes in [deg] point B.
260 :param b_lons: Longitudes in [deg] point B.
261 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
262 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
263 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
264 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
266 :return: Azimuths in [deg].
267 :rtype: :py:class:`numpy.ndarray`, ``(N)``
268 '''
269 if _cosdelta is None:
270 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
272 return r2d*num.arctan2(
273 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
274 num.sin(d2r*(b_lons-a_lons)),
275 num.sin(d2r*b_lats) - num.sin(d2r*a_lats) * _cosdelta)
278def azibazi(*args, **kwargs):
279 '''
280 Azimuth and backazimuth from location A towards B and back.
282 :returns: Azimuth in [deg] from A to B, back azimuth in [deg] from B to A.
283 :rtype: tuple[float, float]
284 '''
286 alat, alon, blat, blon = _latlon_pair(args)
287 if alat == blat and alon == blon:
288 return 0., 180.
290 implementation = kwargs.get('implementation', 'c')
291 assert implementation in ('c', 'python')
292 if implementation == 'c':
293 from pyrocko import orthodrome_ext
294 return orthodrome_ext.azibazi(alat, alon, blat, blon)
296 cd = cosdelta(alat, alon, blat, blon)
297 azi = r2d*math.atan2(
298 math.cos(alat*d2r) * math.cos(blat*d2r) *
299 math.sin(d2r*(blon-alon)),
300 math.sin(d2r*blat) - math.sin(d2r*alat) * cd)
301 bazi = r2d*math.atan2(
302 math.cos(blat*d2r) * math.cos(alat*d2r) *
303 math.sin(d2r*(alon-blon)),
304 math.sin(d2r*alat) - math.sin(d2r*blat) * cd)
306 return azi, bazi
309def azibazi_numpy(a_lats, a_lons, b_lats, b_lons, implementation='c'):
310 '''
311 Azimuth and backazimuth from location A towards B and back.
313 Arguments are given as :py:class:`numpy.ndarray`.
315 :param a_lats: Latitude(s) in [deg] of point A.
316 :type a_lats: :py:class:`numpy.ndarray`
317 :param a_lons: Longitude(s) in [deg] of point A.
318 :type a_lons: :py:class:`numpy.ndarray`
319 :param b_lats: Latitude(s) in [deg] of point B.
320 :type b_lats: :py:class:`numpy.ndarray`
321 :param b_lons: Longitude(s) in [deg] of point B.
322 :type b_lons: :py:class:`numpy.ndarray`
324 :returns: Azimuth(s) in [deg] from A to B,
325 back azimuth(s) in [deg] from B to A.
326 :rtype: :py:class:`numpy.ndarray`, :py:class:`numpy.ndarray`
327 '''
329 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
330 a_lats, a_lons, b_lats, b_lons)
332 assert implementation in ('c', 'python')
333 if implementation == 'c':
334 from pyrocko import orthodrome_ext
335 return orthodrome_ext.azibazi_numpy(a_lats, a_lons, b_lats, b_lons)
337 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
338 azis = azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta)
339 bazis = azimuth_numpy(b_lats, b_lons, a_lats, a_lons, _cosdelta)
341 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
342 ii_eq = num.where(eq)[0]
343 azis[ii_eq] = 0.0
344 bazis[ii_eq] = 180.0
345 return azis, bazis
348def azidist_numpy(*args):
349 '''
350 Calculation of the azimuth (*track angle*) and the distance from locations
351 A towards B on a sphere.
353 The assisting functions used are :py:func:`pyrocko.orthodrome.cosdelta` and
354 :py:func:`pyrocko.orthodrome.azimuth`
356 :param a_lats: Latitudes in [deg] point A.
357 :param a_lons: Longitudes in [deg] point A.
358 :param b_lats: Latitudes in [deg] point B.
359 :param b_lons: Longitudes in [deg] point B.
360 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
361 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
362 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
363 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
365 :return: Azimuths in [deg], distances in [deg].
366 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
367 '''
368 _cosdelta = cosdelta_numpy(*args)
369 _azimuths = azimuth_numpy(_cosdelta=_cosdelta, *args)
370 return _azimuths, r2d*num.arccos(_cosdelta)
373def distance_accurate50m(*args, **kwargs):
374 '''
375 Accurate distance calculation based on a spheroid of rotation.
377 Function returns distance in meter between points A and B, coordinates of
378 which must be given in geographical coordinates and in degrees.
379 The returned distance should be accurate to 50 m using WGS84.
380 Values for the Earth's equator radius and the Earth's oblateness
381 (``f_oblate``) are defined in the pyrocko configuration file
382 :py:class:`pyrocko.config`.
384 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
386 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
387 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
389 .. math::
391 F = \\frac{\\pi}{180}
392 \\frac{(A_{lat} + B_{lat})}{2}, \\quad
393 G = \\frac{\\pi}{180}
394 \\frac{(A_{lat} - B_{lat})}{2}, \\quad
395 l = \\frac{\\pi}{180}
396 \\frac{(A_{lon} - B_{lon})}{2} \\quad
397 \\\\[0.5cm]
398 S = \\sin^2(G) \\cdot \\cos^2(l) +
399 \\cos^2(F) \\cdot \\sin^2(l), \\quad \\quad
400 C = \\cos^2(G) \\cdot \\cos^2(l) +
401 \\sin^2(F) \\cdot \\sin^2(l)
403 .. math::
405 w = \\arctan \\left( \\sqrt{ \\frac{S}{C}} \\right) , \\quad
406 r = \\sqrt{\\frac{S}{C} }
408 The spherical-earth distance D between A and B, can be given with:
410 .. math::
412 D_{sphere} = 2w \\cdot R_{equator}
414 The oblateness of the Earth requires some correction with
415 correction factors h1 and h2:
417 .. math::
419 h_1 = \\frac{3r - 1}{2C}, \\quad
420 h_2 = \\frac{3r +1 }{2S}\\\\[0.5cm]
422 D = D_{\\mathrm{sphere}} \\cdot [ 1 + h_1 \\,f_{\\mathrm{oblate}}
423 \\cdot \\sin^2(F)
424 \\cos^2(G) - h_2\\, f_{\\mathrm{oblate}}
425 \\cdot \\cos^2(F) \\sin^2(G)]
428 :param a: Location point A.
429 :type a: :py:class:`pyrocko.orthodrome.Loc`
430 :param b: Location point B.
431 :type b: :py:class:`pyrocko.orthodrome.Loc`
433 :return: Distance in [m].
434 :rtype: float
435 '''
437 alat, alon, blat, blon = _latlon_pair(args)
439 implementation = kwargs.get('implementation', 'c')
440 assert implementation in ('c', 'python')
441 if implementation == 'c':
442 from pyrocko import orthodrome_ext
443 return orthodrome_ext.distance_accurate50m(alat, alon, blat, blon)
445 f = (alat + blat)*d2r / 2.
446 g = (alat - blat)*d2r / 2.
447 h = (alon - blon)*d2r / 2.
449 s = math.sin(g)**2 * math.cos(h)**2 + math.cos(f)**2 * math.sin(h)**2
450 c = math.cos(g)**2 * math.cos(h)**2 + math.sin(f)**2 * math.sin(h)**2
452 w = math.atan(math.sqrt(s/c))
454 if w == 0.0:
455 return 0.0
457 r = math.sqrt(s*c)/w
458 d = 2.*w*earthradius_equator
459 h1 = (3.*r-1.)/(2.*c)
460 h2 = (3.*r+1.)/(2.*s)
462 return d * (1. +
463 earth_oblateness * h1 * math.sin(f)**2 * math.cos(g)**2 -
464 earth_oblateness * h2 * math.cos(f)**2 * math.sin(g)**2)
467def distance_accurate50m_numpy(
468 a_lats, a_lons, b_lats, b_lons, implementation='c'):
470 '''
471 Accurate distance calculation based on a spheroid of rotation.
473 Function returns distance in meter between points ``a`` and ``b``,
474 coordinates of which must be given in geographical coordinates and in
475 degrees.
476 The returned distance should be accurate to 50 m using WGS84.
477 Values for the Earth's equator radius and the Earth's oblateness
478 (``f_oblate``) are defined in the pyrocko configuration file
479 :py:class:`pyrocko.config`.
481 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
483 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
484 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
486 .. math::
488 F_i = \\frac{\\pi}{180}
489 \\frac{(a_{lat,i} + a_{lat,i})}{2}, \\quad
490 G_i = \\frac{\\pi}{180}
491 \\frac{(a_{lat,i} - b_{lat,i})}{2}, \\quad
492 l_i= \\frac{\\pi}{180}
493 \\frac{(a_{lon,i} - b_{lon,i})}{2} \\\\[0.5cm]
494 S_i = \\sin^2(G_i) \\cdot \\cos^2(l_i) +
495 \\cos^2(F_i) \\cdot \\sin^2(l_i), \\quad \\quad
496 C_i = \\cos^2(G_i) \\cdot \\cos^2(l_i) +
497 \\sin^2(F_i) \\cdot \\sin^2(l_i)
499 .. math::
501 w_i = \\arctan \\left( \\sqrt{\\frac{S_i}{C_i}} \\right), \\quad
502 r_i = \\sqrt{\\frac{S_i}{C_i} }
504 The spherical-earth distance ``D`` between ``a`` and ``b``,
505 can be given with:
507 .. math::
509 D_{\\mathrm{sphere},i} = 2w_i \\cdot R_{\\mathrm{equator}}
511 The oblateness of the Earth requires some correction with
512 correction factors ``h1`` and ``h2``:
514 .. math::
516 h_{1.i} = \\frac{3r - 1}{2C_i}, \\quad
517 h_{2,i} = \\frac{3r +1 }{2S_i}\\\\[0.5cm]
519 D_{AB,i} = D_{\\mathrm{sphere},i} \\cdot [1 + h_{1,i}
520 \\,f_{\\mathrm{oblate}}
521 \\cdot \\sin^2(F_i)
522 \\cos^2(G_i) - h_{2,i}\\, f_{\\mathrm{oblate}}
523 \\cdot \\cos^2(F_i) \\sin^2(G_i)]
525 :param a_lats: Latitudes in [deg] point A.
526 :param a_lons: Longitudes in [deg] point A.
527 :param b_lats: Latitudes in [deg] point B.
528 :param b_lons: Longitudes in [deg] point B.
529 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
530 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
531 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
532 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
534 :return: Distances in [m].
535 :rtype: :py:class:`numpy.ndarray`, ``(N)``
536 '''
538 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
539 a_lats, a_lons, b_lats, b_lons)
541 assert implementation in ('c', 'python')
542 if implementation == 'c':
543 from pyrocko import orthodrome_ext
544 return orthodrome_ext.distance_accurate50m_numpy(
545 a_lats, a_lons, b_lats, b_lons)
547 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
548 ii_neq = num.where(num.logical_not(eq))[0]
550 if num.all(eq):
551 return num.zeros_like(eq, dtype=float)
553 def extr(x):
554 if isinstance(x, num.ndarray) and x.size > 1:
555 return x[ii_neq]
556 else:
557 return x
559 a_lats = extr(a_lats)
560 a_lons = extr(a_lons)
561 b_lats = extr(b_lats)
562 b_lons = extr(b_lons)
564 f = (a_lats + b_lats)*d2r / 2.
565 g = (a_lats - b_lats)*d2r / 2.
566 h = (a_lons - b_lons)*d2r / 2.
568 s = num.sin(g)**2 * num.cos(h)**2 + num.cos(f)**2 * num.sin(h)**2
569 c = num.cos(g)**2 * num.cos(h)**2 + num.sin(f)**2 * num.sin(h)**2
571 w = num.arctan(num.sqrt(s/c))
573 r = num.sqrt(s*c)/w
575 d = 2.*w*earthradius_equator
576 h1 = (3.*r-1.)/(2.*c)
577 h2 = (3.*r+1.)/(2.*s)
579 dists = num.zeros(eq.size, dtype=float)
580 dists[ii_neq] = d * (
581 1. +
582 earth_oblateness * h1 * num.sin(f)**2 * num.cos(g)**2 -
583 earth_oblateness * h2 * num.cos(f)**2 * num.sin(g)**2)
585 return dists
588def ne_to_latlon(lat0, lon0, north_m, east_m):
589 '''
590 Transform local cartesian coordinates to latitude and longitude.
592 From east and north coordinates (``x`` and ``y`` coordinate
593 :py:class:`numpy.ndarray`) relative to a reference differences in
594 longitude and latitude are calculated, which are effectively changes in
595 azimuth and distance, respectively:
597 .. math::
599 \\text{distance change:}\\; \\Delta {\\bf{a}} &= \\sqrt{{\\bf{y}}^2 +
600 {\\bf{x}}^2 }/ \\mathrm{R_E},
602 \\text{azimuth change:}\\; \\Delta \\bf{\\gamma} &= \\arctan( \\bf{x}
603 / \\bf{y}).
605 The projection used preserves the azimuths of the input points.
607 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
608 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
609 :param north_m: Northing distances from origin in [m].
610 :param east_m: Easting distances from origin in [m].
611 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
612 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
613 :type lat0: float
614 :type lon0: float
616 :return: Array with latitudes and longitudes in [deg].
617 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
619 '''
621 a = num.sqrt(north_m**2+east_m**2)/earthradius
622 gamma = num.arctan2(east_m, north_m)
624 return azidist_to_latlon_rad(lat0, lon0, gamma, a)
627def azidist_to_latlon(lat0, lon0, azimuth_deg, distance_deg):
628 '''
629 Absolute latitudes and longitudes are calculated from relative changes.
631 Convenience wrapper to :py:func:`azidist_to_latlon_rad` with azimuth and
632 distance given in degrees.
634 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
635 :type lat0: float
636 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
637 :type lon0: float
638 :param azimuth_deg: Azimuth from origin in [deg].
639 :type azimuth_deg: :py:class:`numpy.ndarray`, ``(N)``
640 :param distance_deg: Distances from origin in [deg].
641 :type distance_deg: :py:class:`numpy.ndarray`, ``(N)``
643 :return: Array with latitudes and longitudes in [deg].
644 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
645 '''
647 return azidist_to_latlon_rad(
648 lat0, lon0, azimuth_deg/180.*num.pi, distance_deg/180.*num.pi)
651def azidist_to_latlon_rad(lat0, lon0, azimuth_rad, distance_rad):
652 '''
653 Absolute latitudes and longitudes are calculated from relative changes.
655 For numerical stability a range between of ``-1.0`` and ``1.0`` is
656 enforced for ``c`` and ``alpha``.
658 .. math::
660 \\Delta {\\bf a}_i \\; \\text{and} \\; \\Delta \\gamma_i \\;
661 \\text{are relative distances and azimuths from lat0 and lon0 for
662 \\textit{i} source points of a finite source.}
664 .. math::
666 \\mathrm{b} &= \\frac{\\pi}{2} -\\frac{\\pi}{180}\\;\\mathrm{lat_0}\\\\
667 {\\bf c}_i &=\\arccos[\\; \\cos(\\Delta {\\bf{a}}_i)
668 \\cos(\\mathrm{b}) + |\\Delta \\gamma_i| \\,
669 \\sin(\\Delta {\\bf a}_i)
670 \\sin(\\mathrm{b})\\; ] \\\\
671 \\mathrm{lat}_i &= \\frac{180}{\\pi}
672 \\left(\\frac{\\pi}{2} - {\\bf c}_i \\right)
674 .. math::
676 \\alpha_i &= \\arcsin \\left[ \\; \\frac{ \\sin(\\Delta {\\bf a}_i )
677 \\sin(|\\Delta \\gamma_i|)}{\\sin({\\bf c}_i)}\\;
678 \\right] \\\\
679 \\alpha_i &= \\begin{cases}
680 \\alpha_i, &\\text{if} \\; \\cos(\\Delta {\\bf a}_i) -
681 \\cos(\\mathrm{b}) \\cos({\\bf{c}}_i) > 0, \\;
682 \\text{else} \\\\
683 \\pi - \\alpha_i, & \\text{if} \\; \\alpha_i > 0,\\;
684 \\text{else}\\\\
685 -\\pi - \\alpha_i, & \\text{if} \\; \\alpha_i < 0.
686 \\end{cases} \\\\
687 \\mathrm{lon}_i &= \\mathrm{lon_0} +
688 \\frac{180}{\\pi} \\,
689 \\frac{\\Delta \\gamma_i }{|\\Delta \\gamma_i|}
690 \\cdot \\alpha_i
691 \\text{, with $\\alpha_i \\in [-\\pi,\\pi]$}
693 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
694 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
695 :param distance_rad: Distances from origin in [rad].
696 :param azimuth_rad: Azimuth from origin in [rad].
697 :type distance_rad: :py:class:`numpy.ndarray`, ``(N)``
698 :type azimuth_rad: :py:class:`numpy.ndarray`, ``(N)``
699 :type lat0: float
700 :type lon0: float
702 :return: Array with latitudes and longitudes in [deg].
703 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
704 '''
706 a = distance_rad
707 gamma = azimuth_rad
709 b = math.pi/2.-lat0*d2r
711 alphasign = 1.
712 alphasign = num.where(gamma < 0, -1., 1.)
713 gamma = num.abs(gamma)
715 c = num.arccos(clip(
716 num.cos(a)*num.cos(b)+num.sin(a)*num.sin(b)*num.cos(gamma), -1., 1.))
718 alpha = num.arcsin(clip(
719 num.sin(a)*num.sin(gamma)/num.sin(c), -1., 1.))
721 alpha = num.where(
722 num.cos(a)-num.cos(b)*num.cos(c) < 0,
723 num.where(alpha > 0, math.pi-alpha, -math.pi-alpha),
724 alpha)
726 lat = r2d * (math.pi/2. - c)
727 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
729 return lat, lon
732def crosstrack_distance(begin_lat, begin_lon, end_lat, end_lon,
733 point_lat, point_lon):
734 '''Calculate distance of a point to a great-circle path.
736 The sign of the results shows side of the path the point is on. Negative
737 distance is right of the path, positive left.
739 .. math ::
741 d_{xt} = \\arcsin \\left( \\sin \\left( \\Delta_{13} \\right) \\cdot
742 \\sin \\left( \\gamma_{13} - \\gamma_{12} \\right) \\right)
744 :param begin_lat: Latitude of the great circle start in [deg].
745 :param begin_lon: Longitude of the great circle start in [deg].
746 :param end_lat: Latitude of the great circle end in [deg].
747 :param end_lon: Longitude of the great circle end in [deg].
748 :param point_lat: Latitude of the point in [deg].
749 :param point_lon: Longitude of the point in [deg].
750 :type begin_lat: float
751 :type begin_lon: float
752 :type end_lat: float
753 :type end_lon: float
754 :type point_lat: float
755 :type point_lon: float
757 :return: Distance of the point to the great-circle path in [deg].
758 :rtype: float
759 '''
760 start = Loc(begin_lat, begin_lon)
761 end = Loc(end_lat, end_lon)
762 point = Loc(point_lat, point_lon)
764 dist_ang = math.acos(cosdelta(start, point))
765 azi_point = azimuth(start, point) * d2r
766 azi_end = azimuth(start, end) * d2r
768 return math.asin(math.sin(dist_ang) * math.sin(azi_point - azi_end)) * r2d
771def alongtrack_distance(begin_lat, begin_lon, end_lat, end_lon,
772 point_lat, point_lon):
773 '''Calculate distance of a point along a great-circle path in [deg].
775 Distance is relative to the beginning of the path.
777 .. math ::
779 d_{At} = \\arccos \\left( \\frac{\\cos \\left( \\Delta_{13} \\right)}
780 {\\cos \\left( \\Delta_{xt} \\right) } \\right)
782 :param begin_lat: Latitude of the great circle start in [deg].
783 :param begin_lon: Longitude of the great circle start in [deg].
784 :param end_lat: Latitude of the great circle end in [deg].
785 :param end_lon: Longitude of the great circle end in [deg].
786 :param point_lat: Latitude of the point in [deg].
787 :param point_lon: Longitude of the point in [deg].
788 :type begin_lat: float
789 :type begin_lon: float
790 :type end_lat: float
791 :type end_lon: float
792 :type point_lat: float
793 :type point_lon: float
795 :return: Distance of the point along the great-circle path in [deg].
796 :rtype: float
797 '''
798 start = Loc(begin_lat, begin_lon)
799 point = Loc(point_lat, point_lon)
800 cos_dist_ang = cosdelta(start, point)
801 dist_rad = crosstrack_distance(
802 begin_lat, begin_lon, end_lat, end_lon, point_lat, point_lon) * d2r
803 return math.acos(cos_dist_ang / math.cos(dist_rad)) * r2d
806def alongtrack_distance_m(begin_lat, begin_lon, end_lat, end_lon,
807 point_lat, point_lon):
808 '''Calculate distance of a point along a great-circle path in [m].
810 Distance is relative to the beginning of the path.
812 .. math ::
814 d_{At} = \\arccos \\left( \\frac{\\cos \\left( \\Delta_{13} \\right)}
815 {\\cos \\left( \\Delta_{xt} \\right) } \\right)
817 :param begin_lat: Latitude of the great circle start in [deg].
818 :param begin_lon: Longitude of the great circle start in [deg].
819 :param end_lat: Latitude of the great circle end in [deg].
820 :param end_lon: Longitude of the great circle end in [deg].
821 :param point_lat: Latitude of the point in [deg].
822 :param point_lon: Longitude of the point in [deg].
823 :type begin_lat: float
824 :type begin_lon: float
825 :type end_lat: float
826 :type end_lon: float
827 :type point_lat: float
828 :type point_lon: float
830 :return: Distance of the point along the great-circle path in [m].
831 :rtype: float
832 '''
833 start = Loc(begin_lat, begin_lon)
834 end = Loc(end_lat, end_lon)
835 azi_end = azimuth(start, end)
836 dist_deg = alongtrack_distance(
837 begin_lat, begin_lon, end_lat, end_lon,
838 point_lat, point_lon)
839 along_point = Loc(
840 *azidist_to_latlon(begin_lat, begin_lon, azi_end, dist_deg))
842 return distance_accurate50m(start, along_point)
845def ne_to_latlon_alternative_method(lat0, lon0, north_m, east_m):
846 '''
847 Transform local cartesian coordinates to latitude and longitude.
849 Like :py:func:`pyrocko.orthodrome.ne_to_latlon`,
850 but this method (implementation below), although it should be numerically
851 more stable, suffers problems at points which are *across the pole*
852 as seen from the cartesian origin.
854 .. math::
856 \\text{distance change:}\\; \\Delta {{\\bf a}_i} &=
857 \\sqrt{{\\bf{y}}^2_i + {\\bf{x}}^2_i }/ \\mathrm{R_E},\\\\
858 \\text{azimuth change:}\\; \\Delta {\\bf \\gamma}_i &=
859 \\arctan( {\\bf x}_i {\\bf y}_i). \\\\
860 \\mathrm{b} &=
861 \\frac{\\pi}{2} -\\frac{\\pi}{180} \\;\\mathrm{lat_0}\\\\
863 .. math::
865 {{\\bf z}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i -
866 \\mathrm{b}}{2} \\right)}
867 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
868 {{\\bf n}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i +
869 \\mathrm{b}}{2} \\right)}
870 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
871 {{\\bf z}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i -
872 \\mathrm{b}}{2} \\right)}
873 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
874 {{\\bf n}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i +
875 \\mathrm{b}}{2} \\right)}
876 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
877 {{\\bf t}_1}_i &= \\arctan{\\left( \\frac{{{\\bf z}_1}_i}
878 {{{\\bf n}_1}_i} \\right) }\\\\
879 {{\\bf t}_2}_i &= \\arctan{\\left( \\frac{{{\\bf z}_2}_i}
880 {{{\\bf n}_2}_i} \\right) } \\\\[0.5cm]
881 c &= \\begin{cases}
882 2 \\cdot \\arccos \\left( {{\\bf z}_1}_i / \\sin({{\\bf t}_1}_i)
883 \\right),\\; \\text{if }
884 |\\sin({{\\bf t}_1}_i)| >
885 |\\sin({{\\bf t}_2}_i)|,\\; \\text{else} \\\\
886 2 \\cdot \\arcsin{\\left( {{\\bf z}_2}_i /
887 \\sin({{\\bf t}_2}_i) \\right)}.
888 \\end{cases}\\\\
890 .. math::
892 {\\bf {lat}}_i &= \\frac{180}{ \\pi } \\left( \\frac{\\pi}{2}
893 - {\\bf {c}}_i \\right) \\\\
894 {\\bf {lon}}_i &= {\\bf {lon}}_0 + \\frac{180}{ \\pi }
895 \\frac{\\gamma_i}{|\\gamma_i|},
896 \\text{ with}\\; \\gamma \\in [-\\pi,\\pi]
898 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
899 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
900 :param north_m: Northing distances from origin in [m].
901 :param east_m: Easting distances from origin in [m].
902 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
903 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
904 :type lat0: float
905 :type lon0: float
907 :return: Array with latitudes and longitudes in [deg].
908 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
909 '''
911 b = math.pi/2.-lat0*d2r
912 a = num.sqrt(north_m**2+east_m**2)/earthradius
914 gamma = num.arctan2(east_m, north_m)
915 alphasign = 1.
916 alphasign = num.where(gamma < 0., -1., 1.)
917 gamma = num.abs(gamma)
919 z1 = num.cos((a-b)/2.)*num.cos(gamma/2.)
920 n1 = num.cos((a+b)/2.)*num.sin(gamma/2.)
921 z2 = num.sin((a-b)/2.)*num.cos(gamma/2.)
922 n2 = num.sin((a+b)/2.)*num.sin(gamma/2.)
923 t1 = num.arctan2(z1, n1)
924 t2 = num.arctan2(z2, n2)
926 alpha = t1 + t2
928 sin_t1 = num.sin(t1)
929 sin_t2 = num.sin(t2)
930 c = num.where(
931 num.abs(sin_t1) > num.abs(sin_t2),
932 num.arccos(z1/sin_t1)*2.,
933 num.arcsin(z2/sin_t2)*2.)
935 lat = r2d * (math.pi/2. - c)
936 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
937 return lat, lon
940def angle_difference(angle_a, angle_b):
941 '''Difference between two angles in [deg].
943 :param angle_a: Angle A [deg].
944 :param angle_b: Angle B [deg].
946 :return: Difference between the two angles in [deg].
947 :rtype: float
948 '''
949 return ((angle_a - angle_b) + 180.) % 360. + 180.
952def latlon_to_ne(*args):
953 '''
954 Relative cartesian coordinates with respect to a reference location.
956 For two locations, a reference location A and another location B, given in
957 geographical coordinates in degrees, the corresponding cartesian
958 coordinates are calculated.
959 Assisting functions are :py:func:`pyrocko.orthodrome.azimuth` and
960 :py:func:`pyrocko.orthodrome.distance_accurate50m`.
962 .. math::
964 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
965 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
966 \\mathrm{)}\\\\[0.3cm]
968 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180}
969 \\varphi_{\\mathrm{azi},AB} )\\\\
970 e &= D_{AB} \\cdot
971 \\sin( \\frac{\\pi }{180} \\varphi_{\\mathrm{azi},AB})
973 :param refloc: Location reference point.
974 :type refloc: :py:class:`pyrocko.orthodrome.Loc`
975 :param loc: Location of interest.
976 :type loc: :py:class:`pyrocko.orthodrome.Loc`
978 :return: Northing and easting from refloc to location in [m].
979 :rtype: tuple[float, float]
981 '''
983 azi = azimuth(*args)
984 dist = distance_accurate50m(*args)
985 n, e = math.cos(azi*d2r)*dist, math.sin(azi*d2r)*dist
986 return n, e
989def latlon_to_ne_numpy(lat0, lon0, lat, lon):
990 '''
991 Relative cartesian coordinates with respect to a reference location.
993 For two locations, a reference location (``lat0``, ``lon0``) and another
994 location B, given in geographical coordinates in degrees,
995 the corresponding cartesian coordinates are calculated.
996 Assisting functions are :py:func:`azimuth`
997 and :py:func:`distance_accurate50m`.
999 :param lat0: Latitude of the reference location in [deg].
1000 :param lon0: Longitude of the reference location in [deg].
1001 :param lat: Latitude of the absolute location in [deg].
1002 :param lon: Longitude of the absolute location in [deg].
1004 :return: ``(n, e)``: relative north and east positions in [m].
1005 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
1007 Implemented formulations:
1009 .. math::
1011 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
1012 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
1013 \\mathrm{)}\\\\[0.3cm]
1015 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180} \\varphi_{
1016 \\mathrm{azi},AB} )\\\\
1017 e &= D_{AB} \\cdot \\sin( \\frac{\\pi }{180} \\varphi_{
1018 \\mathrm{azi},AB} )
1019 '''
1021 azi = azimuth_numpy(lat0, lon0, lat, lon)
1022 dist = distance_accurate50m_numpy(lat0, lon0, lat, lon)
1023 n = num.cos(azi*d2r)*dist
1024 e = num.sin(azi*d2r)*dist
1025 return n, e
1028_wgs84 = None
1031def get_wgs84():
1032 global _wgs84
1033 if _wgs84 is None:
1034 from geographiclib.geodesic import Geodesic
1035 _wgs84 = Geodesic.WGS84
1037 return _wgs84
1040def amap(n):
1042 def wrap(f):
1043 if n == 1:
1044 @wraps(f)
1045 def func(*args):
1046 it = num.nditer(args + (None,))
1047 for ops in it:
1048 ops[-1][...] = f(*ops[:-1])
1050 return it.operands[-1]
1051 elif n == 2:
1052 @wraps(f)
1053 def func(*args):
1054 it = num.nditer(args + (None, None))
1055 for ops in it:
1056 ops[-2][...], ops[-1][...] = f(*ops[:-2])
1058 return it.operands[-2], it.operands[-1]
1059 else:
1060 raise ValueError("Cannot wrap %s" % f.__qualname__)
1062 return func
1063 return wrap
1066@amap(2)
1067def ne_to_latlon2(lat0, lon0, north_m, east_m):
1068 wgs84 = get_wgs84()
1069 az = num.arctan2(east_m, north_m)*r2d
1070 dist = num.sqrt(east_m**2 + north_m**2)
1071 x = wgs84.Direct(lat0, lon0, az, dist)
1072 return x['lat2'], x['lon2']
1075@amap(2)
1076def latlon_to_ne2(lat0, lon0, lat1, lon1):
1077 wgs84 = get_wgs84()
1078 x = wgs84.Inverse(lat0, lon0, lat1, lon1)
1079 dist = x['s12']
1080 az = x['azi1']
1081 n = num.cos(az*d2r)*dist
1082 e = num.sin(az*d2r)*dist
1083 return n, e
1086@amap(1)
1087def distance_accurate15nm(lat1, lon1, lat2, lon2):
1088 wgs84 = get_wgs84()
1089 return wgs84.Inverse(lat1, lon1, lat2, lon2)['s12']
1092def positive_region(region):
1093 '''
1094 Normalize parameterization of a rectangular geographical region.
1096 :param region: ``(west, east, south, north)`` in [deg].
1097 :type region: tuple of float
1099 :returns: ``(west, east, south, north)``, where ``west <= east`` and
1100 where ``west`` and ``east`` are in the range ``[-180., 180.+360.]``.
1101 :rtype: tuple of float
1102 '''
1103 west, east, south, north = [float(x) for x in region]
1105 assert -180. - 360. <= west < 180.
1106 assert -180. < east <= 180. + 360.
1107 assert -90. <= south < 90.
1108 assert -90. < north <= 90.
1110 if east < west:
1111 east += 360.
1113 if west < -180.:
1114 west += 360.
1115 east += 360.
1117 return (west, east, south, north)
1120def points_in_region(p, region):
1121 '''
1122 Check what points are contained in a rectangular geographical region.
1124 :param p: ``(lat, lon)`` pairs in [deg].
1125 :type p: :py:class:`numpy.ndarray` ``(N, 2)``
1126 :param region: ``(west, east, south, north)`` region boundaries in [deg].
1127 :type region: tuple of float
1129 :returns: Mask, returning ``True`` for each point within the region.
1130 :rtype: :py:class:`numpy.ndarray` of bool, shape ``(N)``
1131 '''
1133 w, e, s, n = positive_region(region)
1134 return num.logical_and(
1135 num.logical_and(s <= p[:, 0], p[:, 0] <= n),
1136 num.logical_or(
1137 num.logical_and(w <= p[:, 1], p[:, 1] <= e),
1138 num.logical_and(w-360. <= p[:, 1], p[:, 1] <= e-360.)))
1141def point_in_region(p, region):
1142 '''
1143 Check if a point is contained in a rectangular geographical region.
1145 :param p: ``(lat, lon)`` in [deg].
1146 :type p: tuple of float
1147 :param region: ``(west, east, south, north)`` region boundaries in [deg].
1148 :type region: tuple of float
1150 :returns: ``True``, if point is in region, else ``False``.
1151 :rtype: bool
1152 '''
1154 w, e, s, n = positive_region(region)
1155 return num.logical_and(
1156 num.logical_and(s <= p[0], p[0] <= n),
1157 num.logical_or(
1158 num.logical_and(w <= p[1], p[1] <= e),
1159 num.logical_and(w-360. <= p[1], p[1] <= e-360.)))
1162def radius_to_region(lat, lon, radius):
1163 '''
1164 Get a rectangular region which fully contains a given circular region.
1166 :param lat: Latitude of the center point of circular region in [deg].
1167 :type lat: float
1168 :param lon: Longitude of the center point of circular region in [deg].
1169 :type lon: float
1170 :param radius: Radius of circular region in [m].
1171 :type radius: float
1173 :returns: Rectangular region as ``(east, west, south, north)`` in [deg] or
1174 ``None``.
1175 :rtype: tuple[float, float, float, float]
1176 '''
1177 radius_deg = radius * m2d
1178 if radius_deg < 45.:
1179 lat_min = max(-90., lat - radius_deg)
1180 lat_max = min(90., lat + radius_deg)
1181 absmaxlat = max(abs(lat_min), abs(lat_max))
1182 if absmaxlat > 89:
1183 lon_min = -180.
1184 lon_max = 180.
1185 else:
1186 lon_min = max(
1187 -180. - 360.,
1188 lon - radius_deg / math.cos(absmaxlat*d2r))
1189 lon_max = min(
1190 180. + 360.,
1191 lon + radius_deg / math.cos(absmaxlat*d2r))
1193 lon_min, lon_max, lat_min, lat_max = positive_region(
1194 (lon_min, lon_max, lat_min, lat_max))
1196 return lon_min, lon_max, lat_min, lat_max
1198 else:
1199 return None
1202def geographic_midpoint(lats, lons, weights=None):
1203 '''
1204 Calculate geographic midpoints by finding the center of gravity.
1206 This method suffers from instabilities if points are centered around the
1207 poles.
1209 :param lats: Latitudes in [deg].
1210 :param lons: Longitudes in [deg].
1211 :param weights: Weighting factors.
1212 :type lats: :py:class:`numpy.ndarray`, ``(N)``
1213 :type lons: :py:class:`numpy.ndarray`, ``(N)``
1214 :type weights: optional, :py:class:`numpy.ndarray`, ``(N)``
1216 :return: Latitudes and longitudes of the midpoints in [deg].
1217 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
1218 '''
1219 if not weights:
1220 weights = num.ones(len(lats))
1222 total_weigth = num.sum(weights)
1223 weights /= total_weigth
1224 lats = lats * d2r
1225 lons = lons * d2r
1226 x = num.sum(num.cos(lats) * num.cos(lons) * weights)
1227 y = num.sum(num.cos(lats) * num.sin(lons) * weights)
1228 z = num.sum(num.sin(lats) * weights)
1230 lon = num.arctan2(y, x)
1231 hyp = num.sqrt(x**2 + y**2)
1232 lat = num.arctan2(z, hyp)
1234 return lat/d2r, lon/d2r
1237def geographic_midpoint_locations(locations, weights=None):
1238 coords = num.array([loc.effective_latlon
1239 for loc in locations])
1240 return geographic_midpoint(coords[:, 0], coords[:, 1], weights)
1243def geodetic_to_ecef(lat, lon, alt):
1244 '''
1245 Convert geodetic coordinates to Earth-Centered, Earth-Fixed (ECEF)
1246 Cartesian coordinates. [#1]_ [#2]_
1248 :param lat: Geodetic latitude in [deg].
1249 :param lon: Geodetic longitude in [deg].
1250 :param alt: Geodetic altitude (height) in [m] (positive for points outside
1251 the geoid).
1252 :type lat: float
1253 :type lon: float
1254 :type alt: float
1256 :return: ECEF Cartesian coordinates (X, Y, Z) in [m].
1257 :rtype: tuple[float, float, float]
1259 .. [#1] https://en.wikipedia.org/wiki/ECEF
1260 .. [#2] https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1261 #From_geodetic_to_ECEF_coordinates
1262 '''
1264 f = earth_oblateness
1265 a = earthradius_equator
1266 e2 = 2*f - f**2
1268 lat, lon = num.radians(lat), num.radians(lon)
1269 # Normal (plumb line)
1270 N = a / num.sqrt(1.0 - (e2 * num.sin(lat)**2))
1272 X = (N+alt) * num.cos(lat) * num.cos(lon)
1273 Y = (N+alt) * num.cos(lat) * num.sin(lon)
1274 Z = (N*(1.0-e2) + alt) * num.sin(lat)
1276 return (X, Y, Z)
1279def ecef_to_geodetic(X, Y, Z):
1280 '''
1281 Convert Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates to
1282 geodetic coordinates (Ferrari's solution).
1284 :param X, Y, Z: Cartesian coordinates in ECEF system in [m].
1285 :type X, Y, Z: float
1287 :return: Geodetic coordinates (lat, lon, alt). Latitude and longitude are
1288 in [deg] and altitude is in [m]
1289 (positive for points outside the geoid).
1290 :rtype: tuple, float
1292 .. seealso ::
1293 https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1294 #The_application_of_Ferrari.27s_solution
1295 '''
1296 f = earth_oblateness
1297 a = earthradius_equator
1298 b = a * (1. - f)
1299 e2 = 2.*f - f**2
1301 # usefull
1302 a2 = a**2
1303 b2 = b**2
1304 e4 = e2**2
1305 X2 = X**2
1306 Y2 = Y**2
1307 Z2 = Z**2
1309 r = num.sqrt(X2 + Y2)
1310 r2 = r**2
1312 e_prime2 = (a2 - b2)/b2
1313 E2 = a2 - b2
1314 F = 54. * b2 * Z2
1315 G = r2 + (1.-e2)*Z2 - (e2*E2)
1316 C = (e4 * F * r2) / (G**3)
1317 S = cbrt(1. + C + num.sqrt(C**2 + 2.*C))
1318 P = F / (3. * (S + 1./S + 1.)**2 * G**2)
1319 Q = num.sqrt(1. + (2.*e4*P))
1321 dum1 = -(P*e2*r) / (1.+Q)
1322 dum2 = 0.5 * a2 * (1. + 1./Q)
1323 dum3 = (P * (1.-e2) * Z2) / (Q * (1.+Q))
1324 dum4 = 0.5 * P * r2
1325 r0 = dum1 + num.sqrt(dum2 - dum3 - dum4)
1327 U = num.sqrt((r - e2*r0)**2 + Z2)
1328 V = num.sqrt((r - e2*r0)**2 + (1.-e2)*Z2)
1329 Z0 = (b2*Z) / (a*V)
1331 alt = U * (1. - (b2 / (a*V)))
1332 lat = num.arctan((Z + e_prime2 * Z0)/r)
1333 lon = num.arctan2(Y, X)
1335 return (lat*r2d, lon*r2d, alt)
1338class Farside(Exception):
1339 pass
1342def latlon_to_xyz(latlons):
1343 if latlons.ndim == 1:
1344 return latlon_to_xyz(latlons[num.newaxis, :])[0]
1346 points = num.zeros((latlons.shape[0], 3))
1347 lats = latlons[:, 0]
1348 lons = latlons[:, 1]
1349 points[:, 0] = num.cos(lats*d2r) * num.cos(lons*d2r)
1350 points[:, 1] = num.cos(lats*d2r) * num.sin(lons*d2r)
1351 points[:, 2] = num.sin(lats*d2r)
1352 return points
1355def xyz_to_latlon(xyz):
1356 if xyz.ndim == 1:
1357 return xyz_to_latlon(xyz[num.newaxis, :])[0]
1359 latlons = num.zeros((xyz.shape[0], 2))
1360 latlons[:, 0] = num.arctan2(
1361 xyz[:, 2], num.sqrt(xyz[:, 0]**2 + xyz[:, 1]**2)) * r2d
1362 latlons[:, 1] = num.arctan2(
1363 xyz[:, 1], xyz[:, 0]) * r2d
1364 return latlons
1367def rot_to_00(lat, lon):
1368 rot0 = euler_to_matrix(0., -90.*d2r, 0.).A
1369 rot1 = euler_to_matrix(-d2r*lat, 0., -d2r*lon).A
1370 return num.dot(rot0.T, num.dot(rot1, rot0)).T
1373def distances3d(a, b):
1374 return num.sqrt(num.sum((a-b)**2, axis=a.ndim-1))
1377def circulation(points2):
1378 return num.sum(
1379 (points2[1:, 0] - points2[:-1, 0])
1380 * (points2[1:, 1] + points2[:-1, 1]))
1383def stereographic(points):
1384 dists = distances3d(points[1:, :], points[:-1, :])
1385 if dists.size > 0:
1386 maxdist = num.max(dists)
1387 cutoff = maxdist**2 / 2.
1388 else:
1389 cutoff = 1.0e-5
1391 points = points.copy()
1392 if num.any(points[:, 0] < -1. + cutoff):
1393 raise Farside()
1395 points_out = points[:, 1:].copy()
1396 factor = 1.0 / (1.0 + points[:, 0])
1397 points_out *= factor[:, num.newaxis]
1399 return points_out
1402def stereographic_poly(points):
1403 dists = distances3d(points[1:, :], points[:-1, :])
1404 if dists.size > 0:
1405 maxdist = num.max(dists)
1406 cutoff = maxdist**2 / 2.
1407 else:
1408 cutoff = 1.0e-5
1410 points = points.copy()
1411 if num.any(points[:, 0] < -1. + cutoff):
1412 raise Farside()
1414 points_out = points[:, 1:].copy()
1415 factor = 1.0 / (1.0 + points[:, 0])
1416 points_out *= factor[:, num.newaxis]
1418 if circulation(points_out) >= 0:
1419 raise Farside()
1421 return points_out
1424def gnomonic_x(points, cutoff=0.01):
1425 points_out = points[:, 1:].copy()
1426 if num.any(points[:, 0] < cutoff):
1427 raise Farside()
1429 factor = 1.0 / points[:, 0]
1430 points_out *= factor[:, num.newaxis]
1431 return points_out
1434def cneg(i, x):
1435 if i == 1:
1436 return x
1437 else:
1438 return num.logical_not(x)
1441def contains_points(polygon, points):
1442 '''
1443 Test which points are inside polygon on a sphere.
1445 The inside of the polygon is defined as the area which is to the left hand
1446 side of an observer walking the polygon line, points in order, on the
1447 sphere. Lines between the polygon points are treated as great circle paths.
1448 The polygon may be arbitrarily complex, as long as it does not have any
1449 crossings or thin parts with zero width. The polygon may contain the poles
1450 and is allowed to wrap around the sphere multiple times.
1452 The algorithm works by consecutive cutting of the polygon into (almost)
1453 hemispheres and subsequent Gnomonic projections to perform the
1454 point-in-polygon tests on a 2D plane.
1456 :param polygon: Point coordinates defining the polygon [deg].
1457 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1458 0=lat, 1=lon
1459 :param points: Coordinates of points to test [deg].
1460 :type points: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1461 0=lat, 1=lon
1463 :returns: Boolean mask array.
1464 :rtype: :py:class:`numpy.ndarray` of shape ``(N,)``.
1465 '''
1467 and_ = num.logical_and
1469 points_xyz = latlon_to_xyz(points)
1470 mask_x = 0. <= points_xyz[:, 0]
1471 mask_y = 0. <= points_xyz[:, 1]
1472 mask_z = 0. <= points_xyz[:, 2]
1474 result = num.zeros(points.shape[0], dtype=int)
1476 for ix in [-1, 1]:
1477 for iy in [-1, 1]:
1478 for iz in [-1, 1]:
1479 mask = and_(
1480 and_(cneg(ix, mask_x), cneg(iy, mask_y)),
1481 cneg(iz, mask_z))
1483 center_xyz = num.array([ix, iy, iz], dtype=float)
1485 lat, lon = xyz_to_latlon(center_xyz)
1486 rot = rot_to_00(lat, lon)
1488 points_rot_xyz = num.dot(rot, points_xyz[mask, :].T).T
1489 points_rot_pro = gnomonic_x(points_rot_xyz)
1491 offset = 0.01
1493 poly_xyz = latlon_to_xyz(polygon)
1494 poly_rot_xyz = num.dot(rot, poly_xyz.T).T
1495 poly_rot_xyz[:, 0] -= offset
1496 groups = spoly_cut([poly_rot_xyz], axis=0)
1498 for poly_rot_group_xyz in groups[1]:
1499 poly_rot_group_xyz[:, 0] += offset
1501 poly_rot_group_pro = gnomonic_x(
1502 poly_rot_group_xyz)
1504 if circulation(poly_rot_group_pro) > 0:
1505 result[mask] += path_contains_points(
1506 poly_rot_group_pro, points_rot_pro)
1507 else:
1508 result[mask] -= path_contains_points(
1509 poly_rot_group_pro, points_rot_pro)
1511 return result.astype(num.bool)
1514def contains_point(polygon, point):
1515 '''
1516 Test if point is inside polygon on a sphere.
1518 Convenience wrapper to :py:func:`contains_points` to test a single point.
1520 :param polygon: Point coordinates defining the polygon [deg].
1521 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1522 0=lat, 1=lon
1523 :param point: Coordinates ``(lat, lon)`` of point to test [deg].
1524 :type point: tuple of float
1526 :returns: ``True``, if point is located within polygon, else ``False``.
1527 :rtype: bool
1528 '''
1530 return bool(
1531 contains_points(polygon, num.asarray(point)[num.newaxis, :])[0])