1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import math
7import numpy as num
9from .moment_tensor import euler_to_matrix
10from .config import config
11from .plot.beachball import spoly_cut
13from matplotlib.path import Path
15d2r = math.pi/180.
16r2d = 1./d2r
17earth_oblateness = 1./298.257223563
18earthradius_equator = 6378.14 * 1000.
19earthradius = config().earthradius
20d2m = earthradius_equator*math.pi/180.
21m2d = 1./d2m
23_testpath = Path([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)], closed=True)
25raise_if_slow_path_contains_points = False
28class Slow(Exception):
29 pass
32if hasattr(_testpath, 'contains_points') and num.all(
33 _testpath.contains_points([(0.5, 0.5), (1.5, 0.5)]) == [True, False]):
35 def path_contains_points(verts, points):
36 p = Path(verts, closed=True)
37 return p.contains_points(points).astype(bool)
39else:
40 # work around missing contains_points and bug in matplotlib ~ v1.2.0
42 def path_contains_points(verts, points):
43 if raise_if_slow_path_contains_points:
44 # used by unit test to skip slow gshhg_example.py
45 raise Slow()
47 p = Path(verts, closed=True)
48 result = num.zeros(points.shape[0], dtype=bool)
49 for i in range(result.size):
50 result[i] = p.contains_point(points[i, :])
52 return result
55try:
56 cbrt = num.cbrt
57except AttributeError:
58 def cbrt(x):
59 return x**(1./3.)
62def float_array_broadcast(*args):
63 return num.broadcast_arrays(*[
64 num.asarray(x, dtype=float) for x in args])
67class Loc(object):
68 '''
69 Simple location representation.
71 :attrib lat: Latitude in [deg].
72 :attrib lon: Longitude in [deg].
73 '''
74 def __init__(self, lat, lon):
75 self.lat = lat
76 self.lon = lon
79def clip(x, mi, ma):
80 '''
81 Clip values of an array.
83 :param x: Continunous data to be clipped.
84 :param mi: Clip minimum.
85 :param ma: Clip maximum.
86 :type x: :py:class:`numpy.ndarray`
87 :type mi: float
88 :type ma: float
90 :return: Clipped data.
91 :rtype: :py:class:`numpy.ndarray`
92 '''
93 return num.minimum(num.maximum(mi, x), ma)
96def wrap(x, mi, ma):
97 '''
98 Wrapping continuous data to fundamental phase values.
100 .. math::
101 x_{\\mathrm{wrapped}} = x_{\\mathrm{cont},i} -
102 \\frac{ x_{\\mathrm{cont},i} - r_{\\mathrm{min}} }
103 { r_{\\mathrm{max}} - r_{\\mathrm{min}}}
104 \\cdot ( r_{\\mathrm{max}} - r_{\\mathrm{min}}),\\quad
105 x_{\\mathrm{wrapped}}\\; \\in
106 \\;[ r_{\\mathrm{min}},\\, r_{\\mathrm{max}}].
108 :param x: Continunous data to be wrapped.
109 :param mi: Minimum value of wrapped data.
110 :param ma: Maximum value of wrapped data.
111 :type x: :py:class:`numpy.ndarray`
112 :type mi: float
113 :type ma: float
115 :return: Wrapped data.
116 :rtype: :py:class:`numpy.ndarray`
117 '''
118 return x - num.floor((x-mi)/(ma-mi)) * (ma-mi)
121def _latlon_pair(args):
122 if len(args) == 2:
123 a, b = args
124 return a.lat, a.lon, b.lat, b.lon
126 elif len(args) == 4:
127 return args
130def cosdelta(*args):
131 '''
132 Cosine of the angular distance between two points ``a`` and ``b`` on a
133 sphere.
135 This function (find implementation below) returns the cosine of the
136 distance angle 'delta' between two points ``a`` and ``b``, coordinates of
137 which are expected to be given in geographical coordinates and in degrees.
138 For numerical stability a maximum of 1.0 is enforced.
140 .. math::
142 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
143 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
144 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
145 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\[0.5cm]
147 \\cos(\\Delta) = \\min( 1.0, \\quad \\sin( A_{\\mathrm{lat'}})
148 \\sin( B_{\\mathrm{lat'}} ) +
149 \\cos(A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
150 \\cos( B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )
152 :param a: Location point A.
153 :type a: :py:class:`pyrocko.orthodrome.Loc`
154 :param b: Location point B.
155 :type b: :py:class:`pyrocko.orthodrome.Loc`
157 :return: Cosdelta.
158 :rtype: float
159 '''
161 alat, alon, blat, blon = _latlon_pair(args)
163 return min(
164 1.0,
165 math.sin(alat*d2r) * math.sin(blat*d2r) +
166 math.cos(alat*d2r) * math.cos(blat*d2r) *
167 math.cos(d2r*(blon-alon)))
170def cosdelta_numpy(a_lats, a_lons, b_lats, b_lons):
171 '''
172 Cosine of the angular distance between two points ``a`` and ``b`` on a
173 sphere.
175 This function returns the cosines of the distance
176 angles *delta* between two points ``a`` and ``b`` given as
177 :py:class:`numpy.ndarray`.
178 The coordinates are expected to be given in geographical coordinates
179 and in degrees. For numerical stability a maximum of ``1.0`` is enforced.
181 Please find the details of the implementation in the documentation of
182 the function :py:func:`pyrocko.orthodrome.cosdelta` above.
184 :param a_lats: Latitudes in [deg] point A.
185 :param a_lons: Longitudes in [deg] point A.
186 :param b_lats: Latitudes in [deg] point B.
187 :param b_lons: Longitudes in [deg] point B.
188 :type a_lats: :py:class:`numpy.ndarray`
189 :type a_lons: :py:class:`numpy.ndarray`
190 :type b_lats: :py:class:`numpy.ndarray`
191 :type b_lons: :py:class:`numpy.ndarray`
193 :return: Cosdelta.
194 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
195 '''
196 return num.minimum(
197 1.0,
198 num.sin(a_lats*d2r) * num.sin(b_lats*d2r) +
199 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
200 num.cos(d2r*(b_lons-a_lons)))
203def azimuth(*args):
204 '''
205 Azimuth calculation.
207 This function (find implementation below) returns azimuth ...
208 between points ``a`` and ``b``, coordinates of
209 which are expected to be given in geographical coordinates and in degrees.
211 .. math::
213 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
214 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
215 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
216 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\
218 \\varphi_{\\mathrm{azi},AB} = \\frac{180}{\\pi} \\arctan \\left[
219 \\frac{
220 \\cos( A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
221 \\sin(B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )}
222 {\\sin ( B_{\\mathrm{lat'}} ) - \\sin( A_{\\mathrm{lat'}}
223 cosdelta) } \\right]
225 :param a: Location point A.
226 :type a: :py:class:`pyrocko.orthodrome.Loc`
227 :param b: Location point B.
228 :type b: :py:class:`pyrocko.orthodrome.Loc`
230 :return: Azimuth in degree
231 '''
233 alat, alon, blat, blon = _latlon_pair(args)
235 return r2d*math.atan2(
236 math.cos(alat*d2r) * math.cos(blat*d2r) *
237 math.sin(d2r*(blon-alon)),
238 math.sin(d2r*blat) - math.sin(d2r*alat) * cosdelta(
239 alat, alon, blat, blon))
242def azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta=None):
243 '''
244 Calculation of the azimuth (*track angle*) from a location A towards B.
246 This function returns azimuths (*track angles*) from locations A towards B
247 given in :py:class:`numpy.ndarray`. Coordinates are expected to be given in
248 geographical coordinates and in degrees.
250 Please find the details of the implementation in the documentation of the
251 function :py:func:`pyrocko.orthodrome.azimuth`.
254 :param a_lats: Latitudes in [deg] point A.
255 :param a_lons: Longitudes in [deg] point A.
256 :param b_lats: Latitudes in [deg] point B.
257 :param b_lons: Longitudes in [deg] point B.
258 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
259 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
260 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
261 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
263 :return: Azimuths in [deg].
264 :rtype: :py:class:`numpy.ndarray`, ``(N)``
265 '''
266 if _cosdelta is None:
267 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
269 return r2d*num.arctan2(
270 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
271 num.sin(d2r*(b_lons-a_lons)),
272 num.sin(d2r*b_lats) - num.sin(d2r*a_lats) * _cosdelta)
275def azibazi(*args, **kwargs):
276 '''
277 Azimuth and backazimuth from location A towards B and back.
279 :returns: Azimuth in [deg] from A to B, back azimuth in [deg] from B to A.
280 :rtype: tuple[float, float]
281 '''
283 alat, alon, blat, blon = _latlon_pair(args)
284 if alat == blat and alon == blon:
285 return 0., 180.
287 implementation = kwargs.get('implementation', 'c')
288 assert implementation in ('c', 'python')
289 if implementation == 'c':
290 from pyrocko import orthodrome_ext
291 return orthodrome_ext.azibazi(alat, alon, blat, blon)
293 cd = cosdelta(alat, alon, blat, blon)
294 azi = r2d*math.atan2(
295 math.cos(alat*d2r) * math.cos(blat*d2r) *
296 math.sin(d2r*(blon-alon)),
297 math.sin(d2r*blat) - math.sin(d2r*alat) * cd)
298 bazi = r2d*math.atan2(
299 math.cos(blat*d2r) * math.cos(alat*d2r) *
300 math.sin(d2r*(alon-blon)),
301 math.sin(d2r*alat) - math.sin(d2r*blat) * cd)
303 return azi, bazi
306def azibazi_numpy(a_lats, a_lons, b_lats, b_lons, implementation='c'):
307 '''
308 Azimuth and backazimuth from location A towards B and back.
310 Arguments are given as :py:class:`numpy.ndarray`.
312 :param a_lats: Latitude(s) in [deg] of point A.
313 :type a_lats: :py:class:`numpy.ndarray`
314 :param a_lons: Longitude(s) in [deg] of point A.
315 :type a_lons: :py:class:`numpy.ndarray`
316 :param b_lats: Latitude(s) in [deg] of point B.
317 :type b_lats: :py:class:`numpy.ndarray`
318 :param b_lons: Longitude(s) in [deg] of point B.
319 :type b_lons: :py:class:`numpy.ndarray`
321 :returns: Azimuth(s) in [deg] from A to B,
322 back azimuth(s) in [deg] from B to A.
323 :rtype: :py:class:`numpy.ndarray`, :py:class:`numpy.ndarray`
324 '''
326 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
327 a_lats, a_lons, b_lats, b_lons)
329 assert implementation in ('c', 'python')
330 if implementation == 'c':
331 from pyrocko import orthodrome_ext
332 return orthodrome_ext.azibazi_numpy(a_lats, a_lons, b_lats, b_lons)
334 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
335 azis = azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta)
336 bazis = azimuth_numpy(b_lats, b_lons, a_lats, a_lons, _cosdelta)
338 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
339 ii_eq = num.where(eq)[0]
340 azis[ii_eq] = 0.0
341 bazis[ii_eq] = 180.0
342 return azis, bazis
345def azidist_numpy(*args):
346 '''
347 Calculation of the azimuth (*track angle*) and the distance from locations
348 A towards B on a sphere.
350 The assisting functions used are :py:func:`pyrocko.orthodrome.cosdelta` and
351 :py:func:`pyrocko.orthodrome.azimuth`
353 :param a_lats: Latitudes in [deg] point A.
354 :param a_lons: Longitudes in [deg] point A.
355 :param b_lats: Latitudes in [deg] point B.
356 :param b_lons: Longitudes in [deg] point B.
357 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
358 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
359 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
360 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
362 :return: Azimuths in [deg], distances in [deg].
363 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
364 '''
365 _cosdelta = cosdelta_numpy(*args)
366 _azimuths = azimuth_numpy(_cosdelta=_cosdelta, *args)
367 return _azimuths, r2d*num.arccos(_cosdelta)
370def distance_accurate50m(*args, **kwargs):
371 '''
372 Accurate distance calculation based on a spheroid of rotation.
374 Function returns distance in meter between points A and B, coordinates of
375 which must be given in geographical coordinates and in degrees.
376 The returned distance should be accurate to 50 m using WGS84.
377 Values for the Earth's equator radius and the Earth's oblateness
378 (``f_oblate``) are defined in the pyrocko configuration file
379 :py:class:`pyrocko.config`.
381 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
383 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
384 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
386 .. math::
388 F = \\frac{\\pi}{180}
389 \\frac{(A_{lat} + B_{lat})}{2}, \\quad
390 G = \\frac{\\pi}{180}
391 \\frac{(A_{lat} - B_{lat})}{2}, \\quad
392 l = \\frac{\\pi}{180}
393 \\frac{(A_{lon} - B_{lon})}{2} \\quad
394 \\\\[0.5cm]
395 S = \\sin^2(G) \\cdot \\cos^2(l) +
396 \\cos^2(F) \\cdot \\sin^2(l), \\quad \\quad
397 C = \\cos^2(G) \\cdot \\cos^2(l) +
398 \\sin^2(F) \\cdot \\sin^2(l)
400 .. math::
402 w = \\arctan \\left( \\sqrt{ \\frac{S}{C}} \\right) , \\quad
403 r = \\sqrt{\\frac{S}{C} }
405 The spherical-earth distance D between A and B, can be given with:
407 .. math::
409 D_{sphere} = 2w \\cdot R_{equator}
411 The oblateness of the Earth requires some correction with
412 correction factors h1 and h2:
414 .. math::
416 h_1 = \\frac{3r - 1}{2C}, \\quad
417 h_2 = \\frac{3r +1 }{2S}\\\\[0.5cm]
419 D = D_{\\mathrm{sphere}} \\cdot [ 1 + h_1 \\,f_{\\mathrm{oblate}}
420 \\cdot \\sin^2(F)
421 \\cos^2(G) - h_2\\, f_{\\mathrm{oblate}}
422 \\cdot \\cos^2(F) \\sin^2(G)]
425 :param a: Location point A.
426 :type a: :py:class:`pyrocko.orthodrome.Loc`
427 :param b: Location point B.
428 :type b: :py:class:`pyrocko.orthodrome.Loc`
430 :return: Distance in [m].
431 :rtype: float
432 '''
434 alat, alon, blat, blon = _latlon_pair(args)
436 implementation = kwargs.get('implementation', 'c')
437 assert implementation in ('c', 'python')
438 if implementation == 'c':
439 from pyrocko import orthodrome_ext
440 return orthodrome_ext.distance_accurate50m(alat, alon, blat, blon)
442 f = (alat + blat)*d2r / 2.
443 g = (alat - blat)*d2r / 2.
444 h = (alon - blon)*d2r / 2.
446 s = math.sin(g)**2 * math.cos(h)**2 + math.cos(f)**2 * math.sin(h)**2
447 c = math.cos(g)**2 * math.cos(h)**2 + math.sin(f)**2 * math.sin(h)**2
449 w = math.atan(math.sqrt(s/c))
451 if w == 0.0:
452 return 0.0
454 r = math.sqrt(s*c)/w
455 d = 2.*w*earthradius_equator
456 h1 = (3.*r-1.)/(2.*c)
457 h2 = (3.*r+1.)/(2.*s)
459 return d * (1. +
460 earth_oblateness * h1 * math.sin(f)**2 * math.cos(g)**2 -
461 earth_oblateness * h2 * math.cos(f)**2 * math.sin(g)**2)
464def distance_accurate50m_numpy(
465 a_lats, a_lons, b_lats, b_lons, implementation='c'):
467 '''
468 Accurate distance calculation based on a spheroid of rotation.
470 Function returns distance in meter between points ``a`` and ``b``,
471 coordinates of which must be given in geographical coordinates and in
472 degrees.
473 The returned distance should be accurate to 50 m using WGS84.
474 Values for the Earth's equator radius and the Earth's oblateness
475 (``f_oblate``) are defined in the pyrocko configuration file
476 :py:class:`pyrocko.config`.
478 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
480 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
481 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
483 .. math::
485 F_i = \\frac{\\pi}{180}
486 \\frac{(a_{lat,i} + a_{lat,i})}{2}, \\quad
487 G_i = \\frac{\\pi}{180}
488 \\frac{(a_{lat,i} - b_{lat,i})}{2}, \\quad
489 l_i= \\frac{\\pi}{180}
490 \\frac{(a_{lon,i} - b_{lon,i})}{2} \\\\[0.5cm]
491 S_i = \\sin^2(G_i) \\cdot \\cos^2(l_i) +
492 \\cos^2(F_i) \\cdot \\sin^2(l_i), \\quad \\quad
493 C_i = \\cos^2(G_i) \\cdot \\cos^2(l_i) +
494 \\sin^2(F_i) \\cdot \\sin^2(l_i)
496 .. math::
498 w_i = \\arctan \\left( \\sqrt{\\frac{S_i}{C_i}} \\right), \\quad
499 r_i = \\sqrt{\\frac{S_i}{C_i} }
501 The spherical-earth distance ``D`` between ``a`` and ``b``,
502 can be given with:
504 .. math::
506 D_{\\mathrm{sphere},i} = 2w_i \\cdot R_{\\mathrm{equator}}
508 The oblateness of the Earth requires some correction with
509 correction factors ``h1`` and ``h2``:
511 .. math::
513 h_{1.i} = \\frac{3r - 1}{2C_i}, \\quad
514 h_{2,i} = \\frac{3r +1 }{2S_i}\\\\[0.5cm]
516 D_{AB,i} = D_{\\mathrm{sphere},i} \\cdot [1 + h_{1,i}
517 \\,f_{\\mathrm{oblate}}
518 \\cdot \\sin^2(F_i)
519 \\cos^2(G_i) - h_{2,i}\\, f_{\\mathrm{oblate}}
520 \\cdot \\cos^2(F_i) \\sin^2(G_i)]
522 :param a_lats: Latitudes in [deg] point A.
523 :param a_lons: Longitudes in [deg] point A.
524 :param b_lats: Latitudes in [deg] point B.
525 :param b_lons: Longitudes in [deg] point B.
526 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
527 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
528 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
529 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
531 :return: Distances in [m].
532 :rtype: :py:class:`numpy.ndarray`, ``(N)``
533 '''
535 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
536 a_lats, a_lons, b_lats, b_lons)
538 assert implementation in ('c', 'python')
539 if implementation == 'c':
540 from pyrocko import orthodrome_ext
541 return orthodrome_ext.distance_accurate50m_numpy(
542 a_lats, a_lons, b_lats, b_lons)
544 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
545 ii_neq = num.where(num.logical_not(eq))[0]
547 if num.all(eq):
548 return num.zeros_like(eq, dtype=float)
550 def extr(x):
551 if isinstance(x, num.ndarray) and x.size > 1:
552 return x[ii_neq]
553 else:
554 return x
556 a_lats = extr(a_lats)
557 a_lons = extr(a_lons)
558 b_lats = extr(b_lats)
559 b_lons = extr(b_lons)
561 f = (a_lats + b_lats)*d2r / 2.
562 g = (a_lats - b_lats)*d2r / 2.
563 h = (a_lons - b_lons)*d2r / 2.
565 s = num.sin(g)**2 * num.cos(h)**2 + num.cos(f)**2 * num.sin(h)**2
566 c = num.cos(g)**2 * num.cos(h)**2 + num.sin(f)**2 * num.sin(h)**2
568 w = num.arctan(num.sqrt(s/c))
570 r = num.sqrt(s*c)/w
572 d = 2.*w*earthradius_equator
573 h1 = (3.*r-1.)/(2.*c)
574 h2 = (3.*r+1.)/(2.*s)
576 dists = num.zeros(eq.size, dtype=float)
577 dists[ii_neq] = d * (
578 1. +
579 earth_oblateness * h1 * num.sin(f)**2 * num.cos(g)**2 -
580 earth_oblateness * h2 * num.cos(f)**2 * num.sin(g)**2)
582 return dists
585def ne_to_latlon(lat0, lon0, north_m, east_m):
586 '''
587 Transform local cartesian coordinates to latitude and longitude.
589 From east and north coordinates (``x`` and ``y`` coordinate
590 :py:class:`numpy.ndarray`) relative to a reference differences in
591 longitude and latitude are calculated, which are effectively changes in
592 azimuth and distance, respectively:
594 .. math::
596 \\text{distance change:}\\; \\Delta {\\bf{a}} &= \\sqrt{{\\bf{y}}^2 +
597 {\\bf{x}}^2 }/ \\mathrm{R_E},
599 \\text{azimuth change:}\\; \\Delta \\bf{\\gamma} &= \\arctan( \\bf{x}
600 / \\bf{y}).
602 The projection used preserves the azimuths of the input points.
604 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
605 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
606 :param north_m: Northing distances from origin in [m].
607 :param east_m: Easting distances from origin in [m].
608 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
609 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
610 :type lat0: float
611 :type lon0: float
613 :return: Array with latitudes and longitudes in [deg].
614 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
616 '''
618 a = num.sqrt(north_m**2+east_m**2)/earthradius
619 gamma = num.arctan2(east_m, north_m)
621 return azidist_to_latlon_rad(lat0, lon0, gamma, a)
624def azidist_to_latlon(lat0, lon0, azimuth_deg, distance_deg):
625 '''
626 Absolute latitudes and longitudes are calculated from relative changes.
628 Convenience wrapper to :py:func:`azidist_to_latlon_rad` with azimuth and
629 distance given in degrees.
631 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
632 :type lat0: float
633 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
634 :type lon0: float
635 :param azimuth_deg: Azimuth from origin in [deg].
636 :type azimuth_deg: :py:class:`numpy.ndarray`, ``(N)``
637 :param distance_deg: Distances from origin in [deg].
638 :type distance_deg: :py:class:`numpy.ndarray`, ``(N)``
640 :return: Array with latitudes and longitudes in [deg].
641 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
642 '''
644 return azidist_to_latlon_rad(
645 lat0, lon0, azimuth_deg/180.*num.pi, distance_deg/180.*num.pi)
648def azidist_to_latlon_rad(lat0, lon0, azimuth_rad, distance_rad):
649 '''
650 Absolute latitudes and longitudes are calculated from relative changes.
652 For numerical stability a range between of ``-1.0`` and ``1.0`` is
653 enforced for ``c`` and ``alpha``.
655 .. math::
657 \\Delta {\\bf a}_i \\; \\text{and} \\; \\Delta \\gamma_i \\;
658 \\text{are relative distances and azimuths from lat0 and lon0 for
659 \\textit{i} source points of a finite source.}
661 .. math::
663 \\mathrm{b} &= \\frac{\\pi}{2} -\\frac{\\pi}{180}\\;\\mathrm{lat_0}\\\\
664 {\\bf c}_i &=\\arccos[\\; \\cos(\\Delta {\\bf{a}}_i)
665 \\cos(\\mathrm{b}) + |\\Delta \\gamma_i| \\,
666 \\sin(\\Delta {\\bf a}_i)
667 \\sin(\\mathrm{b})\\; ] \\\\
668 \\mathrm{lat}_i &= \\frac{180}{\\pi}
669 \\left(\\frac{\\pi}{2} - {\\bf c}_i \\right)
671 .. math::
673 \\alpha_i &= \\arcsin \\left[ \\; \\frac{ \\sin(\\Delta {\\bf a}_i )
674 \\sin(|\\Delta \\gamma_i|)}{\\sin({\\bf c}_i)}\\;
675 \\right] \\\\
676 \\alpha_i &= \\begin{cases}
677 \\alpha_i, &\\text{if} \\; \\cos(\\Delta {\\bf a}_i) -
678 \\cos(\\mathrm{b}) \\cos({\\bf{c}}_i) > 0, \\;
679 \\text{else} \\\\
680 \\pi - \\alpha_i, & \\text{if} \\; \\alpha_i > 0,\\;
681 \\text{else}\\\\
682 -\\pi - \\alpha_i, & \\text{if} \\; \\alpha_i < 0.
683 \\end{cases} \\\\
684 \\mathrm{lon}_i &= \\mathrm{lon_0} +
685 \\frac{180}{\\pi} \\,
686 \\frac{\\Delta \\gamma_i }{|\\Delta \\gamma_i|}
687 \\cdot \\alpha_i
688 \\text{, with $\\alpha_i \\in [-\\pi,\\pi]$}
690 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
691 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
692 :param distance_rad: Distances from origin in [rad].
693 :param azimuth_rad: Azimuth from origin in [rad].
694 :type distance_rad: :py:class:`numpy.ndarray`, ``(N)``
695 :type azimuth_rad: :py:class:`numpy.ndarray`, ``(N)``
696 :type lat0: float
697 :type lon0: float
699 :return: Array with latitudes and longitudes in [deg].
700 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
701 '''
703 a = distance_rad
704 gamma = azimuth_rad
706 b = math.pi/2.-lat0*d2r
708 alphasign = 1.
709 alphasign = num.where(gamma < 0, -1., 1.)
710 gamma = num.abs(gamma)
712 c = num.arccos(clip(
713 num.cos(a)*num.cos(b)+num.sin(a)*num.sin(b)*num.cos(gamma), -1., 1.))
715 alpha = num.arcsin(clip(
716 num.sin(a)*num.sin(gamma)/num.sin(c), -1., 1.))
718 alpha = num.where(
719 num.cos(a)-num.cos(b)*num.cos(c) < 0,
720 num.where(alpha > 0, math.pi-alpha, -math.pi-alpha),
721 alpha)
723 lat = r2d * (math.pi/2. - c)
724 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
726 return lat, lon
729def ne_to_latlon_alternative_method(lat0, lon0, north_m, east_m):
730 '''
731 Transform local cartesian coordinates to latitude and longitude.
733 Like :py:func:`pyrocko.orthodrome.ne_to_latlon`,
734 but this method (implementation below), although it should be numerically
735 more stable, suffers problems at points which are *across the pole*
736 as seen from the cartesian origin.
738 .. math::
740 \\text{distance change:}\\; \\Delta {{\\bf a}_i} &=
741 \\sqrt{{\\bf{y}}^2_i + {\\bf{x}}^2_i }/ \\mathrm{R_E},\\\\
742 \\text{azimuth change:}\\; \\Delta {\\bf \\gamma}_i &=
743 \\arctan( {\\bf x}_i {\\bf y}_i). \\\\
744 \\mathrm{b} &=
745 \\frac{\\pi}{2} -\\frac{\\pi}{180} \\;\\mathrm{lat_0}\\\\
747 .. math::
749 {{\\bf z}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i -
750 \\mathrm{b}}{2} \\right)}
751 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
752 {{\\bf n}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i +
753 \\mathrm{b}}{2} \\right)}
754 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
755 {{\\bf z}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i -
756 \\mathrm{b}}{2} \\right)}
757 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
758 {{\\bf n}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i +
759 \\mathrm{b}}{2} \\right)}
760 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
761 {{\\bf t}_1}_i &= \\arctan{\\left( \\frac{{{\\bf z}_1}_i}
762 {{{\\bf n}_1}_i} \\right) }\\\\
763 {{\\bf t}_2}_i &= \\arctan{\\left( \\frac{{{\\bf z}_2}_i}
764 {{{\\bf n}_2}_i} \\right) } \\\\[0.5cm]
765 c &= \\begin{cases}
766 2 \\cdot \\arccos \\left( {{\\bf z}_1}_i / \\sin({{\\bf t}_1}_i)
767 \\right),\\; \\text{if }
768 |\\sin({{\\bf t}_1}_i)| >
769 |\\sin({{\\bf t}_2}_i)|,\\; \\text{else} \\\\
770 2 \\cdot \\arcsin{\\left( {{\\bf z}_2}_i /
771 \\sin({{\\bf t}_2}_i) \\right)}.
772 \\end{cases}\\\\
774 .. math::
776 {\\bf {lat}}_i &= \\frac{180}{ \\pi } \\left( \\frac{\\pi}{2}
777 - {\\bf {c}}_i \\right) \\\\
778 {\\bf {lon}}_i &= {\\bf {lon}}_0 + \\frac{180}{ \\pi }
779 \\frac{\\gamma_i}{|\\gamma_i|},
780 \\text{ with}\\; \\gamma \\in [-\\pi,\\pi]
782 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
783 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
784 :param north_m: Northing distances from origin in [m].
785 :param east_m: Easting distances from origin in [m].
786 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
787 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
788 :type lat0: float
789 :type lon0: float
791 :return: Array with latitudes and longitudes in [deg].
792 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
793 '''
795 b = math.pi/2.-lat0*d2r
796 a = num.sqrt(north_m**2+east_m**2)/earthradius
798 gamma = num.arctan2(east_m, north_m)
799 alphasign = 1.
800 alphasign = num.where(gamma < 0., -1., 1.)
801 gamma = num.abs(gamma)
803 z1 = num.cos((a-b)/2.)*num.cos(gamma/2.)
804 n1 = num.cos((a+b)/2.)*num.sin(gamma/2.)
805 z2 = num.sin((a-b)/2.)*num.cos(gamma/2.)
806 n2 = num.sin((a+b)/2.)*num.sin(gamma/2.)
807 t1 = num.arctan2(z1, n1)
808 t2 = num.arctan2(z2, n2)
810 alpha = t1 + t2
812 sin_t1 = num.sin(t1)
813 sin_t2 = num.sin(t2)
814 c = num.where(
815 num.abs(sin_t1) > num.abs(sin_t2),
816 num.arccos(z1/sin_t1)*2.,
817 num.arcsin(z2/sin_t2)*2.)
819 lat = r2d * (math.pi/2. - c)
820 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
821 return lat, lon
824def latlon_to_ne(*args):
825 '''
826 Relative cartesian coordinates with respect to a reference location.
828 For two locations, a reference location A and another location B, given in
829 geographical coordinates in degrees, the corresponding cartesian
830 coordinates are calculated.
831 Assisting functions are :py:func:`pyrocko.orthodrome.azimuth` and
832 :py:func:`pyrocko.orthodrome.distance_accurate50m`.
834 .. math::
836 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
837 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
838 \\mathrm{)}\\\\[0.3cm]
840 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180}
841 \\varphi_{\\mathrm{azi},AB} )\\\\
842 e &= D_{AB} \\cdot
843 \\sin( \\frac{\\pi }{180} \\varphi_{\\mathrm{azi},AB})
845 :param refloc: Location reference point.
846 :type refloc: :py:class:`pyrocko.orthodrome.Loc`
847 :param loc: Location of interest.
848 :type loc: :py:class:`pyrocko.orthodrome.Loc`
850 :return: Northing and easting from refloc to location in [m].
851 :rtype: tuple[float, float]
853 '''
855 azi = azimuth(*args)
856 dist = distance_accurate50m(*args)
857 n, e = math.cos(azi*d2r)*dist, math.sin(azi*d2r)*dist
858 return n, e
861def latlon_to_ne_numpy(lat0, lon0, lat, lon):
862 '''
863 Relative cartesian coordinates with respect to a reference location.
865 For two locations, a reference location (``lat0``, ``lon0``) and another
866 location B, given in geographical coordinates in degrees,
867 the corresponding cartesian coordinates are calculated.
868 Assisting functions are :py:func:`azimuth`
869 and :py:func:`distance_accurate50m`.
871 :param lat0: Latitude of the reference location in [deg].
872 :param lon0: Longitude of the reference location in [deg].
873 :param lat: Latitude of the absolute location in [deg].
874 :param lon: Longitude of the absolute location in [deg].
876 :return: ``(n, e)``: relative north and east positions in [m].
877 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
879 Implemented formulations:
881 .. math::
883 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
884 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
885 \\mathrm{)}\\\\[0.3cm]
887 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180} \\varphi_{
888 \\mathrm{azi},AB} )\\\\
889 e &= D_{AB} \\cdot \\sin( \\frac{\\pi }{180} \\varphi_{
890 \\mathrm{azi},AB} )
891 '''
893 azi = azimuth_numpy(lat0, lon0, lat, lon)
894 dist = distance_accurate50m_numpy(lat0, lon0, lat, lon)
895 n = num.cos(azi*d2r)*dist
896 e = num.sin(azi*d2r)*dist
897 return n, e
900_wgs84 = None
903def get_wgs84():
904 global _wgs84
905 if _wgs84 is None:
906 from geographiclib.geodesic import Geodesic
907 _wgs84 = Geodesic.WGS84
909 return _wgs84
912def amap(n):
913 def wrap(f):
914 if n == 1:
915 def func(*args):
916 it = num.nditer(args + (None,))
917 for ops in it:
918 ops[-1][...] = f(*ops[:-1])
920 return it.operands[-1]
921 elif n == 2:
922 def func(*args):
923 it = num.nditer(args + (None, None))
924 for ops in it:
925 ops[-2][...], ops[-1][...] = f(*ops[:-2])
927 return it.operands[-2], it.operands[-1]
929 return func
930 return wrap
933@amap(2)
934def ne_to_latlon2(lat0, lon0, north_m, east_m):
935 wgs84 = get_wgs84()
936 az = num.arctan2(east_m, north_m)*r2d
937 dist = num.sqrt(east_m**2 + north_m**2)
938 x = wgs84.Direct(lat0, lon0, az, dist)
939 return x['lat2'], x['lon2']
942@amap(2)
943def latlon_to_ne2(lat0, lon0, lat1, lon1):
944 wgs84 = get_wgs84()
945 x = wgs84.Inverse(lat0, lon0, lat1, lon1)
946 dist = x['s12']
947 az = x['azi1']
948 n = num.cos(az*d2r)*dist
949 e = num.sin(az*d2r)*dist
950 return n, e
953@amap(1)
954def distance_accurate15nm(lat1, lon1, lat2, lon2):
955 wgs84 = get_wgs84()
956 return wgs84.Inverse(lat1, lon1, lat2, lon2)['s12']
959def positive_region(region):
960 '''
961 Normalize parameterization of a rectangular geographical region.
963 :param region: ``(west, east, south, north)`` in [deg].
964 :type region: tuple of float
966 :returns: ``(west, east, south, north)``, where ``west <= east`` and
967 where ``west`` and ``east`` are in the range ``[-180., 180.+360.]``.
968 :rtype: tuple of float
969 '''
970 west, east, south, north = [float(x) for x in region]
972 assert -180. - 360. <= west < 180.
973 assert -180. < east <= 180. + 360.
974 assert -90. <= south < 90.
975 assert -90. < north <= 90.
977 if east < west:
978 east += 360.
980 if west < -180.:
981 west += 360.
982 east += 360.
984 return (west, east, south, north)
987def points_in_region(p, region):
988 '''
989 Check what points are contained in a rectangular geographical region.
991 :param p: ``(lat, lon)`` pairs in [deg].
992 :type p: :py:class:`numpy.ndarray` ``(N, 2)``
993 :param region: ``(west, east, south, north)`` region boundaries in [deg].
994 :type region: tuple of float
996 :returns: Mask, returning ``True`` for each point within the region.
997 :rtype: :py:class:`numpy.ndarray` of bool, shape ``(N)``
998 '''
1000 w, e, s, n = positive_region(region)
1001 return num.logical_and(
1002 num.logical_and(s <= p[:, 0], p[:, 0] <= n),
1003 num.logical_or(
1004 num.logical_and(w <= p[:, 1], p[:, 1] <= e),
1005 num.logical_and(w-360. <= p[:, 1], p[:, 1] <= e-360.)))
1008def point_in_region(p, region):
1009 '''
1010 Check if a point is contained in a rectangular geographical region.
1012 :param p: ``(lat, lon)`` in [deg].
1013 :type p: tuple of float
1014 :param region: ``(west, east, south, north)`` region boundaries in [deg].
1015 :type region: tuple of float
1017 :returns: ``True``, if point is in region, else ``False``.
1018 :rtype: bool
1019 '''
1021 w, e, s, n = positive_region(region)
1022 return num.logical_and(
1023 num.logical_and(s <= p[0], p[0] <= n),
1024 num.logical_or(
1025 num.logical_and(w <= p[1], p[1] <= e),
1026 num.logical_and(w-360. <= p[1], p[1] <= e-360.)))
1029def radius_to_region(lat, lon, radius):
1030 '''
1031 Get a rectangular region which fully contains a given circular region.
1033 :param lat: Latitude of the center point of circular region in [deg].
1034 :type lat: float
1035 :param lon: Longitude of the center point of circular region in [deg].
1036 :type lon: float
1037 :param radius: Radius of circular region in [m].
1038 :type radius: float
1040 :returns: Rectangular region as ``(east, west, south, north)`` in [deg] or
1041 ``None``.
1042 :rtype: tuple[float, float, float, float]
1043 '''
1044 radius_deg = radius * m2d
1045 if radius_deg < 45.:
1046 lat_min = max(-90., lat - radius_deg)
1047 lat_max = min(90., lat + radius_deg)
1048 absmaxlat = max(abs(lat_min), abs(lat_max))
1049 if absmaxlat > 89:
1050 lon_min = -180.
1051 lon_max = 180.
1052 else:
1053 lon_min = max(
1054 -180. - 360.,
1055 lon - radius_deg / math.cos(absmaxlat*d2r))
1056 lon_max = min(
1057 180. + 360.,
1058 lon + radius_deg / math.cos(absmaxlat*d2r))
1060 lon_min, lon_max, lat_min, lat_max = positive_region(
1061 (lon_min, lon_max, lat_min, lat_max))
1063 return lon_min, lon_max, lat_min, lat_max
1065 else:
1066 return None
1069def geographic_midpoint(lats, lons, weights=None):
1070 '''
1071 Calculate geographic midpoints by finding the center of gravity.
1073 This method suffers from instabilities if points are centered around the
1074 poles.
1076 :param lats: Latitudes in [deg].
1077 :param lons: Longitudes in [deg].
1078 :param weights: Weighting factors.
1079 :type lats: :py:class:`numpy.ndarray`, ``(N)``
1080 :type lons: :py:class:`numpy.ndarray`, ``(N)``
1081 :type weights: optional, :py:class:`numpy.ndarray`, ``(N)``
1083 :return: Latitudes and longitudes of the midpoints in [deg].
1084 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
1085 '''
1086 if not weights:
1087 weights = num.ones(len(lats))
1089 total_weigth = num.sum(weights)
1090 weights /= total_weigth
1091 lats = lats * d2r
1092 lons = lons * d2r
1093 x = num.sum(num.cos(lats) * num.cos(lons) * weights)
1094 y = num.sum(num.cos(lats) * num.sin(lons) * weights)
1095 z = num.sum(num.sin(lats) * weights)
1097 lon = num.arctan2(y, x)
1098 hyp = num.sqrt(x**2 + y**2)
1099 lat = num.arctan2(z, hyp)
1101 return lat/d2r, lon/d2r
1104def geographic_midpoint_locations(locations, weights=None):
1105 coords = num.array([loc.effective_latlon
1106 for loc in locations])
1107 return geographic_midpoint(coords[:, 0], coords[:, 1], weights)
1110def geodetic_to_ecef(lat, lon, alt):
1111 '''
1112 Convert geodetic coordinates to Earth-Centered, Earth-Fixed (ECEF)
1113 Cartesian coordinates. [#1]_ [#2]_
1115 :param lat: Geodetic latitude in [deg].
1116 :param lon: Geodetic longitude in [deg].
1117 :param alt: Geodetic altitude (height) in [m] (positive for points outside
1118 the geoid).
1119 :type lat: float
1120 :type lon: float
1121 :type alt: float
1123 :return: ECEF Cartesian coordinates (X, Y, Z) in [m].
1124 :rtype: tuple[float, float, float]
1126 .. [#1] https://en.wikipedia.org/wiki/ECEF
1127 .. [#2] https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1128 #From_geodetic_to_ECEF_coordinates
1129 '''
1131 f = earth_oblateness
1132 a = earthradius_equator
1133 e2 = 2*f - f**2
1135 lat, lon = num.radians(lat), num.radians(lon)
1136 # Normal (plumb line)
1137 N = a / num.sqrt(1.0 - (e2 * num.sin(lat)**2))
1139 X = (N+alt) * num.cos(lat) * num.cos(lon)
1140 Y = (N+alt) * num.cos(lat) * num.sin(lon)
1141 Z = (N*(1.0-e2) + alt) * num.sin(lat)
1143 return (X, Y, Z)
1146def ecef_to_geodetic(X, Y, Z):
1147 '''
1148 Convert Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates to
1149 geodetic coordinates (Ferrari's solution).
1151 :param X, Y, Z: Cartesian coordinates in ECEF system in [m].
1152 :type X, Y, Z: float
1154 :return: Geodetic coordinates (lat, lon, alt). Latitude and longitude are
1155 in [deg] and altitude is in [m]
1156 (positive for points outside the geoid).
1157 :rtype: tuple, float
1159 .. seealso ::
1160 https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1161 #The_application_of_Ferrari.27s_solution
1162 '''
1163 f = earth_oblateness
1164 a = earthradius_equator
1165 b = a * (1. - f)
1166 e2 = 2.*f - f**2
1168 # usefull
1169 a2 = a**2
1170 b2 = b**2
1171 e4 = e2**2
1172 X2 = X**2
1173 Y2 = Y**2
1174 Z2 = Z**2
1176 r = num.sqrt(X2 + Y2)
1177 r2 = r**2
1179 e_prime2 = (a2 - b2)/b2
1180 E2 = a2 - b2
1181 F = 54. * b2 * Z2
1182 G = r2 + (1.-e2)*Z2 - (e2*E2)
1183 C = (e4 * F * r2) / (G**3)
1184 S = cbrt(1. + C + num.sqrt(C**2 + 2.*C))
1185 P = F / (3. * (S + 1./S + 1.)**2 * G**2)
1186 Q = num.sqrt(1. + (2.*e4*P))
1188 dum1 = -(P*e2*r) / (1.+Q)
1189 dum2 = 0.5 * a2 * (1. + 1./Q)
1190 dum3 = (P * (1.-e2) * Z2) / (Q * (1.+Q))
1191 dum4 = 0.5 * P * r2
1192 r0 = dum1 + num.sqrt(dum2 - dum3 - dum4)
1194 U = num.sqrt((r - e2*r0)**2 + Z2)
1195 V = num.sqrt((r - e2*r0)**2 + (1.-e2)*Z2)
1196 Z0 = (b2*Z) / (a*V)
1198 alt = U * (1. - (b2 / (a*V)))
1199 lat = num.arctan((Z + e_prime2 * Z0)/r)
1200 lon = num.arctan2(Y, X)
1202 return (lat*r2d, lon*r2d, alt)
1205class Farside(Exception):
1206 pass
1209def latlon_to_xyz(latlons):
1210 if latlons.ndim == 1:
1211 return latlon_to_xyz(latlons[num.newaxis, :])[0]
1213 points = num.zeros((latlons.shape[0], 3))
1214 lats = latlons[:, 0]
1215 lons = latlons[:, 1]
1216 points[:, 0] = num.cos(lats*d2r) * num.cos(lons*d2r)
1217 points[:, 1] = num.cos(lats*d2r) * num.sin(lons*d2r)
1218 points[:, 2] = num.sin(lats*d2r)
1219 return points
1222def xyz_to_latlon(xyz):
1223 if xyz.ndim == 1:
1224 return xyz_to_latlon(xyz[num.newaxis, :])[0]
1226 latlons = num.zeros((xyz.shape[0], 2))
1227 latlons[:, 0] = num.arctan2(
1228 xyz[:, 2], num.sqrt(xyz[:, 0]**2 + xyz[:, 1]**2)) * r2d
1229 latlons[:, 1] = num.arctan2(
1230 xyz[:, 1], xyz[:, 0]) * r2d
1231 return latlons
1234def rot_to_00(lat, lon):
1235 rot0 = euler_to_matrix(0., -90.*d2r, 0.)
1236 rot1 = euler_to_matrix(-d2r*lat, 0., -d2r*lon)
1237 return num.dot(rot0.T, num.dot(rot1, rot0)).T
1240def distances3d(a, b):
1241 return num.sqrt(num.sum((a-b)**2, axis=a.ndim-1))
1244def circulation(points2):
1245 return num.sum(
1246 (points2[1:, 0] - points2[:-1, 0])
1247 * (points2[1:, 1] + points2[:-1, 1]))
1250def stereographic(points):
1251 dists = distances3d(points[1:, :], points[:-1, :])
1252 if dists.size > 0:
1253 maxdist = num.max(dists)
1254 cutoff = maxdist**2 / 2.
1255 else:
1256 cutoff = 1.0e-5
1258 points = points.copy()
1259 if num.any(points[:, 0] < -1. + cutoff):
1260 raise Farside()
1262 points_out = points[:, 1:].copy()
1263 factor = 1.0 / (1.0 + points[:, 0])
1264 points_out *= factor[:, num.newaxis]
1266 return points_out
1269def stereographic_poly(points):
1270 dists = distances3d(points[1:, :], points[:-1, :])
1271 if dists.size > 0:
1272 maxdist = num.max(dists)
1273 cutoff = maxdist**2 / 2.
1274 else:
1275 cutoff = 1.0e-5
1277 points = points.copy()
1278 if num.any(points[:, 0] < -1. + cutoff):
1279 raise Farside()
1281 points_out = points[:, 1:].copy()
1282 factor = 1.0 / (1.0 + points[:, 0])
1283 points_out *= factor[:, num.newaxis]
1285 if circulation(points_out) >= 0:
1286 raise Farside()
1288 return points_out
1291def gnomonic_x(points, cutoff=0.01):
1292 points_out = points[:, 1:].copy()
1293 if num.any(points[:, 0] < cutoff):
1294 raise Farside()
1296 factor = 1.0 / points[:, 0]
1297 points_out *= factor[:, num.newaxis]
1298 return points_out
1301def cneg(i, x):
1302 if i == 1:
1303 return x
1304 else:
1305 return num.logical_not(x)
1308def contains_points(polygon, points):
1309 '''
1310 Test which points are inside polygon on a sphere.
1312 The inside of the polygon is defined as the area which is to the left hand
1313 side of an observer walking the polygon line, points in order, on the
1314 sphere. Lines between the polygon points are treated as great circle paths.
1315 The polygon may be arbitrarily complex, as long as it does not have any
1316 crossings or thin parts with zero width. The polygon may contain the poles
1317 and is allowed to wrap around the sphere multiple times.
1319 The algorithm works by consecutive cutting of the polygon into (almost)
1320 hemispheres and subsequent Gnomonic projections to perform the
1321 point-in-polygon tests on a 2D plane.
1323 :param polygon: Point coordinates defining the polygon [deg].
1324 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1325 0=lat, 1=lon
1326 :param points: Coordinates of points to test [deg].
1327 :type points: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1328 0=lat, 1=lon
1330 :returns: Boolean mask array.
1331 :rtype: :py:class:`numpy.ndarray` of shape ``(N,)``.
1332 '''
1334 and_ = num.logical_and
1336 points_xyz = latlon_to_xyz(points)
1337 mask_x = 0. <= points_xyz[:, 0]
1338 mask_y = 0. <= points_xyz[:, 1]
1339 mask_z = 0. <= points_xyz[:, 2]
1341 result = num.zeros(points.shape[0], dtype=int)
1343 for ix in [-1, 1]:
1344 for iy in [-1, 1]:
1345 for iz in [-1, 1]:
1346 mask = and_(
1347 and_(cneg(ix, mask_x), cneg(iy, mask_y)),
1348 cneg(iz, mask_z))
1350 center_xyz = num.array([ix, iy, iz], dtype=float)
1352 lat, lon = xyz_to_latlon(center_xyz)
1353 rot = rot_to_00(lat, lon)
1355 points_rot_xyz = num.dot(rot, points_xyz[mask, :].T).T
1356 points_rot_pro = gnomonic_x(points_rot_xyz)
1358 offset = 0.01
1360 poly_xyz = latlon_to_xyz(polygon)
1361 poly_rot_xyz = num.dot(rot, poly_xyz.T).T
1362 poly_rot_xyz[:, 0] -= offset
1363 groups = spoly_cut([poly_rot_xyz], axis=0)
1365 for poly_rot_group_xyz in groups[1]:
1366 poly_rot_group_xyz[:, 0] += offset
1368 poly_rot_group_pro = gnomonic_x(
1369 poly_rot_group_xyz)
1371 if circulation(poly_rot_group_pro) > 0:
1372 result[mask] += path_contains_points(
1373 poly_rot_group_pro, points_rot_pro)
1374 else:
1375 result[mask] -= path_contains_points(
1376 poly_rot_group_pro, points_rot_pro)
1378 return result.astype(bool)
1381def contains_point(polygon, point):
1382 '''
1383 Test if point is inside polygon on a sphere.
1385 Convenience wrapper to :py:func:`contains_points` to test a single point.
1387 :param polygon: Point coordinates defining the polygon [deg].
1388 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1389 0=lat, 1=lon
1390 :param point: Coordinates ``(lat, lon)`` of point to test [deg].
1391 :type point: tuple of float
1393 :returns: ``True``, if point is located within polygon, else ``False``.
1394 :rtype: bool
1395 '''
1397 return bool(
1398 contains_points(polygon, num.asarray(point)[num.newaxis, :])[0])