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