1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import division, absolute_import
7import math
8import numpy as num
10from .moment_tensor import euler_to_matrix
11from .config import config
12from .plot.beachball import spoly_cut
14from matplotlib.path import Path
16d2r = math.pi/180.
17r2d = 1./d2r
18earth_oblateness = 1./298.257223563
19earthradius_equator = 6378.14 * 1000.
20earthradius = config().earthradius
21d2m = earthradius_equator*math.pi/180.
22m2d = 1./d2m
24_testpath = Path([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)], closed=True)
26raise_if_slow_path_contains_points = False
29class Slow(Exception):
30 pass
33if hasattr(_testpath, 'contains_points') and num.all(
34 _testpath.contains_points([(0.5, 0.5), (1.5, 0.5)]) == [True, False]):
36 def path_contains_points(verts, points):
37 p = Path(verts, closed=True)
38 return p.contains_points(points).astype(num.bool)
40else:
41 # work around missing contains_points and bug in matplotlib ~ v1.2.0
43 def path_contains_points(verts, points):
44 if raise_if_slow_path_contains_points:
45 # used by unit test to skip slow gshhg_example.py
46 raise Slow()
48 p = Path(verts, closed=True)
49 result = num.zeros(points.shape[0], dtype=num.bool)
50 for i in range(result.size):
51 result[i] = p.contains_point(points[i, :])
53 return result
56try:
57 cbrt = num.cbrt
58except AttributeError:
59 def cbrt(x):
60 return x**(1./3.)
63def float_array_broadcast(*args):
64 return num.broadcast_arrays(*[
65 num.asarray(x, dtype=float) for x in args])
68class Loc(object):
69 '''
70 Simple location representation.
72 :attrib lat: Latitude degree
73 :attrib lon: Longitude degree
74 '''
75 def __init__(self, lat, lon):
76 self.lat = lat
77 self.lon = lon
80def clip(x, mi, ma):
81 '''
82 Clip values of an array.
84 :param x: Continunous data to be clipped
85 :param mi: Clip minimum
86 :param ma: Clip maximum
87 :type x: :py:class:`numpy.ndarray`
88 :type mi: float
89 :type ma: float
91 :return: Clipped data
92 :rtype: :py:class:`numpy.ndarray`
93 '''
94 return num.minimum(num.maximum(mi, x), ma)
97def wrap(x, mi, ma):
98 '''
99 Wrapping continuous data to fundamental phase values.
101 .. math::
102 x_{\\mathrm{wrapped}} = x_{\\mathrm{cont},i} -
103 \\frac{ x_{\\mathrm{cont},i} - r_{\\mathrm{min}} }
104 { r_{\\mathrm{max}} - r_{\\mathrm{min}}}
105 \\cdot ( r_{\\mathrm{max}} - r_{\\mathrm{min}}),\\quad
106 x_{\\mathrm{wrapped}}\\; \\in
107 \\;[ r_{\\mathrm{min}},\\, r_{\\mathrm{max}}].
109 :param x: Continunous data to be wrapped
110 :param mi: Minimum value of wrapped data
111 :param ma: Maximum value of wrapped data
112 :type x: :py:class:`numpy.ndarray`
113 :type mi: float
114 :type ma: float
116 :return: Wrapped data
117 :rtype: :py:class:`numpy.ndarray`
118 '''
119 return x - num.floor((x-mi)/(ma-mi)) * (ma-mi)
122def _latlon_pair(args):
123 if len(args) == 2:
124 a, b = args
125 return a.lat, a.lon, b.lat, b.lon
127 elif len(args) == 4:
128 return args
131def cosdelta(*args):
132 '''
133 Cosine of the angular distance between two points ``a`` and ``b`` on a
134 sphere.
136 This function (find implementation below) returns the cosine of the
137 distance angle 'delta' between two points ``a`` and ``b``, coordinates of
138 which are expected to be given in geographical coordinates and in degrees.
139 For numerical stability a maximum of 1.0 is enforced.
141 .. math::
143 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
144 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
145 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
146 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\[0.5cm]
148 \\cos(\\Delta) = \\min( 1.0, \\quad \\sin( A_{\\mathrm{lat'}})
149 \\sin( B_{\\mathrm{lat'}} ) +
150 \\cos(A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
151 \\cos( B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )
153 :param a: Location point A
154 :type a: :py:class:`pyrocko.orthodrome.Loc`
155 :param b: Location point B
156 :type b: :py:class:`pyrocko.orthodrome.Loc`
158 :return: cosdelta
159 :rtype: float
160 '''
162 alat, alon, blat, blon = _latlon_pair(args)
164 return min(
165 1.0,
166 math.sin(alat*d2r) * math.sin(blat*d2r) +
167 math.cos(alat*d2r) * math.cos(blat*d2r) *
168 math.cos(d2r*(blon-alon)))
171def cosdelta_numpy(a_lats, a_lons, b_lats, b_lons):
172 '''
173 Cosine of the angular distance between two points ``a`` and ``b`` on a
174 sphere.
176 This function returns the cosines of the distance
177 angles *delta* between two points ``a`` and ``b`` given as
178 :py:class:`numpy.ndarray`.
179 The coordinates are expected to be given in geographical coordinates
180 and in degrees. For numerical stability a maximum of ``1.0`` is enforced.
182 Please find the details of the implementation in the documentation of
183 the function :py:func:`pyrocko.orthodrome.cosdelta` above.
185 :param a_lats: Latitudes (degree) point A
186 :param a_lons: Longitudes (degree) point A
187 :param b_lats: Latitudes (degree) point B
188 :param b_lons: Longitudes (degree) point B
189 :type a_lats: :py:class:`numpy.ndarray`
190 :type a_lons: :py:class:`numpy.ndarray`
191 :type b_lats: :py:class:`numpy.ndarray`
192 :type b_lons: :py:class:`numpy.ndarray`
194 :return: cosdelta
195 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
196 '''
197 return num.minimum(
198 1.0,
199 num.sin(a_lats*d2r) * num.sin(b_lats*d2r) +
200 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
201 num.cos(d2r*(b_lons-a_lons)))
204def azimuth(*args):
205 '''
206 Azimuth calculation.
208 This function (find implementation below) returns azimuth ...
209 between points ``a`` and ``b``, coordinates of
210 which are expected to be given in geographical coordinates and in degrees.
212 .. math::
214 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
215 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
216 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
217 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\
219 \\varphi_{\\mathrm{azi},AB} = \\frac{180}{\\pi} \\arctan \\left[
220 \\frac{
221 \\cos( A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
222 \\sin(B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )}
223 {\\sin ( B_{\\mathrm{lat'}} ) - \\sin( A_{\\mathrm{lat'}}
224 cosdelta) } \\right]
226 :param a: Location point A
227 :type a: :py:class:`pyrocko.orthodrome.Loc`
228 :param b: Location point B
229 :type b: :py:class:`pyrocko.orthodrome.Loc`
231 :return: Azimuth in degree
232 '''
234 alat, alon, blat, blon = _latlon_pair(args)
236 return r2d*math.atan2(
237 math.cos(alat*d2r) * math.cos(blat*d2r) *
238 math.sin(d2r*(blon-alon)),
239 math.sin(d2r*blat) - math.sin(d2r*alat) * cosdelta(
240 alat, alon, blat, blon))
243def azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta=None):
244 '''
245 Calculation of the azimuth (*track angle*) from a location A towards B.
247 This function returns azimuths (*track angles*) from locations A towards B
248 given in :py:class:`numpy.ndarray`. Coordinates are expected to be given in
249 geographical coordinates and in degrees.
251 Please find the details of the implementation in the documentation of the
252 function :py:func:`pyrocko.orthodrome.azimuth`.
255 :param a_lats: Latitudes (degree) point A
256 :param a_lons: Longitudes (degree) point A
257 :param b_lats: Latitudes (degree) point B
258 :param b_lons: Longitudes (degree) point B
259 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
260 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
261 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
262 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
264 :return: Azimuths in degrees
265 :rtype: :py:class:`numpy.ndarray`, ``(N)``
266 '''
267 if _cosdelta is None:
268 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
270 return r2d*num.arctan2(
271 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
272 num.sin(d2r*(b_lons-a_lons)),
273 num.sin(d2r*b_lats) - num.sin(d2r*a_lats) * _cosdelta)
276def azibazi(*args, **kwargs):
277 alat, alon, blat, blon = _latlon_pair(args)
278 if alat == blat and alon == blon:
279 return 0., 180.
281 implementation = kwargs.get('implementation', 'c')
282 assert implementation in ('c', 'python')
283 if implementation == 'c':
284 from pyrocko import orthodrome_ext
285 return orthodrome_ext.azibazi(alat, alon, blat, blon)
287 cd = cosdelta(alat, alon, blat, blon)
288 azi = r2d*math.atan2(
289 math.cos(alat*d2r) * math.cos(blat*d2r) *
290 math.sin(d2r*(blon-alon)),
291 math.sin(d2r*blat) - math.sin(d2r*alat) * cd)
292 bazi = r2d*math.atan2(
293 math.cos(blat*d2r) * math.cos(alat*d2r) *
294 math.sin(d2r*(alon-blon)),
295 math.sin(d2r*alat) - math.sin(d2r*blat) * cd)
297 return azi, bazi
300def azibazi_numpy(a_lats, a_lons, b_lats, b_lons, implementation='c'):
302 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
303 a_lats, a_lons, b_lats, b_lons)
305 assert implementation in ('c', 'python')
306 if implementation == 'c':
307 from pyrocko import orthodrome_ext
308 return orthodrome_ext.azibazi_numpy(a_lats, a_lons, b_lats, b_lons)
310 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
311 azis = azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta)
312 bazis = azimuth_numpy(b_lats, b_lons, a_lats, a_lons, _cosdelta)
314 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
315 ii_eq = num.where(eq)[0]
316 azis[ii_eq] = 0.0
317 bazis[ii_eq] = 180.0
318 return azis, bazis
321def azidist_numpy(*args):
322 '''
323 Calculation of the azimuth (*track angle*) and the distance from locations
324 A towards B on a sphere.
326 The assisting functions used are :py:func:`pyrocko.orthodrome.cosdelta` and
327 :py:func:`pyrocko.orthodrome.azimuth`
329 :param a_lats: Latitudes (degree) point A
330 :param a_lons: Longitudes (degree) point A
331 :param b_lats: Latitudes (degree) point B
332 :param b_lons: Longitudes (degree) point B
333 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
334 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
335 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
336 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
338 :return: Azimuths in degrees, distances in degrees
339 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
340 '''
341 _cosdelta = cosdelta_numpy(*args)
342 _azimuths = azimuth_numpy(_cosdelta=_cosdelta, *args)
343 return _azimuths, r2d*num.arccos(_cosdelta)
346def distance_accurate50m(*args, **kwargs):
347 '''
348 Accurate distance calculation based on a spheroid of rotation.
350 Function returns distance in meter between points A and B, coordinates of
351 which must be given in geographical coordinates and in degrees.
352 The returned distance should be accurate to 50 m using WGS84.
353 Values for the Earth's equator radius and the Earth's oblateness
354 (``f_oblate``) are defined in the pyrocko configuration file
355 :py:class:`pyrocko.config`.
357 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
359 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
360 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
362 .. math::
364 F = \\frac{\\pi}{180}
365 \\frac{(A_{lat} + B_{lat})}{2}, \\quad
366 G = \\frac{\\pi}{180}
367 \\frac{(A_{lat} - B_{lat})}{2}, \\quad
368 l = \\frac{\\pi}{180}
369 \\frac{(A_{lon} - B_{lon})}{2} \\quad
370 \\\\[0.5cm]
371 S = \\sin^2(G) \\cdot \\cos^2(l) +
372 \\cos^2(F) \\cdot \\sin^2(l), \\quad \\quad
373 C = \\cos^2(G) \\cdot \\cos^2(l) +
374 \\sin^2(F) \\cdot \\sin^2(l)
376 .. math::
378 w = \\arctan \\left( \\sqrt{ \\frac{S}{C}} \\right) , \\quad
379 r = \\sqrt{\\frac{S}{C} }
381 The spherical-earth distance D between A and B, can be given with:
383 .. math::
385 D_{sphere} = 2w \\cdot R_{equator}
387 The oblateness of the Earth requires some correction with
388 correction factors h1 and h2:
390 .. math::
392 h_1 = \\frac{3r - 1}{2C}, \\quad
393 h_2 = \\frac{3r +1 }{2S}\\\\[0.5cm]
395 D = D_{\\mathrm{sphere}} \\cdot [ 1 + h_1 \\,f_{\\mathrm{oblate}}
396 \\cdot \\sin^2(F)
397 \\cos^2(G) - h_2\\, f_{\\mathrm{oblate}}
398 \\cdot \\cos^2(F) \\sin^2(G)]
401 :param a: Location point A
402 :type a: :py:class:`pyrocko.orthodrome.Loc`
403 :param b: Location point B
404 :type b: :py:class:`pyrocko.orthodrome.Loc`
406 :return: Distance in meter
407 :rtype: float
408 '''
410 alat, alon, blat, blon = _latlon_pair(args)
412 implementation = kwargs.get('implementation', 'c')
413 assert implementation in ('c', 'python')
414 if implementation == 'c':
415 from pyrocko import orthodrome_ext
416 return orthodrome_ext.distance_accurate50m(alat, alon, blat, blon)
418 f = (alat + blat)*d2r / 2.
419 g = (alat - blat)*d2r / 2.
420 h = (alon - blon)*d2r / 2.
422 s = math.sin(g)**2 * math.cos(h)**2 + math.cos(f)**2 * math.sin(h)**2
423 c = math.cos(g)**2 * math.cos(h)**2 + math.sin(f)**2 * math.sin(h)**2
425 w = math.atan(math.sqrt(s/c))
427 if w == 0.0:
428 return 0.0
430 r = math.sqrt(s*c)/w
431 d = 2.*w*earthradius_equator
432 h1 = (3.*r-1.)/(2.*c)
433 h2 = (3.*r+1.)/(2.*s)
435 return d * (1. +
436 earth_oblateness * h1 * math.sin(f)**2 * math.cos(g)**2 -
437 earth_oblateness * h2 * math.cos(f)**2 * math.sin(g)**2)
440def distance_accurate50m_numpy(
441 a_lats, a_lons, b_lats, b_lons, implementation='c'):
443 '''
444 Accurate distance calculation based on a spheroid of rotation.
446 Function returns distance in meter between points ``a`` and ``b``,
447 coordinates of which must be given in geographical coordinates and in
448 degrees.
449 The returned distance should be accurate to 50 m using WGS84.
450 Values for the Earth's equator radius and the Earth's oblateness
451 (``f_oblate``) are defined in the pyrocko configuration file
452 :py:class:`pyrocko.config`.
454 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
456 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
457 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
459 .. math::
461 F_i = \\frac{\\pi}{180}
462 \\frac{(a_{lat,i} + a_{lat,i})}{2}, \\quad
463 G_i = \\frac{\\pi}{180}
464 \\frac{(a_{lat,i} - b_{lat,i})}{2}, \\quad
465 l_i= \\frac{\\pi}{180}
466 \\frac{(a_{lon,i} - b_{lon,i})}{2} \\\\[0.5cm]
467 S_i = \\sin^2(G_i) \\cdot \\cos^2(l_i) +
468 \\cos^2(F_i) \\cdot \\sin^2(l_i), \\quad \\quad
469 C_i = \\cos^2(G_i) \\cdot \\cos^2(l_i) +
470 \\sin^2(F_i) \\cdot \\sin^2(l_i)
472 .. math::
474 w_i = \\arctan \\left( \\sqrt{\\frac{S_i}{C_i}} \\right), \\quad
475 r_i = \\sqrt{\\frac{S_i}{C_i} }
477 The spherical-earth distance ``D`` between ``a`` and ``b``,
478 can be given with:
480 .. math::
482 D_{\\mathrm{sphere},i} = 2w_i \\cdot R_{\\mathrm{equator}}
484 The oblateness of the Earth requires some correction with
485 correction factors ``h1`` and ``h2``:
487 .. math::
489 h_{1.i} = \\frac{3r - 1}{2C_i}, \\quad
490 h_{2,i} = \\frac{3r +1 }{2S_i}\\\\[0.5cm]
492 D_{AB,i} = D_{\\mathrm{sphere},i} \\cdot [1 + h_{1,i}
493 \\,f_{\\mathrm{oblate}}
494 \\cdot \\sin^2(F_i)
495 \\cos^2(G_i) - h_{2,i}\\, f_{\\mathrm{oblate}}
496 \\cdot \\cos^2(F_i) \\sin^2(G_i)]
498 :param a_lats: Latitudes (degree) point A
499 :param a_lons: Longitudes (degree) point A
500 :param b_lats: Latitudes (degree) point B
501 :param b_lons: Longitudes (degree) point B
502 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
503 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
504 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
505 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
507 :return: Distances in meter
508 :rtype: :py:class:`numpy.ndarray`, ``(N)``
509 '''
511 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
512 a_lats, a_lons, b_lats, b_lons)
514 assert implementation in ('c', 'python')
515 if implementation == 'c':
516 from pyrocko import orthodrome_ext
517 return orthodrome_ext.distance_accurate50m_numpy(
518 a_lats, a_lons, b_lats, b_lons)
520 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
521 ii_neq = num.where(num.logical_not(eq))[0]
523 if num.all(eq):
524 return num.zeros_like(eq, dtype=float)
526 def extr(x):
527 if isinstance(x, num.ndarray) and x.size > 1:
528 return x[ii_neq]
529 else:
530 return x
532 a_lats = extr(a_lats)
533 a_lons = extr(a_lons)
534 b_lats = extr(b_lats)
535 b_lons = extr(b_lons)
537 f = (a_lats + b_lats)*d2r / 2.
538 g = (a_lats - b_lats)*d2r / 2.
539 h = (a_lons - b_lons)*d2r / 2.
541 s = num.sin(g)**2 * num.cos(h)**2 + num.cos(f)**2 * num.sin(h)**2
542 c = num.cos(g)**2 * num.cos(h)**2 + num.sin(f)**2 * num.sin(h)**2
544 w = num.arctan(num.sqrt(s/c))
546 r = num.sqrt(s*c)/w
548 d = 2.*w*earthradius_equator
549 h1 = (3.*r-1.)/(2.*c)
550 h2 = (3.*r+1.)/(2.*s)
552 dists = num.zeros(eq.size, dtype=float)
553 dists[ii_neq] = d * (
554 1. +
555 earth_oblateness * h1 * num.sin(f)**2 * num.cos(g)**2 -
556 earth_oblateness * h2 * num.cos(f)**2 * num.sin(g)**2)
558 return dists
561def ne_to_latlon(lat0, lon0, north_m, east_m):
562 '''
563 Transform local cartesian coordinates to latitude and longitude.
565 From east and north coordinates (``x`` and ``y`` coordinate
566 :py:class:`numpy.ndarray`) relative to a reference differences in
567 longitude and latitude are calculated, which are effectively changes in
568 azimuth and distance, respectively:
570 .. math::
572 \\text{distance change:}\\; \\Delta {\\bf{a}} &= \\sqrt{{\\bf{y}}^2 +
573 {\\bf{x}}^2 }/ \\mathrm{R_E},
575 \\text{azimuth change:}\\; \\Delta \\bf{\\gamma} &= \\arctan( \\bf{x}
576 / \\bf{y}).
578 The projection used preserves the azimuths of the input points.
580 :param lat0: Latitude origin of the cartesian coordinate system.
581 :param lon0: Longitude origin of the cartesian coordinate system.
582 :param north_m: Northing distances from origin in meters.
583 :param east_m: Easting distances from origin in meters.
584 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
585 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
586 :type lat0: float
587 :type lon0: float
589 :return: Array with latitudes and longitudes
590 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
592 '''
594 a = num.sqrt(north_m**2+east_m**2)/earthradius
595 gamma = num.arctan2(east_m, north_m)
597 return azidist_to_latlon_rad(lat0, lon0, gamma, a)
600def azidist_to_latlon(lat0, lon0, azimuth_deg, distance_deg):
601 '''
602 (Durchreichen??).
603 '''
605 return azidist_to_latlon_rad(
606 lat0, lon0, azimuth_deg/180.*num.pi, distance_deg/180.*num.pi)
609def azidist_to_latlon_rad(lat0, lon0, azimuth_rad, distance_rad):
610 '''
611 Absolute latitudes and longitudes are calculated from relative changes.
613 For numerical stability a range between of ``-1.0`` and ``1.0`` is
614 enforced for ``c`` and ``alpha``.
616 .. math::
618 \\Delta {\\bf a}_i \\; \\text{and} \\; \\Delta \\gamma_i \\;
619 \\text{are relative distances and azimuths from lat0 and lon0 for
620 \\textit{i} source points of a finite source.}
622 .. math::
624 \\mathrm{b} &= \\frac{\\pi}{2} -\\frac{\\pi}{180}\\;\\mathrm{lat_0}\\\\
625 {\\bf c}_i &=\\arccos[\\; \\cos(\\Delta {\\bf{a}}_i)
626 \\cos(\\mathrm{b}) + |\\Delta \\gamma_i| \\,
627 \\sin(\\Delta {\\bf a}_i)
628 \\sin(\\mathrm{b})\\; ] \\\\
629 \\mathrm{lat}_i &= \\frac{180}{\\pi}
630 \\left(\\frac{\\pi}{2} - {\\bf c}_i \\right)
632 .. math::
634 \\alpha_i &= \\arcsin \\left[ \\; \\frac{ \\sin(\\Delta {\\bf a}_i )
635 \\sin(|\\Delta \\gamma_i|)}{\\sin({\\bf c}_i)}\\;
636 \\right] \\\\
637 \\alpha_i &= \\begin{cases}
638 \\alpha_i, &\\text{if} \\; \\cos(\\Delta {\\bf a}_i) -
639 \\cos(\\mathrm{b}) \\cos({\\bf{c}}_i) > 0, \\;
640 \\text{else} \\\\
641 \\pi - \\alpha_i, & \\text{if} \\; \\alpha_i > 0,\\;
642 \\text{else}\\\\
643 -\\pi - \\alpha_i, & \\text{if} \\; \\alpha_i < 0.
644 \\end{cases} \\\\
645 \\mathrm{lon}_i &= \\mathrm{lon_0} +
646 \\frac{180}{\\pi} \\,
647 \\frac{\\Delta \\gamma_i }{|\\Delta \\gamma_i|}
648 \\cdot \\alpha_i
649 \\text{, with $\\alpha_i \\in [-\\pi,\\pi]$}
651 :param lat0: Latitude origin of the cartesian coordinate system.
652 :param lon0: Longitude origin of the cartesian coordinate system.
653 :param distance_rad: Distances from origin in radians.
654 :param azimuth_rad: Azimuth from radians.
655 :type distance_rad: :py:class:`numpy.ndarray`, ``(N)``
656 :type azimuth_rad: :py:class:`numpy.ndarray`, ``(N)``
657 :type lat0: float
658 :type lon0: float
660 :return: Array with latitudes and longitudes
661 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
662 '''
664 a = distance_rad
665 gamma = azimuth_rad
667 b = math.pi/2.-lat0*d2r
669 alphasign = 1.
670 alphasign = num.where(gamma < 0, -1., 1.)
671 gamma = num.abs(gamma)
673 c = num.arccos(clip(
674 num.cos(a)*num.cos(b)+num.sin(a)*num.sin(b)*num.cos(gamma), -1., 1.))
676 alpha = num.arcsin(clip(
677 num.sin(a)*num.sin(gamma)/num.sin(c), -1., 1.))
679 alpha = num.where(
680 num.cos(a)-num.cos(b)*num.cos(c) < 0,
681 num.where(alpha > 0, math.pi-alpha, -math.pi-alpha),
682 alpha)
684 lat = r2d * (math.pi/2. - c)
685 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
687 return lat, lon
690def ne_to_latlon_alternative_method(lat0, lon0, north_m, east_m):
691 '''
692 Transform local cartesian coordinates to latitude and longitude.
694 Like :py:func:`pyrocko.orthodrome.ne_to_latlon`,
695 but this method (implementation below), although it should be numerically
696 more stable, suffers problems at points which are *across the pole*
697 as seen from the cartesian origin.
699 .. math::
701 \\text{distance change:}\\; \\Delta {{\\bf a}_i} &=
702 \\sqrt{{\\bf{y}}^2_i + {\\bf{x}}^2_i }/ \\mathrm{R_E},\\\\
703 \\text{azimuth change:}\\; \\Delta {\\bf \\gamma}_i &=
704 \\arctan( {\\bf x}_i {\\bf y}_i). \\\\
705 \\mathrm{b} &=
706 \\frac{\\pi}{2} -\\frac{\\pi}{180} \\;\\mathrm{lat_0}\\\\
708 .. math::
710 {{\\bf z}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i -
711 \\mathrm{b}}{2} \\right)}
712 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
713 {{\\bf n}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i +
714 \\mathrm{b}}{2} \\right)}
715 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
716 {{\\bf z}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i -
717 \\mathrm{b}}{2} \\right)}
718 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
719 {{\\bf n}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i +
720 \\mathrm{b}}{2} \\right)}
721 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
722 {{\\bf t}_1}_i &= \\arctan{\\left( \\frac{{{\\bf z}_1}_i}
723 {{{\\bf n}_1}_i} \\right) }\\\\
724 {{\\bf t}_2}_i &= \\arctan{\\left( \\frac{{{\\bf z}_2}_i}
725 {{{\\bf n}_2}_i} \\right) } \\\\[0.5cm]
726 c &= \\begin{cases}
727 2 \\cdot \\arccos \\left( {{\\bf z}_1}_i / \\sin({{\\bf t}_1}_i)
728 \\right),\\; \\text{if }
729 |\\sin({{\\bf t}_1}_i)| >
730 |\\sin({{\\bf t}_2}_i)|,\\; \\text{else} \\\\
731 2 \\cdot \\arcsin{\\left( {{\\bf z}_2}_i /
732 \\sin({{\\bf t}_2}_i) \\right)}.
733 \\end{cases}\\\\
735 .. math::
737 {\\bf {lat}}_i &= \\frac{180}{ \\pi } \\left( \\frac{\\pi}{2}
738 - {\\bf {c}}_i \\right) \\\\
739 {\\bf {lon}}_i &= {\\bf {lon}}_0 + \\frac{180}{ \\pi }
740 \\frac{\\gamma_i}{|\\gamma_i|},
741 \\text{ with}\\; \\gamma \\in [-\\pi,\\pi]
743 :param lat0: Latitude origin of the cartesian coordinate system.
744 :param lon0: Longitude origin of the cartesian coordinate system.
745 :param north_m: Northing distances from origin in meters.
746 :param east_m: Easting distances from origin in meters.
747 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
748 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
749 :type lat0: float
750 :type lon0: float
752 :return: Array with latitudes and longitudes
753 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
754 '''
756 b = math.pi/2.-lat0*d2r
757 a = num.sqrt(north_m**2+east_m**2)/earthradius
759 gamma = num.arctan2(east_m, north_m)
760 alphasign = 1.
761 alphasign = num.where(gamma < 0., -1., 1.)
762 gamma = num.abs(gamma)
764 z1 = num.cos((a-b)/2.)*num.cos(gamma/2.)
765 n1 = num.cos((a+b)/2.)*num.sin(gamma/2.)
766 z2 = num.sin((a-b)/2.)*num.cos(gamma/2.)
767 n2 = num.sin((a+b)/2.)*num.sin(gamma/2.)
768 t1 = num.arctan2(z1, n1)
769 t2 = num.arctan2(z2, n2)
771 alpha = t1 + t2
773 sin_t1 = num.sin(t1)
774 sin_t2 = num.sin(t2)
775 c = num.where(
776 num.abs(sin_t1) > num.abs(sin_t2),
777 num.arccos(z1/sin_t1)*2.,
778 num.arcsin(z2/sin_t2)*2.)
780 lat = r2d * (math.pi/2. - c)
781 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
782 return lat, lon
785def latlon_to_ne(*args):
786 '''
787 Relative cartesian coordinates with respect to a reference location.
789 For two locations, a reference location A and another location B, given in
790 geographical coordinates in degrees, the corresponding cartesian
791 coordinates are calculated.
792 Assisting functions are :py:func:`pyrocko.orthodrome.azimuth` and
793 :py:func:`pyrocko.orthodrome.distance_accurate50m`.
795 .. math::
797 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
798 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
799 \\mathrm{)}\\\\[0.3cm]
801 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180}
802 \\varphi_{\\mathrm{azi},AB} )\\\\
803 e &= D_{AB} \\cdot
804 \\sin( \\frac{\\pi }{180} \\varphi_{\\mathrm{azi},AB})
806 :param refloc: Location reference point
807 :type refloc: :py:class:`pyrocko.orthodrome.Loc`
808 :param loc: Location of interest
809 :type loc: :py:class:`pyrocko.orthodrome.Loc`
811 :return: Northing and easting from refloc to location
812 :rtype: tuple, float
814 '''
816 azi = azimuth(*args)
817 dist = distance_accurate50m(*args)
818 n, e = math.cos(azi*d2r)*dist, math.sin(azi*d2r)*dist
819 return n, e
822def latlon_to_ne_numpy(lat0, lon0, lat, lon):
823 '''
824 Relative cartesian coordinates with respect to a reference location.
826 For two locations, a reference location (``lat0``, ``lon0``) and another
827 location B, given in geographical coordinates in degrees,
828 the corresponding cartesian coordinates are calculated.
829 Assisting functions are :py:func:`azimuth`
830 and :py:func:`distance_accurate50m`.
832 :param lat0: reference location latitude
833 :param lon0: reference location longitude
834 :param lat: absolute location latitude
835 :param lon: absolute location longitude
837 :return: ``(n, e)``: relative north and east positions
838 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
840 Implemented formulations:
842 .. math::
844 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
845 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
846 \\mathrm{)}\\\\[0.3cm]
848 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180} \\varphi_{
849 \\mathrm{azi},AB} )\\\\
850 e &= D_{AB} \\cdot \\sin( \\frac{\\pi }{180} \\varphi_{
851 \\mathrm{azi},AB} )
852 '''
854 azi = azimuth_numpy(lat0, lon0, lat, lon)
855 dist = distance_accurate50m_numpy(lat0, lon0, lat, lon)
856 n = num.cos(azi*d2r)*dist
857 e = num.sin(azi*d2r)*dist
858 return n, e
861_wgs84 = None
864def get_wgs84():
865 global _wgs84
866 if _wgs84 is None:
867 from geographiclib.geodesic import Geodesic
868 _wgs84 = Geodesic.WGS84
870 return _wgs84
873def amap(n):
874 def wrap(f):
875 if n == 1:
876 def func(*args):
877 it = num.nditer(args + (None,))
878 for ops in it:
879 ops[-1][...] = f(*ops[:-1])
881 return it.operands[-1]
882 elif n == 2:
883 def func(*args):
884 it = num.nditer(args + (None, None))
885 for ops in it:
886 ops[-2][...], ops[-1][...] = f(*ops[:-2])
888 return it.operands[-2], it.operands[-1]
890 return func
891 return wrap
894@amap(2)
895def ne_to_latlon2(lat0, lon0, north_m, east_m):
896 wgs84 = get_wgs84()
897 az = num.arctan2(east_m, north_m)*r2d
898 dist = num.sqrt(east_m**2 + north_m**2)
899 x = wgs84.Direct(lat0, lon0, az, dist)
900 return x['lat2'], x['lon2']
903@amap(2)
904def latlon_to_ne2(lat0, lon0, lat1, lon1):
905 wgs84 = get_wgs84()
906 x = wgs84.Inverse(lat0, lon0, lat1, lon1)
907 dist = x['s12']
908 az = x['azi1']
909 n = num.cos(az*d2r)*dist
910 e = num.sin(az*d2r)*dist
911 return n, e
914@amap(1)
915def distance_accurate15nm(lat1, lon1, lat2, lon2):
916 wgs84 = get_wgs84()
917 return wgs84.Inverse(lat1, lon1, lat2, lon2)['s12']
920def positive_region(region):
921 '''
922 Normalize parameterization of a rectangular geographical region.
924 :param region: ``(west, east, south, north)``
925 :returns: ``(west, east, south, north)``, where ``west <= east`` and
926 where ``west`` and ``east`` are in the range ``[-180., 180.+360.[``
927 '''
928 west, east, south, north = [float(x) for x in region]
930 assert -180. - 360. <= west < 180.
931 assert -180. < east <= 180. + 360.
932 assert -90. <= south < 90.
933 assert -90. < north <= 90.
935 if east < west:
936 east += 360.
938 if west < -180.:
939 west += 360.
940 east += 360.
942 return (west, east, south, north)
945def points_in_region(p, region):
946 '''
947 Check what points are contained in a rectangular geographical region.
949 :param p: NumPy array of shape ``(N, 2)`` where each row is a
950 ``(lat, lon)`` pair [deg]
951 :param region: ``(west, east, south, north)`` [deg]
952 :returns: NumPy array of shape ``(N)``, type ``bool``
953 '''
955 w, e, s, n = positive_region(region)
956 return num.logical_and(
957 num.logical_and(s <= p[:, 0], p[:, 0] <= n),
958 num.logical_or(
959 num.logical_and(w <= p[:, 1], p[:, 1] <= e),
960 num.logical_and(w-360. <= p[:, 1], p[:, 1] <= e-360.)))
963def point_in_region(p, region):
964 '''
965 Check if a point is contained in a rectangular geographical region.
967 :param p: ``(lat, lon)`` [deg]
968 :param region: ``(west, east, south, north)`` [deg]
969 :returns: ``bool``
970 '''
972 w, e, s, n = positive_region(region)
973 return num.logical_and(
974 num.logical_and(s <= p[0], p[0] <= n),
975 num.logical_or(
976 num.logical_and(w <= p[1], p[1] <= e),
977 num.logical_and(w-360. <= p[1], p[1] <= e-360.)))
980def radius_to_region(lat, lon, radius):
981 '''
982 Get a rectangular region which fully contains a given circular region.
984 :param lat,lon: center of circular region [deg]
985 :param radius: radius of circular region [m]
986 :return: rectangular region as ``(east, west, south, north)`` [deg]
987 '''
988 radius_deg = radius * m2d
989 if radius_deg < 45.:
990 lat_min = max(-90., lat - radius_deg)
991 lat_max = min(90., lat + radius_deg)
992 absmaxlat = max(abs(lat_min), abs(lat_max))
993 if absmaxlat > 89:
994 lon_min = -180.
995 lon_max = 180.
996 else:
997 lon_min = max(
998 -180. - 360.,
999 lon - radius_deg / math.cos(absmaxlat*d2r))
1000 lon_max = min(
1001 180. + 360.,
1002 lon + radius_deg / math.cos(absmaxlat*d2r))
1004 lon_min, lon_max, lat_min, lat_max = positive_region(
1005 (lon_min, lon_max, lat_min, lat_max))
1007 return lon_min, lon_max, lat_min, lat_max
1009 else:
1010 return None
1013def geographic_midpoint(lats, lons, weights=None):
1014 '''
1015 Calculate geographic midpoints by finding the center of gravity.
1017 This method suffers from instabilities if points are centered around the
1018 poles.
1020 :param lats: array of latitudes
1021 :param lons: array of longitudes
1022 :param weights: array weighting factors (optional)
1023 :type lats: :py:class:`numpy.ndarray`, ``(N)``
1024 :type lons: :py:class:`numpy.ndarray`, ``(N)``
1025 :type weights: :py:class:`numpy.ndarray`, ``(N)``
1027 :return: Latitudes and longitudes of the modpoints
1028 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
1029 '''
1030 if not weights:
1031 weights = num.ones(len(lats))
1033 total_weigth = num.sum(weights)
1034 weights /= total_weigth
1035 lats = lats * d2r
1036 lons = lons * d2r
1037 x = num.sum(num.cos(lats) * num.cos(lons) * weights)
1038 y = num.sum(num.cos(lats) * num.sin(lons) * weights)
1039 z = num.sum(num.sin(lats) * weights)
1041 lon = num.arctan2(y, x)
1042 hyp = num.sqrt(x**2 + y**2)
1043 lat = num.arctan2(z, hyp)
1045 return lat/d2r, lon/d2r
1048def geographic_midpoint_locations(locations, weights=None):
1049 coords = num.array([loc.effective_latlon
1050 for loc in locations])
1051 return geographic_midpoint(coords[:, 0], coords[:, 1], weights)
1054def geodetic_to_ecef(lat, lon, alt):
1055 '''
1056 Convert geodetic coordinates to Earth-Centered, Earth-Fixed (ECEF)
1057 Cartesian coordinates. [#1]_ [#2]_
1059 :param lat: Geodetic latitude in [deg].
1060 :param lon: Geodetic longitude in [deg].
1061 :param alt: Geodetic altitude (height) in [m] (positive for points outside
1062 the geoid).
1063 :type lat: float
1064 :type lon: float
1065 :type alt: float
1067 :return: ECEF Cartesian coordinates (X, Y, Z) in [m].
1068 :rtype: tuple, float
1070 .. [#1] https://en.wikipedia.org/wiki/ECEF
1071 .. [#2] https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1072 #From_geodetic_to_ECEF_coordinates
1073 '''
1075 f = earth_oblateness
1076 a = earthradius_equator
1077 e2 = 2*f - f**2
1079 lat, lon = num.radians(lat), num.radians(lon)
1080 # Normal (plumb line)
1081 N = a / num.sqrt(1.0 - (e2 * num.sin(lat)**2))
1083 X = (N+alt) * num.cos(lat) * num.cos(lon)
1084 Y = (N+alt) * num.cos(lat) * num.sin(lon)
1085 Z = (N*(1.0-e2) + alt) * num.sin(lat)
1087 return (X, Y, Z)
1090def ecef_to_geodetic(X, Y, Z):
1091 '''
1092 Convert Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates to
1093 geodetic coordinates (Ferrari's solution).
1095 :param X, Y, Z: Cartesian coordinates in ECEF system in [m].
1096 :type X, Y, Z: float
1098 :return: Geodetic coordinates (lat, lon, alt). Latitude and longitude are
1099 in [deg] and altitude is in [m]
1100 (positive for points outside the geoid).
1101 :rtype: tuple, float
1103 .. seealso ::
1104 https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1105 #The_application_of_Ferrari.27s_solution
1106 '''
1107 f = earth_oblateness
1108 a = earthradius_equator
1109 b = a * (1. - f)
1110 e2 = 2.*f - f**2
1112 # usefull
1113 a2 = a**2
1114 b2 = b**2
1115 e4 = e2**2
1116 X2 = X**2
1117 Y2 = Y**2
1118 Z2 = Z**2
1120 r = num.sqrt(X2 + Y2)
1121 r2 = r**2
1123 e_prime2 = (a2 - b2)/b2
1124 E2 = a2 - b2
1125 F = 54. * b2 * Z2
1126 G = r2 + (1.-e2)*Z2 - (e2*E2)
1127 C = (e4 * F * r2) / (G**3)
1128 S = cbrt(1. + C + num.sqrt(C**2 + 2.*C))
1129 P = F / (3. * (S + 1./S + 1.)**2 * G**2)
1130 Q = num.sqrt(1. + (2.*e4*P))
1132 dum1 = -(P*e2*r) / (1.+Q)
1133 dum2 = 0.5 * a2 * (1. + 1./Q)
1134 dum3 = (P * (1.-e2) * Z2) / (Q * (1.+Q))
1135 dum4 = 0.5 * P * r2
1136 r0 = dum1 + num.sqrt(dum2 - dum3 - dum4)
1138 U = num.sqrt((r - e2*r0)**2 + Z2)
1139 V = num.sqrt((r - e2*r0)**2 + (1.-e2)*Z2)
1140 Z0 = (b2*Z) / (a*V)
1142 alt = U * (1. - (b2 / (a*V)))
1143 lat = num.arctan((Z + e_prime2 * Z0)/r)
1144 lon = num.arctan2(Y, X)
1146 return (lat*r2d, lon*r2d, alt)
1149class Farside(Exception):
1150 pass
1153def latlon_to_xyz(latlons):
1154 if latlons.ndim == 1:
1155 return latlon_to_xyz(latlons[num.newaxis, :])[0]
1157 points = num.zeros((latlons.shape[0], 3))
1158 lats = latlons[:, 0]
1159 lons = latlons[:, 1]
1160 points[:, 0] = num.cos(lats*d2r) * num.cos(lons*d2r)
1161 points[:, 1] = num.cos(lats*d2r) * num.sin(lons*d2r)
1162 points[:, 2] = num.sin(lats*d2r)
1163 return points
1166def xyz_to_latlon(xyz):
1167 if xyz.ndim == 1:
1168 return xyz_to_latlon(xyz[num.newaxis, :])[0]
1170 latlons = num.zeros((xyz.shape[0], 2))
1171 latlons[:, 0] = num.arctan2(
1172 xyz[:, 2], num.sqrt(xyz[:, 0]**2 + xyz[:, 1]**2)) * r2d
1173 latlons[:, 1] = num.arctan2(
1174 xyz[:, 1], xyz[:, 0]) * r2d
1175 return latlons
1178def rot_to_00(lat, lon):
1179 rot0 = euler_to_matrix(0., -90.*d2r, 0.).A
1180 rot1 = euler_to_matrix(-d2r*lat, 0., -d2r*lon).A
1181 return num.dot(rot0.T, num.dot(rot1, rot0)).T
1184def distances3d(a, b):
1185 return num.sqrt(num.sum((a-b)**2, axis=a.ndim-1))
1188def circulation(points2):
1189 return num.sum(
1190 (points2[1:, 0] - points2[:-1, 0])
1191 * (points2[1:, 1] + points2[:-1, 1]))
1194def stereographic(points):
1195 dists = distances3d(points[1:, :], points[:-1, :])
1196 if dists.size > 0:
1197 maxdist = num.max(dists)
1198 cutoff = maxdist**2 / 2.
1199 else:
1200 cutoff = 1.0e-5
1202 points = points.copy()
1203 if num.any(points[:, 0] < -1. + cutoff):
1204 raise Farside()
1206 points_out = points[:, 1:].copy()
1207 factor = 1.0 / (1.0 + points[:, 0])
1208 points_out *= factor[:, num.newaxis]
1210 return points_out
1213def stereographic_poly(points):
1214 dists = distances3d(points[1:, :], points[:-1, :])
1215 if dists.size > 0:
1216 maxdist = num.max(dists)
1217 cutoff = maxdist**2 / 2.
1218 else:
1219 cutoff = 1.0e-5
1221 points = points.copy()
1222 if num.any(points[:, 0] < -1. + cutoff):
1223 raise Farside()
1225 points_out = points[:, 1:].copy()
1226 factor = 1.0 / (1.0 + points[:, 0])
1227 points_out *= factor[:, num.newaxis]
1229 if circulation(points_out) >= 0:
1230 raise Farside()
1232 return points_out
1235def gnomonic_x(points, cutoff=0.01):
1236 points_out = points[:, 1:].copy()
1237 if num.any(points[:, 0] < cutoff):
1238 raise Farside()
1240 factor = 1.0 / points[:, 0]
1241 points_out *= factor[:, num.newaxis]
1242 return points_out
1245def cneg(i, x):
1246 if i == 1:
1247 return x
1248 else:
1249 return num.logical_not(x)
1252def contains_points(polygon, points):
1253 '''
1254 Test which points are inside polygon on a sphere.
1256 :param polygon: Point coordinates defining the polygon [deg].
1257 :type polygon: :py:class:`numpy.ndarray` of shape (N, 2), second index
1258 0=lat, 1=lon
1259 :param points: Coordinates of points to test [deg].
1260 :type points: :py:class:`numpy.ndarray` of shape (N, 2), second index
1261 0=lat, 1=lon
1263 :returns: Boolean mask array.
1264 :rtype: :py:class:`numpy.ndarray` of shape (N,)
1266 The inside of the polygon is defined as the area which is to the left hand
1267 side of an observer walking the polygon line, points in order, on the
1268 sphere. Lines between the polygon points are treated as great circle paths.
1269 The polygon may be arbitrarily complex, as long as it does not have any
1270 crossings or thin parts with zero width. The polygon may contain the poles
1271 and is allowed to wrap around the sphere multiple times.
1273 The algorithm works by consecutive cutting of the polygon into (almost)
1274 hemispheres and subsequent Gnomonic projections to perform the
1275 point-in-polygon tests on a 2D plane.
1276 '''
1278 and_ = num.logical_and
1280 points_xyz = latlon_to_xyz(points)
1281 mask_x = 0. <= points_xyz[:, 0]
1282 mask_y = 0. <= points_xyz[:, 1]
1283 mask_z = 0. <= points_xyz[:, 2]
1285 result = num.zeros(points.shape[0], dtype=int)
1287 for ix in [-1, 1]:
1288 for iy in [-1, 1]:
1289 for iz in [-1, 1]:
1290 mask = and_(
1291 and_(cneg(ix, mask_x), cneg(iy, mask_y)),
1292 cneg(iz, mask_z))
1294 center_xyz = num.array([ix, iy, iz], dtype=float)
1296 lat, lon = xyz_to_latlon(center_xyz)
1297 rot = rot_to_00(lat, lon)
1299 points_rot_xyz = num.dot(rot, points_xyz[mask, :].T).T
1300 points_rot_pro = gnomonic_x(points_rot_xyz)
1302 offset = 0.01
1304 poly_xyz = latlon_to_xyz(polygon)
1305 poly_rot_xyz = num.dot(rot, poly_xyz.T).T
1306 poly_rot_xyz[:, 0] -= offset
1307 groups = spoly_cut([poly_rot_xyz], axis=0)
1309 for poly_rot_group_xyz in groups[1]:
1310 poly_rot_group_xyz[:, 0] += offset
1312 poly_rot_group_pro = gnomonic_x(
1313 poly_rot_group_xyz)
1315 if circulation(poly_rot_group_pro) > 0:
1316 result[mask] += path_contains_points(
1317 poly_rot_group_pro, points_rot_pro)
1318 else:
1319 result[mask] -= path_contains_points(
1320 poly_rot_group_pro, points_rot_pro)
1322 return result.astype(num.bool)
1325def contains_point(polygon, point):
1326 '''
1327 Test if point is inside polygon on a sphere.
1329 :param polygon: Point coordinates defining the polygon [deg].
1330 :type polygon: :py:class:`numpy.ndarray` of shape (N, 2), second index
1331 0=lat, 1=lon
1332 :param point: Coordinates ``(lat, lon)`` of point to test [deg].
1334 Convenience wrapper to :py:func:`contains_points` to test a single point.
1335 '''
1337 return bool(
1338 contains_points(polygon, num.asarray(point)[num.newaxis, :])[0])