Coverage for /usr/local/lib/python3.13/dist-packages/pyrocko/orthodrome.py: 77%
464 statements
« prev ^ index » next coverage.py v7.6.0, created at 2025-12-04 10:41 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7Some basic geodetic functions.
8'''
10import math
11from functools import wraps
13import numpy as num
14from matplotlib.path import Path
16from .config import config
17from .moment_tensor import euler_to_matrix
18from .plot.beachball import spoly_cut
20d2r = math.pi/180.
21r2d = 1./d2r
22earth_oblateness = 1./298.257223563
23earthradius_equator = 6378.14 * 1000.
24earthradius = config().earthradius
25d2m = earthradius_equator*math.pi/180.
26m2d = 1./d2m
28_testpath = Path([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)], closed=True)
30raise_if_slow_path_contains_points = False
33class Slow(Exception):
34 pass
37if hasattr(_testpath, 'contains_points') and num.all(
38 _testpath.contains_points([(0.5, 0.5), (1.5, 0.5)]) == [True, False]):
40 def path_contains_points(verts, points):
41 p = Path(verts, closed=True)
42 return p.contains_points(points).astype(bool)
44else:
45 # work around missing contains_points and bug in matplotlib ~ v1.2.0
47 def path_contains_points(verts, points):
48 if raise_if_slow_path_contains_points:
49 # used by unit test to skip slow gshhg_example.py
50 raise Slow()
52 p = Path(verts, closed=True)
53 result = num.zeros(points.shape[0], dtype=bool)
54 for i in range(result.size):
55 result[i] = p.contains_point(points[i, :])
57 return result
60try:
61 cbrt = num.cbrt
62except AttributeError:
63 def cbrt(x):
64 return x**(1./3.)
67def float_array_broadcast(*args):
68 return num.broadcast_arrays(*[
69 num.asarray(x, dtype=float) for x in args])
72class Loc(object):
73 '''
74 Simple location representation.
76 :attrib lat: Latitude in [deg].
77 :attrib lon: Longitude in [deg].
78 '''
79 __slots__ = ['lat', 'lon']
81 def __init__(self, lat, lon):
82 self.lat = lat
83 self.lon = lon
86def clip(x, mi, ma):
87 '''
88 Clip values of an array.
90 :param x: Continunous data to be clipped.
91 :param mi: Clip minimum.
92 :param ma: Clip maximum.
93 :type x: :py:class:`numpy.ndarray`
94 :type mi: float
95 :type ma: float
97 :return: Clipped data.
98 :rtype: :py:class:`numpy.ndarray`
99 '''
100 return num.minimum(num.maximum(mi, x), ma)
103def wrap(x, mi, ma):
104 '''
105 Wrapping continuous data to fundamental phase values.
107 .. math::
108 x_{\\mathrm{wrapped}} = x_{\\mathrm{cont},i} -
109 \\frac{ x_{\\mathrm{cont},i} - r_{\\mathrm{min}} }
110 { r_{\\mathrm{max}} - r_{\\mathrm{min}}}
111 \\cdot ( r_{\\mathrm{max}} - r_{\\mathrm{min}}),\\quad
112 x_{\\mathrm{wrapped}}\\; \\in
113 \\;[ r_{\\mathrm{min}},\\, r_{\\mathrm{max}}].
115 :param x: Continunous data to be wrapped.
116 :param mi: Minimum value of wrapped data.
117 :param ma: Maximum value of wrapped data.
118 :type x: :py:class:`numpy.ndarray`
119 :type mi: float
120 :type ma: float
122 :return: Wrapped data.
123 :rtype: :py:class:`numpy.ndarray`
124 '''
125 return x - num.floor((x-mi)/(ma-mi)) * (ma-mi)
128def _latlon_pair(args):
129 if len(args) == 2:
130 a, b = args
131 return a.lat, a.lon, b.lat, b.lon
133 elif len(args) == 4:
134 return args
137def cosdelta(*args):
138 '''
139 Cosine of the angular distance between two points ``a`` and ``b`` on a
140 sphere.
142 This function (find implementation below) returns the cosine of the
143 distance angle 'delta' between two points ``a`` and ``b``, coordinates of
144 which are expected to be given in geographical coordinates and in degrees.
145 For numerical stability a maximum of 1.0 is enforced.
147 .. math::
149 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
150 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
151 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
152 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\[0.5cm]
154 \\cos(\\Delta) = \\min( 1.0, \\quad \\sin( A_{\\mathrm{lat'}})
155 \\sin( B_{\\mathrm{lat'}} ) +
156 \\cos(A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
157 \\cos( B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )
159 :param a: Location point A.
160 :type a: :py:class:`pyrocko.orthodrome.Loc`
161 :param b: Location point B.
162 :type b: :py:class:`pyrocko.orthodrome.Loc`
164 :return: Cosdelta.
165 :rtype: float
166 '''
168 alat, alon, blat, blon = _latlon_pair(args)
170 return min(
171 1.0,
172 math.sin(alat*d2r) * math.sin(blat*d2r) +
173 math.cos(alat*d2r) * math.cos(blat*d2r) *
174 math.cos(d2r*(blon-alon)))
177def cosdelta_numpy(a_lats, a_lons, b_lats, b_lons):
178 '''
179 Cosine of the angular distance between two points ``a`` and ``b`` on a
180 sphere.
182 This function returns the cosines of the distance
183 angles *delta* between two points ``a`` and ``b`` given as
184 :py:class:`numpy.ndarray`.
185 The coordinates are expected to be given in geographical coordinates
186 and in degrees. For numerical stability a maximum of ``1.0`` is enforced.
188 Please find the details of the implementation in the documentation of
189 the function :py:func:`pyrocko.orthodrome.cosdelta` above.
191 :param a_lats: Latitudes in [deg] point A.
192 :param a_lons: Longitudes in [deg] point A.
193 :param b_lats: Latitudes in [deg] point B.
194 :param b_lons: Longitudes in [deg] point B.
195 :type a_lats: :py:class:`numpy.ndarray`
196 :type a_lons: :py:class:`numpy.ndarray`
197 :type b_lats: :py:class:`numpy.ndarray`
198 :type b_lons: :py:class:`numpy.ndarray`
200 :return: Cosdelta.
201 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
202 '''
203 return num.minimum(
204 1.0,
205 num.sin(a_lats*d2r) * num.sin(b_lats*d2r) +
206 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
207 num.cos(d2r*(b_lons-a_lons)))
210def azimuth(*args):
211 '''
212 Azimuth calculation.
214 This function (find implementation below) returns azimuth ...
215 between points ``a`` and ``b``, coordinates of
216 which are expected to be given in geographical coordinates and in degrees.
218 .. math::
220 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
221 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
222 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
223 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\
225 \\varphi_{\\mathrm{azi},AB} = \\frac{180}{\\pi} \\arctan \\left[
226 \\frac{
227 \\cos( A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
228 \\sin(B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )}
229 {\\sin ( B_{\\mathrm{lat'}} ) - \\sin( A_{\\mathrm{lat'}}
230 cosdelta) } \\right]
232 :param a: Location point A.
233 :type a: :py:class:`pyrocko.orthodrome.Loc`
234 :param b: Location point B.
235 :type b: :py:class:`pyrocko.orthodrome.Loc`
237 :return: Azimuth in degree
238 '''
240 alat, alon, blat, blon = _latlon_pair(args)
242 return r2d*math.atan2(
243 math.cos(alat*d2r) * math.cos(blat*d2r) *
244 math.sin(d2r*(blon-alon)),
245 math.sin(d2r*blat) - math.sin(d2r*alat) * cosdelta(
246 alat, alon, blat, blon))
249def azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta=None):
250 '''
251 Calculation of the azimuth (*track angle*) from a location A towards B.
253 This function returns azimuths (*track angles*) from locations A towards B
254 given in :py:class:`numpy.ndarray`. Coordinates are expected to be given in
255 geographical coordinates and in degrees.
257 Please find the details of the implementation in the documentation of the
258 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)
281def azibazi(*args, **kwargs):
282 '''
283 Azimuth and backazimuth from location A towards B and back.
285 :returns: Azimuth in [deg] from A to B, back azimuth in [deg] from B to A.
286 :rtype: tuple[float, float]
287 '''
289 alat, alon, blat, blon = _latlon_pair(args)
290 if alat == blat and alon == blon:
291 return 0., 180.
293 implementation = kwargs.get('implementation', 'c')
294 assert implementation in ('c', 'python')
295 if implementation == 'c':
296 from pyrocko import orthodrome_ext
297 return orthodrome_ext.azibazi(alat, alon, blat, blon)
299 cd = cosdelta(alat, alon, blat, blon)
300 azi = r2d*math.atan2(
301 math.cos(alat*d2r) * math.cos(blat*d2r) *
302 math.sin(d2r*(blon-alon)),
303 math.sin(d2r*blat) - math.sin(d2r*alat) * cd)
304 bazi = r2d*math.atan2(
305 math.cos(blat*d2r) * math.cos(alat*d2r) *
306 math.sin(d2r*(alon-blon)),
307 math.sin(d2r*alat) - math.sin(d2r*blat) * cd)
309 return azi, bazi
312def azibazi_numpy(a_lats, a_lons, b_lats, b_lons, implementation='c'):
313 '''
314 Azimuth and backazimuth from location A towards B and back.
316 Arguments are given as :py:class:`numpy.ndarray`.
318 :param a_lats: Latitude(s) in [deg] of point A.
319 :type a_lats: :py:class:`numpy.ndarray`
320 :param a_lons: Longitude(s) in [deg] of point A.
321 :type a_lons: :py:class:`numpy.ndarray`
322 :param b_lats: Latitude(s) in [deg] of point B.
323 :type b_lats: :py:class:`numpy.ndarray`
324 :param b_lons: Longitude(s) in [deg] of point B.
325 :type b_lons: :py:class:`numpy.ndarray`
327 :returns: Azimuth(s) in [deg] from A to B,
328 back azimuth(s) in [deg] from B to A.
329 :rtype: :py:class:`numpy.ndarray`, :py:class:`numpy.ndarray`
330 '''
332 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
333 a_lats, a_lons, b_lats, b_lons)
335 assert implementation in ('c', 'python')
336 if implementation == 'c':
337 from pyrocko import orthodrome_ext
338 return orthodrome_ext.azibazi_numpy(a_lats, a_lons, b_lats, b_lons)
340 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
341 azis = azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta)
342 bazis = azimuth_numpy(b_lats, b_lons, a_lats, a_lons, _cosdelta)
344 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
345 ii_eq = num.where(eq)[0]
346 azis[ii_eq] = 0.0
347 bazis[ii_eq] = 180.0
348 return azis, bazis
351def azidist_numpy(*args):
352 '''
353 Calculation of the azimuth (*track angle*) and the distance from locations
354 A towards B on a sphere.
356 The assisting functions used are :py:func:`pyrocko.orthodrome.cosdelta` and
357 :py:func:`pyrocko.orthodrome.azimuth`
359 :param a_lats: Latitudes in [deg] point A.
360 :param a_lons: Longitudes in [deg] point A.
361 :param b_lats: Latitudes in [deg] point B.
362 :param b_lons: Longitudes in [deg] point B.
363 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
364 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
365 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
366 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
368 :return: Azimuths in [deg], distances in [deg].
369 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
370 '''
371 _cosdelta = cosdelta_numpy(*args)
372 _azimuths = azimuth_numpy(_cosdelta=_cosdelta, *args)
373 return _azimuths, r2d*num.arccos(_cosdelta)
376def distance_accurate50m(*args, **kwargs):
377 '''
378 Accurate distance calculation based on a spheroid of rotation.
380 Function returns distance in meter between points A and B, coordinates of
381 which must be given in geographical coordinates and in degrees.
382 The returned distance should be accurate to 50 m using WGS84.
383 Values for the Earth's equator radius and the Earth's oblateness
384 (``f_oblate``) are defined in the pyrocko configuration file
385 :py:class:`pyrocko.config`.
387 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
389 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
390 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
392 .. math::
394 F = \\frac{\\pi}{180}
395 \\frac{(A_{lat} + B_{lat})}{2}, \\quad
396 G = \\frac{\\pi}{180}
397 \\frac{(A_{lat} - B_{lat})}{2}, \\quad
398 l = \\frac{\\pi}{180}
399 \\frac{(A_{lon} - B_{lon})}{2} \\quad
400 \\\\[0.5cm]
401 S = \\sin^2(G) \\cdot \\cos^2(l) +
402 \\cos^2(F) \\cdot \\sin^2(l), \\quad \\quad
403 C = \\cos^2(G) \\cdot \\cos^2(l) +
404 \\sin^2(F) \\cdot \\sin^2(l)
406 .. math::
408 w = \\arctan \\left( \\sqrt{ \\frac{S}{C}} \\right) , \\quad
409 r = \\sqrt{\\frac{S}{C} }
411 The spherical-earth distance D between A and B, can be given with:
413 .. math::
415 D_{sphere} = 2w \\cdot R_{equator}
417 The oblateness of the Earth requires some correction with
418 correction factors h1 and h2:
420 .. math::
422 h_1 = \\frac{3r - 1}{2C}, \\quad
423 h_2 = \\frac{3r +1 }{2S}\\\\[0.5cm]
425 D = D_{\\mathrm{sphere}} \\cdot [ 1 + h_1 \\,f_{\\mathrm{oblate}}
426 \\cdot \\sin^2(F)
427 \\cos^2(G) - h_2\\, f_{\\mathrm{oblate}}
428 \\cdot \\cos^2(F) \\sin^2(G)]
431 :param a: Location point A.
432 :type a: :py:class:`pyrocko.orthodrome.Loc`
433 :param b: Location point B.
434 :type b: :py:class:`pyrocko.orthodrome.Loc`
436 :return: Distance in [m].
437 :rtype: float
438 '''
440 alat, alon, blat, blon = _latlon_pair(args)
442 implementation = kwargs.get('implementation', 'c')
443 assert implementation in ('c', 'python')
444 if implementation == 'c':
445 from pyrocko import orthodrome_ext
446 return orthodrome_ext.distance_accurate50m(alat, alon, blat, blon)
448 f = (alat + blat)*d2r / 2.
449 g = (alat - blat)*d2r / 2.
450 h = (alon - blon)*d2r / 2.
452 s = math.sin(g)**2 * math.cos(h)**2 + math.cos(f)**2 * math.sin(h)**2
453 c = math.cos(g)**2 * math.cos(h)**2 + math.sin(f)**2 * math.sin(h)**2
455 w = math.atan(math.sqrt(s/c))
457 if w == 0.0:
458 return 0.0
460 r = math.sqrt(s*c)/w
461 d = 2.*w*earthradius_equator
462 h1 = (3.*r-1.)/(2.*c)
463 h2 = (3.*r+1.)/(2.*s)
465 return d * (1. +
466 earth_oblateness * h1 * math.sin(f)**2 * math.cos(g)**2 -
467 earth_oblateness * h2 * math.cos(f)**2 * math.sin(g)**2)
470def distance_accurate50m_numpy(
471 a_lats, a_lons, b_lats, b_lons, implementation='c'):
473 '''
474 Accurate distance calculation based on a spheroid of rotation.
476 Function returns distance in meter between points ``a`` and ``b``,
477 coordinates of which must be given in geographical coordinates and in
478 degrees.
479 The returned distance should be accurate to 50 m using WGS84.
480 Values for the Earth's equator radius and the Earth's oblateness
481 (``f_oblate``) are defined in the pyrocko configuration file
482 :py:class:`pyrocko.config`.
484 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
486 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
487 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
489 .. math::
491 F_i = \\frac{\\pi}{180}
492 \\frac{(a_{lat,i} + a_{lat,i})}{2}, \\quad
493 G_i = \\frac{\\pi}{180}
494 \\frac{(a_{lat,i} - b_{lat,i})}{2}, \\quad
495 l_i= \\frac{\\pi}{180}
496 \\frac{(a_{lon,i} - b_{lon,i})}{2} \\\\[0.5cm]
497 S_i = \\sin^2(G_i) \\cdot \\cos^2(l_i) +
498 \\cos^2(F_i) \\cdot \\sin^2(l_i), \\quad \\quad
499 C_i = \\cos^2(G_i) \\cdot \\cos^2(l_i) +
500 \\sin^2(F_i) \\cdot \\sin^2(l_i)
502 .. math::
504 w_i = \\arctan \\left( \\sqrt{\\frac{S_i}{C_i}} \\right), \\quad
505 r_i = \\sqrt{\\frac{S_i}{C_i} }
507 The spherical-earth distance ``D`` between ``a`` and ``b``,
508 can be given with:
510 .. math::
512 D_{\\mathrm{sphere},i} = 2w_i \\cdot R_{\\mathrm{equator}}
514 The oblateness of the Earth requires some correction with
515 correction factors ``h1`` and ``h2``:
517 .. math::
519 h_{1.i} = \\frac{3r - 1}{2C_i}, \\quad
520 h_{2,i} = \\frac{3r +1 }{2S_i}\\\\[0.5cm]
522 D_{AB,i} = D_{\\mathrm{sphere},i} \\cdot [1 + h_{1,i}
523 \\,f_{\\mathrm{oblate}}
524 \\cdot \\sin^2(F_i)
525 \\cos^2(G_i) - h_{2,i}\\, f_{\\mathrm{oblate}}
526 \\cdot \\cos^2(F_i) \\sin^2(G_i)]
528 :param a_lats: Latitudes in [deg] point A.
529 :param a_lons: Longitudes in [deg] point A.
530 :param b_lats: Latitudes in [deg] point B.
531 :param b_lons: Longitudes in [deg] point B.
532 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
533 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
534 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
535 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
537 :return: Distances in [m].
538 :rtype: :py:class:`numpy.ndarray`, ``(N)``
539 '''
541 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
542 a_lats, a_lons, b_lats, b_lons)
544 assert implementation in ('c', 'python')
545 if implementation == 'c':
546 from pyrocko import orthodrome_ext
547 return orthodrome_ext.distance_accurate50m_numpy(
548 a_lats, a_lons, b_lats, b_lons)
550 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
551 ii_neq = num.where(num.logical_not(eq))[0]
553 if num.all(eq):
554 return num.zeros_like(eq, dtype=float)
556 def extr(x):
557 if isinstance(x, num.ndarray) and x.size > 1:
558 return x[ii_neq]
559 else:
560 return x
562 a_lats = extr(a_lats)
563 a_lons = extr(a_lons)
564 b_lats = extr(b_lats)
565 b_lons = extr(b_lons)
567 f = (a_lats + b_lats)*d2r / 2.
568 g = (a_lats - b_lats)*d2r / 2.
569 h = (a_lons - b_lons)*d2r / 2.
571 s = num.sin(g)**2 * num.cos(h)**2 + num.cos(f)**2 * num.sin(h)**2
572 c = num.cos(g)**2 * num.cos(h)**2 + num.sin(f)**2 * num.sin(h)**2
574 w = num.arctan(num.sqrt(s/c))
576 r = num.sqrt(s*c)/w
578 d = 2.*w*earthradius_equator
579 h1 = (3.*r-1.)/(2.*c)
580 h2 = (3.*r+1.)/(2.*s)
582 dists = num.zeros(eq.size, dtype=float)
583 dists[ii_neq] = d * (
584 1. +
585 earth_oblateness * h1 * num.sin(f)**2 * num.cos(g)**2 -
586 earth_oblateness * h2 * num.cos(f)**2 * num.sin(g)**2)
588 return dists
591def ne_to_latlon(lat0, lon0, north_m, east_m):
592 '''
593 Transform local cartesian coordinates to latitude and longitude.
595 From east and north coordinates (``x`` and ``y`` coordinate
596 :py:class:`numpy.ndarray`) relative to a reference differences in
597 longitude and latitude are calculated, which are effectively changes in
598 azimuth and distance, respectively:
600 .. math::
602 \\text{distance change:}\\; \\Delta {\\bf{a}} &= \\sqrt{{\\bf{y}}^2 +
603 {\\bf{x}}^2 }/ \\mathrm{R_E},
605 \\text{azimuth change:}\\; \\Delta \\bf{\\gamma} &= \\arctan( \\bf{x}
606 / \\bf{y}).
608 The projection used preserves the azimuths of the input points.
610 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
611 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
612 :param north_m: Northing distances from origin in [m].
613 :param east_m: Easting distances from origin in [m].
614 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
615 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
616 :type lat0: float
617 :type lon0: float
619 :return: Array with latitudes and longitudes in [deg].
620 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
622 '''
624 a = num.sqrt(north_m**2+east_m**2)/earthradius
625 gamma = num.arctan2(east_m, north_m)
627 return azidist_to_latlon_rad(lat0, lon0, gamma, a)
630def azidist_to_latlon(lat0, lon0, azimuth_deg, distance_deg):
631 '''
632 Absolute latitudes and longitudes are calculated from relative changes.
634 Convenience wrapper to :py:func:`azidist_to_latlon_rad` with azimuth and
635 distance given in degrees.
637 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
638 :type lat0: float
639 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
640 :type lon0: float
641 :param azimuth_deg: Azimuth from origin in [deg].
642 :type azimuth_deg: :py:class:`numpy.ndarray`, ``(N)``
643 :param distance_deg: Distances from origin in [deg].
644 :type distance_deg: :py:class:`numpy.ndarray`, ``(N)``
646 :return: Array with latitudes and longitudes in [deg].
647 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
648 '''
650 return azidist_to_latlon_rad(
651 lat0, lon0, azimuth_deg/180.*num.pi, distance_deg/180.*num.pi)
654def azidist_to_latlon_rad(lat0, lon0, azimuth_rad, distance_rad):
655 '''
656 Absolute latitudes and longitudes are calculated from relative changes.
658 For numerical stability a range between of ``-1.0`` and ``1.0`` is
659 enforced for ``c`` and ``alpha``.
661 .. math::
663 \\Delta {\\bf a}_i \\; \\text{and} \\; \\Delta \\gamma_i \\;
664 \\text{are relative distances and azimuths from lat0 and lon0 for
665 \\textit{i} source points of a finite source.}
667 .. math::
669 \\mathrm{b} &= \\frac{\\pi}{2} -\\frac{\\pi}{180}\\;\\mathrm{lat_0}\\\\
670 {\\bf c}_i &=\\arccos[\\; \\cos(\\Delta {\\bf{a}}_i)
671 \\cos(\\mathrm{b}) + |\\Delta \\gamma_i| \\,
672 \\sin(\\Delta {\\bf a}_i)
673 \\sin(\\mathrm{b})\\; ] \\\\
674 \\mathrm{lat}_i &= \\frac{180}{\\pi}
675 \\left(\\frac{\\pi}{2} - {\\bf c}_i \\right)
677 .. math::
679 \\alpha_i &= \\arcsin \\left[ \\; \\frac{ \\sin(\\Delta {\\bf a}_i )
680 \\sin(|\\Delta \\gamma_i|)}{\\sin({\\bf c}_i)}\\;
681 \\right] \\\\
682 \\alpha_i &= \\begin{cases}
683 \\alpha_i, &\\text{if} \\; \\cos(\\Delta {\\bf a}_i) -
684 \\cos(\\mathrm{b}) \\cos({\\bf{c}}_i) > 0, \\;
685 \\text{else} \\\\
686 \\pi - \\alpha_i, & \\text{if} \\; \\alpha_i > 0,\\;
687 \\text{else}\\\\
688 -\\pi - \\alpha_i, & \\text{if} \\; \\alpha_i < 0.
689 \\end{cases} \\\\
690 \\mathrm{lon}_i &= \\mathrm{lon_0} +
691 \\frac{180}{\\pi} \\,
692 \\frac{\\Delta \\gamma_i }{|\\Delta \\gamma_i|}
693 \\cdot \\alpha_i
694 \\text{, with $\\alpha_i \\in [-\\pi,\\pi]$}
696 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
697 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
698 :param distance_rad: Distances from origin in [rad].
699 :param azimuth_rad: Azimuth from origin in [rad].
700 :type distance_rad: :py:class:`numpy.ndarray`, ``(N)``
701 :type azimuth_rad: :py:class:`numpy.ndarray`, ``(N)``
702 :type lat0: float
703 :type lon0: float
705 :return: Array with latitudes and longitudes in [deg].
706 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
707 '''
709 a = distance_rad
710 gamma = azimuth_rad
712 b = math.pi/2.-lat0*d2r
714 alphasign = 1.
715 alphasign = num.where(gamma < 0, -1., 1.)
716 gamma = num.abs(gamma)
718 c = num.arccos(clip(
719 num.cos(a)*num.cos(b)+num.sin(a)*num.sin(b)*num.cos(gamma), -1., 1.))
721 alpha = num.arcsin(clip(
722 num.sin(a)*num.sin(gamma)/num.sin(c), -1., 1.))
724 alpha = num.where(
725 num.cos(a)-num.cos(b)*num.cos(c) < 0,
726 num.where(alpha > 0, math.pi-alpha, -math.pi-alpha),
727 alpha)
729 lat = r2d * (math.pi/2. - c)
730 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
732 return lat, lon
735def crosstrack_distance(begin_lat, begin_lon, end_lat, end_lon,
736 point_lat, point_lon):
737 '''Calculate distance of a point to a great-circle path.
739 The sign of the results shows side of the path the point is on. Negative
740 distance is right of the path, positive left.
742 .. math ::
744 d_{xt} = \\arcsin \\left( \\sin \\left( \\Delta_{13} \\right) \\cdot
745 \\sin \\left( \\gamma_{13} - \\gamma_{12} \\right) \\right)
747 :param begin_lat: Latitude of the great circle start in [deg].
748 :param begin_lon: Longitude of the great circle start in [deg].
749 :param end_lat: Latitude of the great circle end in [deg].
750 :param end_lon: Longitude of the great circle end in [deg].
751 :param point_lat: Latitude of the point in [deg].
752 :param point_lon: Longitude of the point in [deg].
753 :type begin_lat: float
754 :type begin_lon: float
755 :type end_lat: float
756 :type end_lon: float
757 :type point_lat: float
758 :type point_lon: float
760 :return: Distance of the point to the great-circle path in [deg].
761 :rtype: float
762 '''
763 start = Loc(begin_lat, begin_lon)
764 end = Loc(end_lat, end_lon)
765 point = Loc(point_lat, point_lon)
767 dist_ang = math.acos(cosdelta(start, point))
768 azi_point = azimuth(start, point) * d2r
769 azi_end = azimuth(start, end) * d2r
771 return math.asin(math.sin(dist_ang) * math.sin(azi_point - azi_end)) * r2d
774def alongtrack_distance(begin_lat, begin_lon, end_lat, end_lon,
775 point_lat, point_lon):
776 '''Calculate distance of a point along a great-circle path in [deg].
778 Distance is relative to the beginning of the path. Negative distance is
779 before the beginning of the path.
781 .. math ::
783 d_{At} = \\arccos \\left( \\frac{\\cos \\left( \\Delta_{13} \\right)}
784 {\\cos \\left( \\Delta_{xt} \\right) } \\right)
786 :param begin_lat: Latitude of the great circle start in [deg].
787 :param begin_lon: Longitude of the great circle start in [deg].
788 :param end_lat: Latitude of the great circle end in [deg].
789 :param end_lon: Longitude of the great circle end in [deg].
790 :param point_lat: Latitude of the point in [deg].
791 :param point_lon: Longitude of the point in [deg].
792 :type begin_lat: float
793 :type begin_lon: float
794 :type end_lat: float
795 :type end_lon: float
796 :type point_lat: float
797 :type point_lon: float
799 :return: Distance of the point along the great-circle path in [deg].
800 :rtype: float
801 '''
802 begin = Loc(begin_lat, begin_lon)
803 end = Loc(end_lat, end_lon)
804 point = Loc(point_lat, point_lon)
806 def along_distance(begin, end):
807 cos_dist_ang = cosdelta(begin, point)
808 dist_rad = crosstrack_distance(
809 begin.lat, begin.lon, end.lat, end.lon, point.lat, point.lon) * d2r
810 return math.acos(cos_dist_ang / math.cos(dist_rad)) * r2d
812 distance_profile = math.acos(cosdelta(begin, end)) * r2d
813 distance_begin = along_distance(begin, end)
814 distance_end = along_distance(end, begin)
816 if distance_end > distance_profile:
817 return -distance_begin
818 return distance_begin
821def alongtrack_distance_m(begin_lat, begin_lon, end_lat, end_lon,
822 point_lat, point_lon):
823 '''Calculate distance of a point along a great-circle path in [m].
825 Distance is relative to the beginning of the path.
827 .. math ::
829 d_{At} = \\arccos \\left( \\frac{\\cos \\left( \\Delta_{13} \\right)}
830 {\\cos \\left( \\Delta_{xt} \\right) } \\right)
832 :param begin_lat: Latitude of the great circle start in [deg].
833 :param begin_lon: Longitude of the great circle start in [deg].
834 :param end_lat: Latitude of the great circle end in [deg].
835 :param end_lon: Longitude of the great circle end in [deg].
836 :param point_lat: Latitude of the point in [deg].
837 :param point_lon: Longitude of the point in [deg].
838 :type begin_lat: float
839 :type begin_lon: float
840 :type end_lat: float
841 :type end_lon: float
842 :type point_lat: float
843 :type point_lon: float
845 :return: Distance of the point along the great-circle path in [m].
846 :rtype: float
847 '''
848 start = Loc(begin_lat, begin_lon)
849 end = Loc(end_lat, end_lon)
850 azi_end = azimuth(start, end)
851 dist_deg = alongtrack_distance(
852 begin_lat, begin_lon, end_lat, end_lon,
853 point_lat, point_lon)
854 along_point = Loc(
855 *azidist_to_latlon(begin_lat, begin_lon, azi_end, dist_deg))
857 distance = distance_accurate50m(start, along_point)
858 if dist_deg < 0:
859 distance = -distance
860 return distance
863def ne_to_latlon_alternative_method(lat0, lon0, north_m, east_m):
864 '''
865 Transform local cartesian coordinates to latitude and longitude.
867 Like :py:func:`pyrocko.orthodrome.ne_to_latlon`,
868 but this method (implementation below), although it should be numerically
869 more stable, suffers problems at points which are *across the pole*
870 as seen from the cartesian origin.
872 .. math::
874 \\text{distance change:}\\; \\Delta {{\\bf a}_i} &=
875 \\sqrt{{\\bf{y}}^2_i + {\\bf{x}}^2_i }/ \\mathrm{R_E},\\\\
876 \\text{azimuth change:}\\; \\Delta {\\bf \\gamma}_i &=
877 \\arctan( {\\bf x}_i {\\bf y}_i). \\\\
878 \\mathrm{b} &=
879 \\frac{\\pi}{2} -\\frac{\\pi}{180} \\;\\mathrm{lat_0}\\\\
881 .. math::
883 {{\\bf z}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i -
884 \\mathrm{b}}{2} \\right)}
885 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
886 {{\\bf n}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i +
887 \\mathrm{b}}{2} \\right)}
888 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
889 {{\\bf z}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i -
890 \\mathrm{b}}{2} \\right)}
891 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
892 {{\\bf n}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i +
893 \\mathrm{b}}{2} \\right)}
894 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
895 {{\\bf t}_1}_i &= \\arctan{\\left( \\frac{{{\\bf z}_1}_i}
896 {{{\\bf n}_1}_i} \\right) }\\\\
897 {{\\bf t}_2}_i &= \\arctan{\\left( \\frac{{{\\bf z}_2}_i}
898 {{{\\bf n}_2}_i} \\right) } \\\\[0.5cm]
899 c &= \\begin{cases}
900 2 \\cdot \\arccos \\left( {{\\bf z}_1}_i / \\sin({{\\bf t}_1}_i)
901 \\right),\\; \\text{if }
902 |\\sin({{\\bf t}_1}_i)| >
903 |\\sin({{\\bf t}_2}_i)|,\\; \\text{else} \\\\
904 2 \\cdot \\arcsin{\\left( {{\\bf z}_2}_i /
905 \\sin({{\\bf t}_2}_i) \\right)}.
906 \\end{cases}\\\\
908 .. math::
910 {\\bf {lat}}_i &= \\frac{180}{ \\pi } \\left( \\frac{\\pi}{2}
911 - {\\bf {c}}_i \\right) \\\\
912 {\\bf {lon}}_i &= {\\bf {lon}}_0 + \\frac{180}{ \\pi }
913 \\frac{\\gamma_i}{|\\gamma_i|},
914 \\text{ with}\\; \\gamma \\in [-\\pi,\\pi]
916 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
917 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
918 :param north_m: Northing distances from origin in [m].
919 :param east_m: Easting distances from origin in [m].
920 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
921 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
922 :type lat0: float
923 :type lon0: float
925 :return: Array with latitudes and longitudes in [deg].
926 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
927 '''
929 b = math.pi/2.-lat0*d2r
930 a = num.sqrt(north_m**2+east_m**2)/earthradius
932 gamma = num.arctan2(east_m, north_m)
933 alphasign = 1.
934 alphasign = num.where(gamma < 0., -1., 1.)
935 gamma = num.abs(gamma)
937 z1 = num.cos((a-b)/2.)*num.cos(gamma/2.)
938 n1 = num.cos((a+b)/2.)*num.sin(gamma/2.)
939 z2 = num.sin((a-b)/2.)*num.cos(gamma/2.)
940 n2 = num.sin((a+b)/2.)*num.sin(gamma/2.)
941 t1 = num.arctan2(z1, n1)
942 t2 = num.arctan2(z2, n2)
944 alpha = t1 + t2
946 sin_t1 = num.sin(t1)
947 sin_t2 = num.sin(t2)
948 c = num.where(
949 num.abs(sin_t1) > num.abs(sin_t2),
950 num.arccos(z1/sin_t1)*2.,
951 num.arcsin(z2/sin_t2)*2.)
953 lat = r2d * (math.pi/2. - c)
954 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
955 return lat, lon
958def angle_difference(angle_a, angle_b):
959 '''
960 Difference between two angles in [deg].
962 :param angle_a: Angle A [deg].
963 :param angle_b: Angle B [deg].
965 :return: Difference between the two angles in [deg].
966 :rtype: float
967 '''
968 return ((angle_a - angle_b) + 180.) % 360. - 180.
971def latlon_to_ne(*args):
972 '''
973 Relative cartesian coordinates with respect to a reference location.
975 For two locations, a reference location A and another location B, given in
976 geographical coordinates in degrees, the corresponding cartesian
977 coordinates are calculated.
978 Assisting functions are :py:func:`pyrocko.orthodrome.azimuth` and
979 :py:func:`pyrocko.orthodrome.distance_accurate50m`.
981 .. math::
983 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
984 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
985 \\mathrm{)}\\\\[0.3cm]
987 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180}
988 \\varphi_{\\mathrm{azi},AB} )\\\\
989 e &= D_{AB} \\cdot
990 \\sin( \\frac{\\pi }{180} \\varphi_{\\mathrm{azi},AB})
992 :param refloc: Location reference point.
993 :type refloc: :py:class:`pyrocko.orthodrome.Loc`
994 :param loc: Location of interest.
995 :type loc: :py:class:`pyrocko.orthodrome.Loc`
997 :return: Northing and easting from refloc to location in [m].
998 :rtype: tuple[float, float]
1000 '''
1002 azi = azimuth(*args)
1003 dist = distance_accurate50m(*args)
1004 n, e = math.cos(azi*d2r)*dist, math.sin(azi*d2r)*dist
1005 return n, e
1008def latlon_to_ne_numpy(lat0, lon0, lat, lon):
1009 '''
1010 Relative cartesian coordinates with respect to a reference location.
1012 For two locations, a reference location (``lat0``, ``lon0``) and another
1013 location B, given in geographical coordinates in degrees,
1014 the corresponding cartesian coordinates are calculated.
1015 Assisting functions are :py:func:`azimuth`
1016 and :py:func:`distance_accurate50m`.
1018 :param lat0: Latitude of the reference location in [deg].
1019 :param lon0: Longitude of the reference location in [deg].
1020 :param lat: Latitude of the absolute location in [deg].
1021 :param lon: Longitude of the absolute location in [deg].
1023 :return: ``(n, e)``: relative north and east positions in [m].
1024 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
1026 Implemented formulations:
1028 .. math::
1030 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
1031 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
1032 \\mathrm{)}\\\\[0.3cm]
1034 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180} \\varphi_{
1035 \\mathrm{azi},AB} )\\\\
1036 e &= D_{AB} \\cdot \\sin( \\frac{\\pi }{180} \\varphi_{
1037 \\mathrm{azi},AB} )
1038 '''
1040 azi = azimuth_numpy(lat0, lon0, lat, lon)
1041 dist = distance_accurate50m_numpy(lat0, lon0, lat, lon)
1042 n = num.cos(azi*d2r)*dist
1043 e = num.sin(azi*d2r)*dist
1044 return n, e
1047_wgs84 = None
1050def get_wgs84():
1051 global _wgs84
1052 if _wgs84 is None:
1053 from geographiclib.geodesic import Geodesic
1054 _wgs84 = Geodesic.WGS84
1056 return _wgs84
1059def amap(n):
1061 def wrap(f):
1062 if n == 1:
1063 @wraps(f)
1064 def func(*args):
1065 it = num.nditer(args + (None,))
1066 for ops in it:
1067 ops[-1][...] = f(*ops[:-1])
1069 return it.operands[-1]
1070 elif n == 2:
1071 @wraps(f)
1072 def func(*args):
1073 it = num.nditer(args + (None, None))
1074 for ops in it:
1075 ops[-2][...], ops[-1][...] = f(*ops[:-2])
1077 return it.operands[-2], it.operands[-1]
1078 else:
1079 raise ValueError('Cannot wrap %s' % f.__qualname__)
1081 return func
1082 return wrap
1085@amap(2)
1086def ne_to_latlon2(lat0, lon0, north_m, east_m):
1087 wgs84 = get_wgs84()
1088 az = num.arctan2(east_m, north_m)*r2d
1089 dist = num.sqrt(east_m**2 + north_m**2)
1090 x = wgs84.Direct(lat0, lon0, az, dist)
1091 return x['lat2'], x['lon2']
1094@amap(2)
1095def latlon_to_ne2(lat0, lon0, lat1, lon1):
1096 wgs84 = get_wgs84()
1097 x = wgs84.Inverse(lat0, lon0, lat1, lon1)
1098 dist = x['s12']
1099 az = x['azi1']
1100 n = num.cos(az*d2r)*dist
1101 e = num.sin(az*d2r)*dist
1102 return n, e
1105@amap(1)
1106def distance_accurate15nm(lat1, lon1, lat2, lon2):
1107 wgs84 = get_wgs84()
1108 return wgs84.Inverse(lat1, lon1, lat2, lon2)['s12']
1111def positive_region(region):
1112 '''
1113 Normalize parameterization of a rectangular geographical region.
1115 :param region: ``(west, east, south, north)`` in [deg].
1116 :type region: :py:class:`tuple` of :py:class:`float`
1118 :returns: ``(west, east, south, north)``, where ``west <= east`` and
1119 where ``west`` and ``east`` are in the range ``[-180., 180.+360.]``.
1120 :rtype: :py:class:`tuple` of :py:class:`float`
1121 '''
1122 west, east, south, north = [float(x) for x in region]
1124 assert -180. - 360. <= west < 180.
1125 assert -180. < east <= 180. + 360.
1126 assert -90. <= south < 90.
1127 assert -90. < north <= 90.
1129 if east < west:
1130 east += 360.
1132 if west < -180.:
1133 west += 360.
1134 east += 360.
1136 return (west, east, south, north)
1139def points_in_region(p, region):
1140 '''
1141 Check what points are contained in a rectangular geographical region.
1143 :param p: ``(lat, lon)`` pairs in [deg].
1144 :type p: :py:class:`numpy.ndarray` ``(N, 2)``
1145 :param region: ``(west, east, south, north)`` region boundaries in [deg].
1146 :type region: :py:class:`tuple` of :py:class:`float`
1148 :returns: Mask, returning ``True`` for each point within the region.
1149 :rtype: :py:class:`numpy.ndarray` of bool, shape ``(N)``
1150 '''
1152 w, e, s, n = positive_region(region)
1153 return num.logical_and(
1154 num.logical_and(s <= p[:, 0], p[:, 0] <= n),
1155 num.logical_or(
1156 num.logical_and(w <= p[:, 1], p[:, 1] <= e),
1157 num.logical_and(w-360. <= p[:, 1], p[:, 1] <= e-360.)))
1160def point_in_region(p, region):
1161 '''
1162 Check if a point is contained in a rectangular geographical region.
1164 :param p: ``(lat, lon)`` in [deg].
1165 :type p: :py:class:`tuple` of :py:class:`float`
1166 :param region: ``(west, east, south, north)`` region boundaries in [deg].
1167 :type region: :py:class:`tuple` of :py:class:`float`
1169 :returns: ``True``, if point is in region, else ``False``.
1170 :rtype: bool
1171 '''
1173 w, e, s, n = positive_region(region)
1174 return num.logical_and(
1175 num.logical_and(s <= p[0], p[0] <= n),
1176 num.logical_or(
1177 num.logical_and(w <= p[1], p[1] <= e),
1178 num.logical_and(w-360. <= p[1], p[1] <= e-360.)))
1181def radius_to_region(lat, lon, radius):
1182 '''
1183 Get a rectangular region which fully contains a given circular region.
1185 :param lat: Latitude of the center point of circular region in [deg].
1186 :type lat: float
1187 :param lon: Longitude of the center point of circular region in [deg].
1188 :type lon: float
1189 :param radius: Radius of circular region in [m].
1190 :type radius: float
1192 :returns: Rectangular region as ``(east, west, south, north)`` in [deg] or
1193 ``None``.
1194 :rtype: tuple[float, float, float, float]
1195 '''
1196 radius_deg = radius * m2d
1197 if radius_deg < 45.:
1198 lat_min = max(-90., lat - radius_deg)
1199 lat_max = min(90., lat + radius_deg)
1200 absmaxlat = max(abs(lat_min), abs(lat_max))
1201 if absmaxlat > 89:
1202 lon_min = -180.
1203 lon_max = 180.
1204 else:
1205 lon_min = max(
1206 -180. - 360.,
1207 lon - radius_deg / math.cos(absmaxlat*d2r))
1208 lon_max = min(
1209 180. + 360.,
1210 lon + radius_deg / math.cos(absmaxlat*d2r))
1212 lon_min, lon_max, lat_min, lat_max = positive_region(
1213 (lon_min, lon_max, lat_min, lat_max))
1215 return lon_min, lon_max, lat_min, lat_max
1217 else:
1218 return None
1221def geographic_midpoint(
1222 lats, lons, weights=None, depths=None, earthradius=earthradius):
1224 '''
1225 Calculate geographic midpoints by finding the center of gravity.
1227 This method suffers from instabilities if points are centered around the
1228 poles.
1230 :param lats: Latitudes in [deg].
1231 :param lons: Longitudes in [deg].
1232 :param weights: Weighting factors.
1233 :param depths: Depths in [m].
1234 :type lats: :py:class:`numpy.ndarray`, ``(N)``
1235 :type lons: :py:class:`numpy.ndarray`, ``(N)``
1236 :type weights: optional, :py:class:`numpy.ndarray`, ``(N)``
1237 :type depths: optional, :py:class:`numpy.ndarray`, ``(N)``
1239 :return: Latitudes and longitudes of the midpoint in [deg] (and depth [m]
1240 if depths are given).
1241 :rtype: ``(lat, lon)`` or ``(lat, lon, depth)``
1242 '''
1243 if not weights:
1244 weights = num.ones(len(lats))
1246 total_weigth = num.sum(weights)
1247 weights /= total_weigth
1248 lats = lats * d2r
1249 lons = lons * d2r
1250 if depths is not None:
1251 radii = (earthradius - depths) / earthradius
1252 else:
1253 radii = 1.0
1255 x = num.sum(num.cos(lats) * num.cos(lons) * weights * radii)
1256 y = num.sum(num.cos(lats) * num.sin(lons) * weights * radii)
1257 z = num.sum(num.sin(lats) * weights * radii)
1259 lon = num.arctan2(y, x)
1260 hyp = num.sqrt(x**2 + y**2)
1261 lat = num.arctan2(z, hyp)
1262 depth = earthradius - num.sqrt(x**2 + y**2 + z**2) * earthradius
1264 if depths is None:
1265 return lat/d2r, lon/d2r
1266 else:
1267 return lat/d2r, lon/d2r, depth
1270def geographic_midpoint_locations(
1271 locations, weights=None, include_depth=False):
1273 if not include_depth:
1274 coords = num.array([loc.effective_latlon
1275 for loc in locations])
1276 return geographic_midpoint(
1277 coords[:, 0], coords[:, 1], weights)
1278 else:
1279 coords = num.array([loc.effective_latlon + (loc.depth,)
1280 for loc in locations])
1281 return geographic_midpoint(
1282 coords[:, 0], coords[:, 1], weights=weights, depths=coords[:, 2])
1285def geodetic_to_ecef(lat, lon, alt):
1286 '''
1287 Convert geodetic coordinates to Earth-Centered, Earth-Fixed (ECEF)
1288 Cartesian coordinates. [#1]_ [#2]_
1290 :param lat: Geodetic latitude in [deg].
1291 :param lon: Geodetic longitude in [deg].
1292 :param alt: Geodetic altitude (height) in [m] (positive for points outside
1293 the geoid).
1294 :type lat: float
1295 :type lon: float
1296 :type alt: float
1298 :return: ECEF Cartesian coordinates (X, Y, Z) in [m].
1299 :rtype: tuple[float, float, float]
1301 .. [#1] https://en.wikipedia.org/wiki/ECEF
1302 .. [#2] https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1303 #From_geodetic_to_ECEF_coordinates
1304 '''
1306 f = earth_oblateness
1307 a = earthradius_equator
1308 e2 = 2*f - f**2
1310 lat, lon = num.radians(lat), num.radians(lon)
1311 # Normal (plumb line)
1312 N = a / num.sqrt(1.0 - (e2 * num.sin(lat)**2))
1314 X = (N+alt) * num.cos(lat) * num.cos(lon)
1315 Y = (N+alt) * num.cos(lat) * num.sin(lon)
1316 Z = (N*(1.0-e2) + alt) * num.sin(lat)
1318 return (X, Y, Z)
1321def ecef_to_geodetic(X, Y, Z):
1322 '''
1323 Convert Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates to
1324 geodetic coordinates (Ferrari's solution).
1326 :param X: ECEF X coordinate [m].
1327 :type X: float
1328 :param Y: ECEF Y coordinate [m].
1329 :type Y: float
1330 :param Z: ECEF Z coordinate [m].
1331 :type Z: float
1333 :return: Geodetic coordinates (lat, lon, alt). Latitude and longitude are
1334 in [deg] and altitude is in [m]
1335 (positive for points outside the geoid).
1336 :rtype: :py:class:`tuple` of :py:class:`float`
1338 .. seealso ::
1339 https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1340 #The_application_of_Ferrari.27s_solution
1341 '''
1342 f = earth_oblateness
1343 a = earthradius_equator
1344 b = a * (1. - f)
1345 e2 = 2.*f - f**2
1347 # usefull
1348 a2 = a**2
1349 b2 = b**2
1350 e4 = e2**2
1351 X2 = X**2
1352 Y2 = Y**2
1353 Z2 = Z**2
1355 r = num.sqrt(X2 + Y2)
1356 r2 = r**2
1358 e_prime2 = (a2 - b2)/b2
1359 E2 = a2 - b2
1360 F = 54. * b2 * Z2
1361 G = r2 + (1.-e2)*Z2 - (e2*E2)
1362 C = (e4 * F * r2) / (G**3)
1363 S = cbrt(1. + C + num.sqrt(C**2 + 2.*C))
1364 P = F / (3. * (S + 1./S + 1.)**2 * G**2)
1365 Q = num.sqrt(1. + (2.*e4*P))
1367 dum1 = -(P*e2*r) / (1.+Q)
1368 dum2 = 0.5 * a2 * (1. + 1./Q)
1369 dum3 = (P * (1.-e2) * Z2) / (Q * (1.+Q))
1370 dum4 = 0.5 * P * r2
1371 r0 = dum1 + num.sqrt(dum2 - dum3 - dum4)
1373 U = num.sqrt((r - e2*r0)**2 + Z2)
1374 V = num.sqrt((r - e2*r0)**2 + (1.-e2)*Z2)
1375 Z0 = (b2*Z) / (a*V)
1377 alt = U * (1. - (b2 / (a*V)))
1378 lat = num.arctan((Z + e_prime2 * Z0)/r)
1379 lon = num.arctan2(Y, X)
1381 return (lat*r2d, lon*r2d, alt)
1384class Farside(Exception):
1385 pass
1388def latlon_to_xyz(latlons):
1389 if latlons.ndim == 1:
1390 return latlon_to_xyz(latlons[num.newaxis, :])[0]
1392 points = num.zeros((latlons.shape[0], 3))
1393 lats = latlons[:, 0]
1394 lons = latlons[:, 1]
1395 points[:, 0] = num.cos(lats*d2r) * num.cos(lons*d2r)
1396 points[:, 1] = num.cos(lats*d2r) * num.sin(lons*d2r)
1397 points[:, 2] = num.sin(lats*d2r)
1398 return points
1401def xyz_to_latlon(xyz):
1402 if xyz.ndim == 1:
1403 return xyz_to_latlon(xyz[num.newaxis, :])[0]
1405 latlons = num.zeros((xyz.shape[0], 2))
1406 latlons[:, 0] = num.arctan2(
1407 xyz[:, 2], num.sqrt(xyz[:, 0]**2 + xyz[:, 1]**2)) * r2d
1408 latlons[:, 1] = num.arctan2(
1409 xyz[:, 1], xyz[:, 0]) * r2d
1410 return latlons
1413def rot_to_00(lat, lon):
1414 rot0 = euler_to_matrix(0., -90.*d2r, 0.)
1415 rot1 = euler_to_matrix(-d2r*lat, 0., -d2r*lon)
1416 return num.dot(rot0.T, num.dot(rot1, rot0)).T
1419def distances3d(a, b):
1420 return num.sqrt(num.sum((a-b)**2, axis=a.ndim-1))
1423def circulation(points2):
1424 return num.sum(
1425 (points2[1:, 0] - points2[:-1, 0])
1426 * (points2[1:, 1] + points2[:-1, 1]))
1429def stereographic(points):
1430 dists = distances3d(points[1:, :], points[:-1, :])
1431 if dists.size > 0:
1432 maxdist = num.max(dists)
1433 cutoff = maxdist**2 / 2.
1434 else:
1435 cutoff = 1.0e-5
1437 points = points.copy()
1438 if num.any(points[:, 0] < -1. + cutoff):
1439 raise Farside()
1441 points_out = points[:, 1:].copy()
1442 factor = 1.0 / (1.0 + points[:, 0])
1443 points_out *= factor[:, num.newaxis]
1445 return points_out
1448def stereographic_poly(points):
1449 dists = distances3d(points[1:, :], points[:-1, :])
1450 if dists.size > 0:
1451 maxdist = num.max(dists)
1452 cutoff = maxdist**2 / 2.
1453 else:
1454 cutoff = 1.0e-5
1456 points = points.copy()
1457 if num.any(points[:, 0] < -1. + cutoff):
1458 raise Farside()
1460 points_out = points[:, 1:].copy()
1461 factor = 1.0 / (1.0 + points[:, 0])
1462 points_out *= factor[:, num.newaxis]
1464 if circulation(points_out) >= 0:
1465 raise Farside()
1467 return points_out
1470def gnomonic_x(points, cutoff=0.01):
1471 points_out = points[:, 1:].copy()
1472 if num.any(points[:, 0] < cutoff):
1473 raise Farside()
1475 factor = 1.0 / points[:, 0]
1476 points_out *= factor[:, num.newaxis]
1477 return points_out
1480def cneg(i, x):
1481 if i == 1:
1482 return x
1483 else:
1484 return num.logical_not(x)
1487def contains_points(polygon, points):
1488 '''
1489 Test which points are inside polygon on a sphere.
1491 The inside of the polygon is defined as the area which is to the left hand
1492 side of an observer walking the polygon line, points in order, on the
1493 sphere. Lines between the polygon points are treated as great circle paths.
1494 The polygon may be arbitrarily complex, as long as it does not have any
1495 crossings or thin parts with zero width. The polygon may contain the poles
1496 and is allowed to wrap around the sphere multiple times.
1498 The algorithm works by consecutive cutting of the polygon into (almost)
1499 hemispheres and subsequent Gnomonic projections to perform the
1500 point-in-polygon tests on a 2D plane.
1502 :param polygon: Point coordinates defining the polygon [deg].
1503 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1504 0=lat, 1=lon
1505 :param points: Coordinates of points to test [deg].
1506 :type points: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1507 0=lat, 1=lon
1509 :returns: Boolean mask array.
1510 :rtype: :py:class:`numpy.ndarray` of shape ``(N,)``.
1511 '''
1513 and_ = num.logical_and
1515 points_xyz = latlon_to_xyz(points)
1516 mask_x = 0. <= points_xyz[:, 0]
1517 mask_y = 0. <= points_xyz[:, 1]
1518 mask_z = 0. <= points_xyz[:, 2]
1520 result = num.zeros(points.shape[0], dtype=int)
1522 for ix in [-1, 1]:
1523 for iy in [-1, 1]:
1524 for iz in [-1, 1]:
1525 mask = and_(
1526 and_(cneg(ix, mask_x), cneg(iy, mask_y)),
1527 cneg(iz, mask_z))
1529 center_xyz = num.array([ix, iy, iz], dtype=float)
1531 lat, lon = xyz_to_latlon(center_xyz)
1532 rot = rot_to_00(lat, lon)
1534 points_rot_xyz = num.dot(rot, points_xyz[mask, :].T).T
1535 points_rot_pro = gnomonic_x(points_rot_xyz)
1537 offset = 0.01
1539 poly_xyz = latlon_to_xyz(polygon)
1540 poly_rot_xyz = num.dot(rot, poly_xyz.T).T
1541 poly_rot_xyz[:, 0] -= offset
1542 groups = spoly_cut([poly_rot_xyz], axis=0)
1544 for poly_rot_group_xyz in groups[1]:
1545 poly_rot_group_xyz[:, 0] += offset
1547 poly_rot_group_pro = gnomonic_x(
1548 poly_rot_group_xyz)
1550 if circulation(poly_rot_group_pro) > 0:
1551 result[mask] += path_contains_points(
1552 poly_rot_group_pro, points_rot_pro)
1553 else:
1554 result[mask] -= path_contains_points(
1555 poly_rot_group_pro, points_rot_pro)
1557 return result.astype(bool)
1560def contains_point(polygon, point):
1561 '''
1562 Test if point is inside polygon on a sphere.
1564 Convenience wrapper to :py:func:`contains_points` to test a single point.
1566 :param polygon: Point coordinates defining the polygon [deg].
1567 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1568 0=lat, 1=lon
1569 :param point: Coordinates ``(lat, lon)`` of point to test [deg].
1570 :type point: :py:class:`tuple` of :py:class:`float`
1572 :returns: ``True``, if point is located within polygon, else ``False``.
1573 :rtype: bool
1574 '''
1576 return bool(
1577 contains_points(polygon, num.asarray(point)[num.newaxis, :])[0])