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 Absolute latitudes and longitudes are calculated from relative changes.
604 Convenience wrapper to :py:func:`azidist_to_latlon_rad` with azimuth and
605 distance given in degrees.
607 :param lat0: Latitude origin of the cartesian coordinate system.
608 :type lat0: float
609 :param lon0: Longitude origin of the cartesian coordinate system.
610 :type lon0: float
611 :param azimuth_deg: Azimuth from origin in degrees.
612 :type azimuth_deg: :py:class:`numpy.ndarray`, ``(N)``
613 :param distance_deg: Distances from origin in degrees.
614 :type distance_deg: :py:class:`numpy.ndarray`, ``(N)``
616 :return: Array with latitudes and longitudes
617 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
618 '''
620 return azidist_to_latlon_rad(
621 lat0, lon0, azimuth_deg/180.*num.pi, distance_deg/180.*num.pi)
624def azidist_to_latlon_rad(lat0, lon0, azimuth_rad, distance_rad):
625 '''
626 Absolute latitudes and longitudes are calculated from relative changes.
628 For numerical stability a range between of ``-1.0`` and ``1.0`` is
629 enforced for ``c`` and ``alpha``.
631 .. math::
633 \\Delta {\\bf a}_i \\; \\text{and} \\; \\Delta \\gamma_i \\;
634 \\text{are relative distances and azimuths from lat0 and lon0 for
635 \\textit{i} source points of a finite source.}
637 .. math::
639 \\mathrm{b} &= \\frac{\\pi}{2} -\\frac{\\pi}{180}\\;\\mathrm{lat_0}\\\\
640 {\\bf c}_i &=\\arccos[\\; \\cos(\\Delta {\\bf{a}}_i)
641 \\cos(\\mathrm{b}) + |\\Delta \\gamma_i| \\,
642 \\sin(\\Delta {\\bf a}_i)
643 \\sin(\\mathrm{b})\\; ] \\\\
644 \\mathrm{lat}_i &= \\frac{180}{\\pi}
645 \\left(\\frac{\\pi}{2} - {\\bf c}_i \\right)
647 .. math::
649 \\alpha_i &= \\arcsin \\left[ \\; \\frac{ \\sin(\\Delta {\\bf a}_i )
650 \\sin(|\\Delta \\gamma_i|)}{\\sin({\\bf c}_i)}\\;
651 \\right] \\\\
652 \\alpha_i &= \\begin{cases}
653 \\alpha_i, &\\text{if} \\; \\cos(\\Delta {\\bf a}_i) -
654 \\cos(\\mathrm{b}) \\cos({\\bf{c}}_i) > 0, \\;
655 \\text{else} \\\\
656 \\pi - \\alpha_i, & \\text{if} \\; \\alpha_i > 0,\\;
657 \\text{else}\\\\
658 -\\pi - \\alpha_i, & \\text{if} \\; \\alpha_i < 0.
659 \\end{cases} \\\\
660 \\mathrm{lon}_i &= \\mathrm{lon_0} +
661 \\frac{180}{\\pi} \\,
662 \\frac{\\Delta \\gamma_i }{|\\Delta \\gamma_i|}
663 \\cdot \\alpha_i
664 \\text{, with $\\alpha_i \\in [-\\pi,\\pi]$}
666 :param lat0: Latitude origin of the cartesian coordinate system.
667 :param lon0: Longitude origin of the cartesian coordinate system.
668 :param distance_rad: Distances from origin in radians.
669 :param azimuth_rad: Azimuth from origin radians.
670 :type distance_rad: :py:class:`numpy.ndarray`, ``(N)``
671 :type azimuth_rad: :py:class:`numpy.ndarray`, ``(N)``
672 :type lat0: float
673 :type lon0: float
675 :return: Array with latitudes and longitudes
676 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
677 '''
679 a = distance_rad
680 gamma = azimuth_rad
682 b = math.pi/2.-lat0*d2r
684 alphasign = 1.
685 alphasign = num.where(gamma < 0, -1., 1.)
686 gamma = num.abs(gamma)
688 c = num.arccos(clip(
689 num.cos(a)*num.cos(b)+num.sin(a)*num.sin(b)*num.cos(gamma), -1., 1.))
691 alpha = num.arcsin(clip(
692 num.sin(a)*num.sin(gamma)/num.sin(c), -1., 1.))
694 alpha = num.where(
695 num.cos(a)-num.cos(b)*num.cos(c) < 0,
696 num.where(alpha > 0, math.pi-alpha, -math.pi-alpha),
697 alpha)
699 lat = r2d * (math.pi/2. - c)
700 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
702 return lat, lon
705def ne_to_latlon_alternative_method(lat0, lon0, north_m, east_m):
706 '''
707 Transform local cartesian coordinates to latitude and longitude.
709 Like :py:func:`pyrocko.orthodrome.ne_to_latlon`,
710 but this method (implementation below), although it should be numerically
711 more stable, suffers problems at points which are *across the pole*
712 as seen from the cartesian origin.
714 .. math::
716 \\text{distance change:}\\; \\Delta {{\\bf a}_i} &=
717 \\sqrt{{\\bf{y}}^2_i + {\\bf{x}}^2_i }/ \\mathrm{R_E},\\\\
718 \\text{azimuth change:}\\; \\Delta {\\bf \\gamma}_i &=
719 \\arctan( {\\bf x}_i {\\bf y}_i). \\\\
720 \\mathrm{b} &=
721 \\frac{\\pi}{2} -\\frac{\\pi}{180} \\;\\mathrm{lat_0}\\\\
723 .. math::
725 {{\\bf z}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i -
726 \\mathrm{b}}{2} \\right)}
727 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
728 {{\\bf n}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i +
729 \\mathrm{b}}{2} \\right)}
730 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
731 {{\\bf z}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i -
732 \\mathrm{b}}{2} \\right)}
733 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
734 {{\\bf n}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i +
735 \\mathrm{b}}{2} \\right)}
736 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
737 {{\\bf t}_1}_i &= \\arctan{\\left( \\frac{{{\\bf z}_1}_i}
738 {{{\\bf n}_1}_i} \\right) }\\\\
739 {{\\bf t}_2}_i &= \\arctan{\\left( \\frac{{{\\bf z}_2}_i}
740 {{{\\bf n}_2}_i} \\right) } \\\\[0.5cm]
741 c &= \\begin{cases}
742 2 \\cdot \\arccos \\left( {{\\bf z}_1}_i / \\sin({{\\bf t}_1}_i)
743 \\right),\\; \\text{if }
744 |\\sin({{\\bf t}_1}_i)| >
745 |\\sin({{\\bf t}_2}_i)|,\\; \\text{else} \\\\
746 2 \\cdot \\arcsin{\\left( {{\\bf z}_2}_i /
747 \\sin({{\\bf t}_2}_i) \\right)}.
748 \\end{cases}\\\\
750 .. math::
752 {\\bf {lat}}_i &= \\frac{180}{ \\pi } \\left( \\frac{\\pi}{2}
753 - {\\bf {c}}_i \\right) \\\\
754 {\\bf {lon}}_i &= {\\bf {lon}}_0 + \\frac{180}{ \\pi }
755 \\frac{\\gamma_i}{|\\gamma_i|},
756 \\text{ with}\\; \\gamma \\in [-\\pi,\\pi]
758 :param lat0: Latitude origin of the cartesian coordinate system.
759 :param lon0: Longitude origin of the cartesian coordinate system.
760 :param north_m: Northing distances from origin in meters.
761 :param east_m: Easting distances from origin in meters.
762 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
763 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
764 :type lat0: float
765 :type lon0: float
767 :return: Array with latitudes and longitudes
768 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
769 '''
771 b = math.pi/2.-lat0*d2r
772 a = num.sqrt(north_m**2+east_m**2)/earthradius
774 gamma = num.arctan2(east_m, north_m)
775 alphasign = 1.
776 alphasign = num.where(gamma < 0., -1., 1.)
777 gamma = num.abs(gamma)
779 z1 = num.cos((a-b)/2.)*num.cos(gamma/2.)
780 n1 = num.cos((a+b)/2.)*num.sin(gamma/2.)
781 z2 = num.sin((a-b)/2.)*num.cos(gamma/2.)
782 n2 = num.sin((a+b)/2.)*num.sin(gamma/2.)
783 t1 = num.arctan2(z1, n1)
784 t2 = num.arctan2(z2, n2)
786 alpha = t1 + t2
788 sin_t1 = num.sin(t1)
789 sin_t2 = num.sin(t2)
790 c = num.where(
791 num.abs(sin_t1) > num.abs(sin_t2),
792 num.arccos(z1/sin_t1)*2.,
793 num.arcsin(z2/sin_t2)*2.)
795 lat = r2d * (math.pi/2. - c)
796 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
797 return lat, lon
800def latlon_to_ne(*args):
801 '''
802 Relative cartesian coordinates with respect to a reference location.
804 For two locations, a reference location A and another location B, given in
805 geographical coordinates in degrees, the corresponding cartesian
806 coordinates are calculated.
807 Assisting functions are :py:func:`pyrocko.orthodrome.azimuth` and
808 :py:func:`pyrocko.orthodrome.distance_accurate50m`.
810 .. math::
812 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
813 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
814 \\mathrm{)}\\\\[0.3cm]
816 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180}
817 \\varphi_{\\mathrm{azi},AB} )\\\\
818 e &= D_{AB} \\cdot
819 \\sin( \\frac{\\pi }{180} \\varphi_{\\mathrm{azi},AB})
821 :param refloc: Location reference point
822 :type refloc: :py:class:`pyrocko.orthodrome.Loc`
823 :param loc: Location of interest
824 :type loc: :py:class:`pyrocko.orthodrome.Loc`
826 :return: Northing and easting from refloc to location
827 :rtype: tuple, float
829 '''
831 azi = azimuth(*args)
832 dist = distance_accurate50m(*args)
833 n, e = math.cos(azi*d2r)*dist, math.sin(azi*d2r)*dist
834 return n, e
837def latlon_to_ne_numpy(lat0, lon0, lat, lon):
838 '''
839 Relative cartesian coordinates with respect to a reference location.
841 For two locations, a reference location (``lat0``, ``lon0``) and another
842 location B, given in geographical coordinates in degrees,
843 the corresponding cartesian coordinates are calculated.
844 Assisting functions are :py:func:`azimuth`
845 and :py:func:`distance_accurate50m`.
847 :param lat0: reference location latitude
848 :param lon0: reference location longitude
849 :param lat: absolute location latitude
850 :param lon: absolute location longitude
852 :return: ``(n, e)``: relative north and east positions
853 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
855 Implemented formulations:
857 .. math::
859 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
860 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
861 \\mathrm{)}\\\\[0.3cm]
863 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180} \\varphi_{
864 \\mathrm{azi},AB} )\\\\
865 e &= D_{AB} \\cdot \\sin( \\frac{\\pi }{180} \\varphi_{
866 \\mathrm{azi},AB} )
867 '''
869 azi = azimuth_numpy(lat0, lon0, lat, lon)
870 dist = distance_accurate50m_numpy(lat0, lon0, lat, lon)
871 n = num.cos(azi*d2r)*dist
872 e = num.sin(azi*d2r)*dist
873 return n, e
876_wgs84 = None
879def get_wgs84():
880 global _wgs84
881 if _wgs84 is None:
882 from geographiclib.geodesic import Geodesic
883 _wgs84 = Geodesic.WGS84
885 return _wgs84
888def amap(n):
889 def wrap(f):
890 if n == 1:
891 def func(*args):
892 it = num.nditer(args + (None,))
893 for ops in it:
894 ops[-1][...] = f(*ops[:-1])
896 return it.operands[-1]
897 elif n == 2:
898 def func(*args):
899 it = num.nditer(args + (None, None))
900 for ops in it:
901 ops[-2][...], ops[-1][...] = f(*ops[:-2])
903 return it.operands[-2], it.operands[-1]
905 return func
906 return wrap
909@amap(2)
910def ne_to_latlon2(lat0, lon0, north_m, east_m):
911 wgs84 = get_wgs84()
912 az = num.arctan2(east_m, north_m)*r2d
913 dist = num.sqrt(east_m**2 + north_m**2)
914 x = wgs84.Direct(lat0, lon0, az, dist)
915 return x['lat2'], x['lon2']
918@amap(2)
919def latlon_to_ne2(lat0, lon0, lat1, lon1):
920 wgs84 = get_wgs84()
921 x = wgs84.Inverse(lat0, lon0, lat1, lon1)
922 dist = x['s12']
923 az = x['azi1']
924 n = num.cos(az*d2r)*dist
925 e = num.sin(az*d2r)*dist
926 return n, e
929@amap(1)
930def distance_accurate15nm(lat1, lon1, lat2, lon2):
931 wgs84 = get_wgs84()
932 return wgs84.Inverse(lat1, lon1, lat2, lon2)['s12']
935def positive_region(region):
936 '''
937 Normalize parameterization of a rectangular geographical region.
939 :param region: ``(west, east, south, north)``
940 :returns: ``(west, east, south, north)``, where ``west <= east`` and
941 where ``west`` and ``east`` are in the range ``[-180., 180.+360.[``
942 '''
943 west, east, south, north = [float(x) for x in region]
945 assert -180. - 360. <= west < 180.
946 assert -180. < east <= 180. + 360.
947 assert -90. <= south < 90.
948 assert -90. < north <= 90.
950 if east < west:
951 east += 360.
953 if west < -180.:
954 west += 360.
955 east += 360.
957 return (west, east, south, north)
960def points_in_region(p, region):
961 '''
962 Check what points are contained in a rectangular geographical region.
964 :param p: NumPy array of shape ``(N, 2)`` where each row is a
965 ``(lat, lon)`` pair [deg]
966 :param region: ``(west, east, south, north)`` [deg]
967 :returns: NumPy array of shape ``(N)``, type ``bool``
968 '''
970 w, e, s, n = positive_region(region)
971 return num.logical_and(
972 num.logical_and(s <= p[:, 0], p[:, 0] <= n),
973 num.logical_or(
974 num.logical_and(w <= p[:, 1], p[:, 1] <= e),
975 num.logical_and(w-360. <= p[:, 1], p[:, 1] <= e-360.)))
978def point_in_region(p, region):
979 '''
980 Check if a point is contained in a rectangular geographical region.
982 :param p: ``(lat, lon)`` [deg]
983 :param region: ``(west, east, south, north)`` [deg]
984 :returns: ``bool``
985 '''
987 w, e, s, n = positive_region(region)
988 return num.logical_and(
989 num.logical_and(s <= p[0], p[0] <= n),
990 num.logical_or(
991 num.logical_and(w <= p[1], p[1] <= e),
992 num.logical_and(w-360. <= p[1], p[1] <= e-360.)))
995def radius_to_region(lat, lon, radius):
996 '''
997 Get a rectangular region which fully contains a given circular region.
999 :param lat,lon: center of circular region [deg]
1000 :param radius: radius of circular region [m]
1001 :return: rectangular region as ``(east, west, south, north)`` [deg]
1002 '''
1003 radius_deg = radius * m2d
1004 if radius_deg < 45.:
1005 lat_min = max(-90., lat - radius_deg)
1006 lat_max = min(90., lat + radius_deg)
1007 absmaxlat = max(abs(lat_min), abs(lat_max))
1008 if absmaxlat > 89:
1009 lon_min = -180.
1010 lon_max = 180.
1011 else:
1012 lon_min = max(
1013 -180. - 360.,
1014 lon - radius_deg / math.cos(absmaxlat*d2r))
1015 lon_max = min(
1016 180. + 360.,
1017 lon + radius_deg / math.cos(absmaxlat*d2r))
1019 lon_min, lon_max, lat_min, lat_max = positive_region(
1020 (lon_min, lon_max, lat_min, lat_max))
1022 return lon_min, lon_max, lat_min, lat_max
1024 else:
1025 return None
1028def geographic_midpoint(lats, lons, weights=None):
1029 '''
1030 Calculate geographic midpoints by finding the center of gravity.
1032 This method suffers from instabilities if points are centered around the
1033 poles.
1035 :param lats: array of latitudes
1036 :param lons: array of longitudes
1037 :param weights: array weighting factors (optional)
1038 :type lats: :py:class:`numpy.ndarray`, ``(N)``
1039 :type lons: :py:class:`numpy.ndarray`, ``(N)``
1040 :type weights: :py:class:`numpy.ndarray`, ``(N)``
1042 :return: Latitudes and longitudes of the modpoints
1043 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
1044 '''
1045 if not weights:
1046 weights = num.ones(len(lats))
1048 total_weigth = num.sum(weights)
1049 weights /= total_weigth
1050 lats = lats * d2r
1051 lons = lons * d2r
1052 x = num.sum(num.cos(lats) * num.cos(lons) * weights)
1053 y = num.sum(num.cos(lats) * num.sin(lons) * weights)
1054 z = num.sum(num.sin(lats) * weights)
1056 lon = num.arctan2(y, x)
1057 hyp = num.sqrt(x**2 + y**2)
1058 lat = num.arctan2(z, hyp)
1060 return lat/d2r, lon/d2r
1063def geographic_midpoint_locations(locations, weights=None):
1064 coords = num.array([loc.effective_latlon
1065 for loc in locations])
1066 return geographic_midpoint(coords[:, 0], coords[:, 1], weights)
1069def geodetic_to_ecef(lat, lon, alt):
1070 '''
1071 Convert geodetic coordinates to Earth-Centered, Earth-Fixed (ECEF)
1072 Cartesian coordinates. [#1]_ [#2]_
1074 :param lat: Geodetic latitude in [deg].
1075 :param lon: Geodetic longitude in [deg].
1076 :param alt: Geodetic altitude (height) in [m] (positive for points outside
1077 the geoid).
1078 :type lat: float
1079 :type lon: float
1080 :type alt: float
1082 :return: ECEF Cartesian coordinates (X, Y, Z) in [m].
1083 :rtype: tuple, float
1085 .. [#1] https://en.wikipedia.org/wiki/ECEF
1086 .. [#2] https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1087 #From_geodetic_to_ECEF_coordinates
1088 '''
1090 f = earth_oblateness
1091 a = earthradius_equator
1092 e2 = 2*f - f**2
1094 lat, lon = num.radians(lat), num.radians(lon)
1095 # Normal (plumb line)
1096 N = a / num.sqrt(1.0 - (e2 * num.sin(lat)**2))
1098 X = (N+alt) * num.cos(lat) * num.cos(lon)
1099 Y = (N+alt) * num.cos(lat) * num.sin(lon)
1100 Z = (N*(1.0-e2) + alt) * num.sin(lat)
1102 return (X, Y, Z)
1105def ecef_to_geodetic(X, Y, Z):
1106 '''
1107 Convert Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates to
1108 geodetic coordinates (Ferrari's solution).
1110 :param X, Y, Z: Cartesian coordinates in ECEF system in [m].
1111 :type X, Y, Z: float
1113 :return: Geodetic coordinates (lat, lon, alt). Latitude and longitude are
1114 in [deg] and altitude is in [m]
1115 (positive for points outside the geoid).
1116 :rtype: tuple, float
1118 .. seealso ::
1119 https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1120 #The_application_of_Ferrari.27s_solution
1121 '''
1122 f = earth_oblateness
1123 a = earthradius_equator
1124 b = a * (1. - f)
1125 e2 = 2.*f - f**2
1127 # usefull
1128 a2 = a**2
1129 b2 = b**2
1130 e4 = e2**2
1131 X2 = X**2
1132 Y2 = Y**2
1133 Z2 = Z**2
1135 r = num.sqrt(X2 + Y2)
1136 r2 = r**2
1138 e_prime2 = (a2 - b2)/b2
1139 E2 = a2 - b2
1140 F = 54. * b2 * Z2
1141 G = r2 + (1.-e2)*Z2 - (e2*E2)
1142 C = (e4 * F * r2) / (G**3)
1143 S = cbrt(1. + C + num.sqrt(C**2 + 2.*C))
1144 P = F / (3. * (S + 1./S + 1.)**2 * G**2)
1145 Q = num.sqrt(1. + (2.*e4*P))
1147 dum1 = -(P*e2*r) / (1.+Q)
1148 dum2 = 0.5 * a2 * (1. + 1./Q)
1149 dum3 = (P * (1.-e2) * Z2) / (Q * (1.+Q))
1150 dum4 = 0.5 * P * r2
1151 r0 = dum1 + num.sqrt(dum2 - dum3 - dum4)
1153 U = num.sqrt((r - e2*r0)**2 + Z2)
1154 V = num.sqrt((r - e2*r0)**2 + (1.-e2)*Z2)
1155 Z0 = (b2*Z) / (a*V)
1157 alt = U * (1. - (b2 / (a*V)))
1158 lat = num.arctan((Z + e_prime2 * Z0)/r)
1159 lon = num.arctan2(Y, X)
1161 return (lat*r2d, lon*r2d, alt)
1164class Farside(Exception):
1165 pass
1168def latlon_to_xyz(latlons):
1169 if latlons.ndim == 1:
1170 return latlon_to_xyz(latlons[num.newaxis, :])[0]
1172 points = num.zeros((latlons.shape[0], 3))
1173 lats = latlons[:, 0]
1174 lons = latlons[:, 1]
1175 points[:, 0] = num.cos(lats*d2r) * num.cos(lons*d2r)
1176 points[:, 1] = num.cos(lats*d2r) * num.sin(lons*d2r)
1177 points[:, 2] = num.sin(lats*d2r)
1178 return points
1181def xyz_to_latlon(xyz):
1182 if xyz.ndim == 1:
1183 return xyz_to_latlon(xyz[num.newaxis, :])[0]
1185 latlons = num.zeros((xyz.shape[0], 2))
1186 latlons[:, 0] = num.arctan2(
1187 xyz[:, 2], num.sqrt(xyz[:, 0]**2 + xyz[:, 1]**2)) * r2d
1188 latlons[:, 1] = num.arctan2(
1189 xyz[:, 1], xyz[:, 0]) * r2d
1190 return latlons
1193def rot_to_00(lat, lon):
1194 rot0 = euler_to_matrix(0., -90.*d2r, 0.).A
1195 rot1 = euler_to_matrix(-d2r*lat, 0., -d2r*lon).A
1196 return num.dot(rot0.T, num.dot(rot1, rot0)).T
1199def distances3d(a, b):
1200 return num.sqrt(num.sum((a-b)**2, axis=a.ndim-1))
1203def circulation(points2):
1204 return num.sum(
1205 (points2[1:, 0] - points2[:-1, 0])
1206 * (points2[1:, 1] + points2[:-1, 1]))
1209def stereographic(points):
1210 dists = distances3d(points[1:, :], points[:-1, :])
1211 if dists.size > 0:
1212 maxdist = num.max(dists)
1213 cutoff = maxdist**2 / 2.
1214 else:
1215 cutoff = 1.0e-5
1217 points = points.copy()
1218 if num.any(points[:, 0] < -1. + cutoff):
1219 raise Farside()
1221 points_out = points[:, 1:].copy()
1222 factor = 1.0 / (1.0 + points[:, 0])
1223 points_out *= factor[:, num.newaxis]
1225 return points_out
1228def stereographic_poly(points):
1229 dists = distances3d(points[1:, :], points[:-1, :])
1230 if dists.size > 0:
1231 maxdist = num.max(dists)
1232 cutoff = maxdist**2 / 2.
1233 else:
1234 cutoff = 1.0e-5
1236 points = points.copy()
1237 if num.any(points[:, 0] < -1. + cutoff):
1238 raise Farside()
1240 points_out = points[:, 1:].copy()
1241 factor = 1.0 / (1.0 + points[:, 0])
1242 points_out *= factor[:, num.newaxis]
1244 if circulation(points_out) >= 0:
1245 raise Farside()
1247 return points_out
1250def gnomonic_x(points, cutoff=0.01):
1251 points_out = points[:, 1:].copy()
1252 if num.any(points[:, 0] < cutoff):
1253 raise Farside()
1255 factor = 1.0 / points[:, 0]
1256 points_out *= factor[:, num.newaxis]
1257 return points_out
1260def cneg(i, x):
1261 if i == 1:
1262 return x
1263 else:
1264 return num.logical_not(x)
1267def contains_points(polygon, points):
1268 '''
1269 Test which points are inside polygon on a sphere.
1271 :param polygon: Point coordinates defining the polygon [deg].
1272 :type polygon: :py:class:`numpy.ndarray` of shape (N, 2), second index
1273 0=lat, 1=lon
1274 :param points: Coordinates of points to test [deg].
1275 :type points: :py:class:`numpy.ndarray` of shape (N, 2), second index
1276 0=lat, 1=lon
1278 :returns: Boolean mask array.
1279 :rtype: :py:class:`numpy.ndarray` of shape (N,)
1281 The inside of the polygon is defined as the area which is to the left hand
1282 side of an observer walking the polygon line, points in order, on the
1283 sphere. Lines between the polygon points are treated as great circle paths.
1284 The polygon may be arbitrarily complex, as long as it does not have any
1285 crossings or thin parts with zero width. The polygon may contain the poles
1286 and is allowed to wrap around the sphere multiple times.
1288 The algorithm works by consecutive cutting of the polygon into (almost)
1289 hemispheres and subsequent Gnomonic projections to perform the
1290 point-in-polygon tests on a 2D plane.
1291 '''
1293 and_ = num.logical_and
1295 points_xyz = latlon_to_xyz(points)
1296 mask_x = 0. <= points_xyz[:, 0]
1297 mask_y = 0. <= points_xyz[:, 1]
1298 mask_z = 0. <= points_xyz[:, 2]
1300 result = num.zeros(points.shape[0], dtype=int)
1302 for ix in [-1, 1]:
1303 for iy in [-1, 1]:
1304 for iz in [-1, 1]:
1305 mask = and_(
1306 and_(cneg(ix, mask_x), cneg(iy, mask_y)),
1307 cneg(iz, mask_z))
1309 center_xyz = num.array([ix, iy, iz], dtype=float)
1311 lat, lon = xyz_to_latlon(center_xyz)
1312 rot = rot_to_00(lat, lon)
1314 points_rot_xyz = num.dot(rot, points_xyz[mask, :].T).T
1315 points_rot_pro = gnomonic_x(points_rot_xyz)
1317 offset = 0.01
1319 poly_xyz = latlon_to_xyz(polygon)
1320 poly_rot_xyz = num.dot(rot, poly_xyz.T).T
1321 poly_rot_xyz[:, 0] -= offset
1322 groups = spoly_cut([poly_rot_xyz], axis=0)
1324 for poly_rot_group_xyz in groups[1]:
1325 poly_rot_group_xyz[:, 0] += offset
1327 poly_rot_group_pro = gnomonic_x(
1328 poly_rot_group_xyz)
1330 if circulation(poly_rot_group_pro) > 0:
1331 result[mask] += path_contains_points(
1332 poly_rot_group_pro, points_rot_pro)
1333 else:
1334 result[mask] -= path_contains_points(
1335 poly_rot_group_pro, points_rot_pro)
1337 return result.astype(num.bool)
1340def contains_point(polygon, point):
1341 '''
1342 Test if point is inside polygon on a sphere.
1344 :param polygon: Point coordinates defining the polygon [deg].
1345 :type polygon: :py:class:`numpy.ndarray` of shape (N, 2), second index
1346 0=lat, 1=lon
1347 :param point: Coordinates ``(lat, lon)`` of point to test [deg].
1349 Convenience wrapper to :py:func:`contains_points` to test a single point.
1350 '''
1352 return bool(
1353 contains_points(polygon, num.asarray(point)[num.newaxis, :])[0])