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 in [deg].
73 :attrib lon: Longitude in [deg].
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 in [deg] point A.
186 :param a_lons: Longitudes in [deg] point A.
187 :param b_lats: Latitudes in [deg] point B.
188 :param b_lons: Longitudes in [deg] 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 in [deg] point A.
256 :param a_lons: Longitudes in [deg] point A.
257 :param b_lats: Latitudes in [deg] point B.
258 :param b_lons: Longitudes in [deg] 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 [deg].
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 '''
278 Azimuth and backazimuth from location A towards B and back.
280 :returns: Azimuth in [deg] from A to B, back azimuth in [deg] from B to A.
281 :rtype: tuple[float, float]
282 '''
284 alat, alon, blat, blon = _latlon_pair(args)
285 if alat == blat and alon == blon:
286 return 0., 180.
288 implementation = kwargs.get('implementation', 'c')
289 assert implementation in ('c', 'python')
290 if implementation == 'c':
291 from pyrocko import orthodrome_ext
292 return orthodrome_ext.azibazi(alat, alon, blat, blon)
294 cd = cosdelta(alat, alon, blat, blon)
295 azi = r2d*math.atan2(
296 math.cos(alat*d2r) * math.cos(blat*d2r) *
297 math.sin(d2r*(blon-alon)),
298 math.sin(d2r*blat) - math.sin(d2r*alat) * cd)
299 bazi = r2d*math.atan2(
300 math.cos(blat*d2r) * math.cos(alat*d2r) *
301 math.sin(d2r*(alon-blon)),
302 math.sin(d2r*alat) - math.sin(d2r*blat) * cd)
304 return azi, bazi
307def azibazi_numpy(a_lats, a_lons, b_lats, b_lons, implementation='c'):
308 '''
309 Azimuth and backazimuth from location A towards B and back.
311 Arguments are given as :py:class:`numpy.ndarray`.
313 :param a_lats: Latitude(s) in [deg] of point A.
314 :type a_lats: :py:class:`numpy.ndarray`
315 :param a_lons: Longitude(s) in [deg] of point A.
316 :type a_lons: :py:class:`numpy.ndarray`
317 :param b_lats: Latitude(s) in [deg] of point B.
318 :type b_lats: :py:class:`numpy.ndarray`
319 :param b_lons: Longitude(s) in [deg] of point B.
320 :type b_lons: :py:class:`numpy.ndarray`
322 :returns: Azimuth(s) in [deg] from A to B,
323 back azimuth(s) in [deg] from B to A.
324 :rtype: :py:class:`numpy.ndarray`, :py:class:`numpy.ndarray`
325 '''
327 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
328 a_lats, a_lons, b_lats, b_lons)
330 assert implementation in ('c', 'python')
331 if implementation == 'c':
332 from pyrocko import orthodrome_ext
333 return orthodrome_ext.azibazi_numpy(a_lats, a_lons, b_lats, b_lons)
335 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
336 azis = azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta)
337 bazis = azimuth_numpy(b_lats, b_lons, a_lats, a_lons, _cosdelta)
339 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
340 ii_eq = num.where(eq)[0]
341 azis[ii_eq] = 0.0
342 bazis[ii_eq] = 180.0
343 return azis, bazis
346def azidist_numpy(*args):
347 '''
348 Calculation of the azimuth (*track angle*) and the distance from locations
349 A towards B on a sphere.
351 The assisting functions used are :py:func:`pyrocko.orthodrome.cosdelta` and
352 :py:func:`pyrocko.orthodrome.azimuth`
354 :param a_lats: Latitudes in [deg] point A.
355 :param a_lons: Longitudes in [deg] point A.
356 :param b_lats: Latitudes in [deg] point B.
357 :param b_lons: Longitudes in [deg] point B.
358 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
359 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
360 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
361 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
363 :return: Azimuths in [deg], distances in [deg].
364 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
365 '''
366 _cosdelta = cosdelta_numpy(*args)
367 _azimuths = azimuth_numpy(_cosdelta=_cosdelta, *args)
368 return _azimuths, r2d*num.arccos(_cosdelta)
371def distance_accurate50m(*args, **kwargs):
372 '''
373 Accurate distance calculation based on a spheroid of rotation.
375 Function returns distance in meter between points A and B, coordinates of
376 which must be given in geographical coordinates and in degrees.
377 The returned distance should be accurate to 50 m using WGS84.
378 Values for the Earth's equator radius and the Earth's oblateness
379 (``f_oblate``) are defined in the pyrocko configuration file
380 :py:class:`pyrocko.config`.
382 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
384 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
385 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
387 .. math::
389 F = \\frac{\\pi}{180}
390 \\frac{(A_{lat} + B_{lat})}{2}, \\quad
391 G = \\frac{\\pi}{180}
392 \\frac{(A_{lat} - B_{lat})}{2}, \\quad
393 l = \\frac{\\pi}{180}
394 \\frac{(A_{lon} - B_{lon})}{2} \\quad
395 \\\\[0.5cm]
396 S = \\sin^2(G) \\cdot \\cos^2(l) +
397 \\cos^2(F) \\cdot \\sin^2(l), \\quad \\quad
398 C = \\cos^2(G) \\cdot \\cos^2(l) +
399 \\sin^2(F) \\cdot \\sin^2(l)
401 .. math::
403 w = \\arctan \\left( \\sqrt{ \\frac{S}{C}} \\right) , \\quad
404 r = \\sqrt{\\frac{S}{C} }
406 The spherical-earth distance D between A and B, can be given with:
408 .. math::
410 D_{sphere} = 2w \\cdot R_{equator}
412 The oblateness of the Earth requires some correction with
413 correction factors h1 and h2:
415 .. math::
417 h_1 = \\frac{3r - 1}{2C}, \\quad
418 h_2 = \\frac{3r +1 }{2S}\\\\[0.5cm]
420 D = D_{\\mathrm{sphere}} \\cdot [ 1 + h_1 \\,f_{\\mathrm{oblate}}
421 \\cdot \\sin^2(F)
422 \\cos^2(G) - h_2\\, f_{\\mathrm{oblate}}
423 \\cdot \\cos^2(F) \\sin^2(G)]
426 :param a: Location point A.
427 :type a: :py:class:`pyrocko.orthodrome.Loc`
428 :param b: Location point B.
429 :type b: :py:class:`pyrocko.orthodrome.Loc`
431 :return: Distance in [m].
432 :rtype: float
433 '''
435 alat, alon, blat, blon = _latlon_pair(args)
437 implementation = kwargs.get('implementation', 'c')
438 assert implementation in ('c', 'python')
439 if implementation == 'c':
440 from pyrocko import orthodrome_ext
441 return orthodrome_ext.distance_accurate50m(alat, alon, blat, blon)
443 f = (alat + blat)*d2r / 2.
444 g = (alat - blat)*d2r / 2.
445 h = (alon - blon)*d2r / 2.
447 s = math.sin(g)**2 * math.cos(h)**2 + math.cos(f)**2 * math.sin(h)**2
448 c = math.cos(g)**2 * math.cos(h)**2 + math.sin(f)**2 * math.sin(h)**2
450 w = math.atan(math.sqrt(s/c))
452 if w == 0.0:
453 return 0.0
455 r = math.sqrt(s*c)/w
456 d = 2.*w*earthradius_equator
457 h1 = (3.*r-1.)/(2.*c)
458 h2 = (3.*r+1.)/(2.*s)
460 return d * (1. +
461 earth_oblateness * h1 * math.sin(f)**2 * math.cos(g)**2 -
462 earth_oblateness * h2 * math.cos(f)**2 * math.sin(g)**2)
465def distance_accurate50m_numpy(
466 a_lats, a_lons, b_lats, b_lons, implementation='c'):
468 '''
469 Accurate distance calculation based on a spheroid of rotation.
471 Function returns distance in meter between points ``a`` and ``b``,
472 coordinates of which must be given in geographical coordinates and in
473 degrees.
474 The returned distance should be accurate to 50 m using WGS84.
475 Values for the Earth's equator radius and the Earth's oblateness
476 (``f_oblate``) are defined in the pyrocko configuration file
477 :py:class:`pyrocko.config`.
479 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:
481 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
482 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``
484 .. math::
486 F_i = \\frac{\\pi}{180}
487 \\frac{(a_{lat,i} + a_{lat,i})}{2}, \\quad
488 G_i = \\frac{\\pi}{180}
489 \\frac{(a_{lat,i} - b_{lat,i})}{2}, \\quad
490 l_i= \\frac{\\pi}{180}
491 \\frac{(a_{lon,i} - b_{lon,i})}{2} \\\\[0.5cm]
492 S_i = \\sin^2(G_i) \\cdot \\cos^2(l_i) +
493 \\cos^2(F_i) \\cdot \\sin^2(l_i), \\quad \\quad
494 C_i = \\cos^2(G_i) \\cdot \\cos^2(l_i) +
495 \\sin^2(F_i) \\cdot \\sin^2(l_i)
497 .. math::
499 w_i = \\arctan \\left( \\sqrt{\\frac{S_i}{C_i}} \\right), \\quad
500 r_i = \\sqrt{\\frac{S_i}{C_i} }
502 The spherical-earth distance ``D`` between ``a`` and ``b``,
503 can be given with:
505 .. math::
507 D_{\\mathrm{sphere},i} = 2w_i \\cdot R_{\\mathrm{equator}}
509 The oblateness of the Earth requires some correction with
510 correction factors ``h1`` and ``h2``:
512 .. math::
514 h_{1.i} = \\frac{3r - 1}{2C_i}, \\quad
515 h_{2,i} = \\frac{3r +1 }{2S_i}\\\\[0.5cm]
517 D_{AB,i} = D_{\\mathrm{sphere},i} \\cdot [1 + h_{1,i}
518 \\,f_{\\mathrm{oblate}}
519 \\cdot \\sin^2(F_i)
520 \\cos^2(G_i) - h_{2,i}\\, f_{\\mathrm{oblate}}
521 \\cdot \\cos^2(F_i) \\sin^2(G_i)]
523 :param a_lats: Latitudes in [deg] point A.
524 :param a_lons: Longitudes in [deg] point A.
525 :param b_lats: Latitudes in [deg] point B.
526 :param b_lons: Longitudes in [deg] point B.
527 :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
528 :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
529 :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
530 :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
532 :return: Distances in [m].
533 :rtype: :py:class:`numpy.ndarray`, ``(N)``
534 '''
536 a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
537 a_lats, a_lons, b_lats, b_lons)
539 assert implementation in ('c', 'python')
540 if implementation == 'c':
541 from pyrocko import orthodrome_ext
542 return orthodrome_ext.distance_accurate50m_numpy(
543 a_lats, a_lons, b_lats, b_lons)
545 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
546 ii_neq = num.where(num.logical_not(eq))[0]
548 if num.all(eq):
549 return num.zeros_like(eq, dtype=float)
551 def extr(x):
552 if isinstance(x, num.ndarray) and x.size > 1:
553 return x[ii_neq]
554 else:
555 return x
557 a_lats = extr(a_lats)
558 a_lons = extr(a_lons)
559 b_lats = extr(b_lats)
560 b_lons = extr(b_lons)
562 f = (a_lats + b_lats)*d2r / 2.
563 g = (a_lats - b_lats)*d2r / 2.
564 h = (a_lons - b_lons)*d2r / 2.
566 s = num.sin(g)**2 * num.cos(h)**2 + num.cos(f)**2 * num.sin(h)**2
567 c = num.cos(g)**2 * num.cos(h)**2 + num.sin(f)**2 * num.sin(h)**2
569 w = num.arctan(num.sqrt(s/c))
571 r = num.sqrt(s*c)/w
573 d = 2.*w*earthradius_equator
574 h1 = (3.*r-1.)/(2.*c)
575 h2 = (3.*r+1.)/(2.*s)
577 dists = num.zeros(eq.size, dtype=float)
578 dists[ii_neq] = d * (
579 1. +
580 earth_oblateness * h1 * num.sin(f)**2 * num.cos(g)**2 -
581 earth_oblateness * h2 * num.cos(f)**2 * num.sin(g)**2)
583 return dists
586def ne_to_latlon(lat0, lon0, north_m, east_m):
587 '''
588 Transform local cartesian coordinates to latitude and longitude.
590 From east and north coordinates (``x`` and ``y`` coordinate
591 :py:class:`numpy.ndarray`) relative to a reference differences in
592 longitude and latitude are calculated, which are effectively changes in
593 azimuth and distance, respectively:
595 .. math::
597 \\text{distance change:}\\; \\Delta {\\bf{a}} &= \\sqrt{{\\bf{y}}^2 +
598 {\\bf{x}}^2 }/ \\mathrm{R_E},
600 \\text{azimuth change:}\\; \\Delta \\bf{\\gamma} &= \\arctan( \\bf{x}
601 / \\bf{y}).
603 The projection used preserves the azimuths of the input points.
605 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
606 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
607 :param north_m: Northing distances from origin in [m].
608 :param east_m: Easting distances from origin in [m].
609 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
610 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
611 :type lat0: float
612 :type lon0: float
614 :return: Array with latitudes and longitudes in [deg].
615 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
617 '''
619 a = num.sqrt(north_m**2+east_m**2)/earthradius
620 gamma = num.arctan2(east_m, north_m)
622 return azidist_to_latlon_rad(lat0, lon0, gamma, a)
625def azidist_to_latlon(lat0, lon0, azimuth_deg, distance_deg):
626 '''
627 Absolute latitudes and longitudes are calculated from relative changes.
629 Convenience wrapper to :py:func:`azidist_to_latlon_rad` with azimuth and
630 distance given in degrees.
632 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
633 :type lat0: float
634 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
635 :type lon0: float
636 :param azimuth_deg: Azimuth from origin in [deg].
637 :type azimuth_deg: :py:class:`numpy.ndarray`, ``(N)``
638 :param distance_deg: Distances from origin in [deg].
639 :type distance_deg: :py:class:`numpy.ndarray`, ``(N)``
641 :return: Array with latitudes and longitudes in [deg].
642 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
643 '''
645 return azidist_to_latlon_rad(
646 lat0, lon0, azimuth_deg/180.*num.pi, distance_deg/180.*num.pi)
649def azidist_to_latlon_rad(lat0, lon0, azimuth_rad, distance_rad):
650 '''
651 Absolute latitudes and longitudes are calculated from relative changes.
653 For numerical stability a range between of ``-1.0`` and ``1.0`` is
654 enforced for ``c`` and ``alpha``.
656 .. math::
658 \\Delta {\\bf a}_i \\; \\text{and} \\; \\Delta \\gamma_i \\;
659 \\text{are relative distances and azimuths from lat0 and lon0 for
660 \\textit{i} source points of a finite source.}
662 .. math::
664 \\mathrm{b} &= \\frac{\\pi}{2} -\\frac{\\pi}{180}\\;\\mathrm{lat_0}\\\\
665 {\\bf c}_i &=\\arccos[\\; \\cos(\\Delta {\\bf{a}}_i)
666 \\cos(\\mathrm{b}) + |\\Delta \\gamma_i| \\,
667 \\sin(\\Delta {\\bf a}_i)
668 \\sin(\\mathrm{b})\\; ] \\\\
669 \\mathrm{lat}_i &= \\frac{180}{\\pi}
670 \\left(\\frac{\\pi}{2} - {\\bf c}_i \\right)
672 .. math::
674 \\alpha_i &= \\arcsin \\left[ \\; \\frac{ \\sin(\\Delta {\\bf a}_i )
675 \\sin(|\\Delta \\gamma_i|)}{\\sin({\\bf c}_i)}\\;
676 \\right] \\\\
677 \\alpha_i &= \\begin{cases}
678 \\alpha_i, &\\text{if} \\; \\cos(\\Delta {\\bf a}_i) -
679 \\cos(\\mathrm{b}) \\cos({\\bf{c}}_i) > 0, \\;
680 \\text{else} \\\\
681 \\pi - \\alpha_i, & \\text{if} \\; \\alpha_i > 0,\\;
682 \\text{else}\\\\
683 -\\pi - \\alpha_i, & \\text{if} \\; \\alpha_i < 0.
684 \\end{cases} \\\\
685 \\mathrm{lon}_i &= \\mathrm{lon_0} +
686 \\frac{180}{\\pi} \\,
687 \\frac{\\Delta \\gamma_i }{|\\Delta \\gamma_i|}
688 \\cdot \\alpha_i
689 \\text{, with $\\alpha_i \\in [-\\pi,\\pi]$}
691 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
692 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
693 :param distance_rad: Distances from origin in [rad].
694 :param azimuth_rad: Azimuth from origin in [rad].
695 :type distance_rad: :py:class:`numpy.ndarray`, ``(N)``
696 :type azimuth_rad: :py:class:`numpy.ndarray`, ``(N)``
697 :type lat0: float
698 :type lon0: float
700 :return: Array with latitudes and longitudes in [deg].
701 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
702 '''
704 a = distance_rad
705 gamma = azimuth_rad
707 b = math.pi/2.-lat0*d2r
709 alphasign = 1.
710 alphasign = num.where(gamma < 0, -1., 1.)
711 gamma = num.abs(gamma)
713 c = num.arccos(clip(
714 num.cos(a)*num.cos(b)+num.sin(a)*num.sin(b)*num.cos(gamma), -1., 1.))
716 alpha = num.arcsin(clip(
717 num.sin(a)*num.sin(gamma)/num.sin(c), -1., 1.))
719 alpha = num.where(
720 num.cos(a)-num.cos(b)*num.cos(c) < 0,
721 num.where(alpha > 0, math.pi-alpha, -math.pi-alpha),
722 alpha)
724 lat = r2d * (math.pi/2. - c)
725 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
727 return lat, lon
730def ne_to_latlon_alternative_method(lat0, lon0, north_m, east_m):
731 '''
732 Transform local cartesian coordinates to latitude and longitude.
734 Like :py:func:`pyrocko.orthodrome.ne_to_latlon`,
735 but this method (implementation below), although it should be numerically
736 more stable, suffers problems at points which are *across the pole*
737 as seen from the cartesian origin.
739 .. math::
741 \\text{distance change:}\\; \\Delta {{\\bf a}_i} &=
742 \\sqrt{{\\bf{y}}^2_i + {\\bf{x}}^2_i }/ \\mathrm{R_E},\\\\
743 \\text{azimuth change:}\\; \\Delta {\\bf \\gamma}_i &=
744 \\arctan( {\\bf x}_i {\\bf y}_i). \\\\
745 \\mathrm{b} &=
746 \\frac{\\pi}{2} -\\frac{\\pi}{180} \\;\\mathrm{lat_0}\\\\
748 .. math::
750 {{\\bf z}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i -
751 \\mathrm{b}}{2} \\right)}
752 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
753 {{\\bf n}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i +
754 \\mathrm{b}}{2} \\right)}
755 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
756 {{\\bf z}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i -
757 \\mathrm{b}}{2} \\right)}
758 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
759 {{\\bf n}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i +
760 \\mathrm{b}}{2} \\right)}
761 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\
762 {{\\bf t}_1}_i &= \\arctan{\\left( \\frac{{{\\bf z}_1}_i}
763 {{{\\bf n}_1}_i} \\right) }\\\\
764 {{\\bf t}_2}_i &= \\arctan{\\left( \\frac{{{\\bf z}_2}_i}
765 {{{\\bf n}_2}_i} \\right) } \\\\[0.5cm]
766 c &= \\begin{cases}
767 2 \\cdot \\arccos \\left( {{\\bf z}_1}_i / \\sin({{\\bf t}_1}_i)
768 \\right),\\; \\text{if }
769 |\\sin({{\\bf t}_1}_i)| >
770 |\\sin({{\\bf t}_2}_i)|,\\; \\text{else} \\\\
771 2 \\cdot \\arcsin{\\left( {{\\bf z}_2}_i /
772 \\sin({{\\bf t}_2}_i) \\right)}.
773 \\end{cases}\\\\
775 .. math::
777 {\\bf {lat}}_i &= \\frac{180}{ \\pi } \\left( \\frac{\\pi}{2}
778 - {\\bf {c}}_i \\right) \\\\
779 {\\bf {lon}}_i &= {\\bf {lon}}_0 + \\frac{180}{ \\pi }
780 \\frac{\\gamma_i}{|\\gamma_i|},
781 \\text{ with}\\; \\gamma \\in [-\\pi,\\pi]
783 :param lat0: Latitude origin of the cartesian coordinate system in [deg].
784 :param lon0: Longitude origin of the cartesian coordinate system in [deg].
785 :param north_m: Northing distances from origin in [m].
786 :param east_m: Easting distances from origin in [m].
787 :type north_m: :py:class:`numpy.ndarray`, ``(N)``
788 :type east_m: :py:class:`numpy.ndarray`, ``(N)``
789 :type lat0: float
790 :type lon0: float
792 :return: Array with latitudes and longitudes in [deg].
793 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
794 '''
796 b = math.pi/2.-lat0*d2r
797 a = num.sqrt(north_m**2+east_m**2)/earthradius
799 gamma = num.arctan2(east_m, north_m)
800 alphasign = 1.
801 alphasign = num.where(gamma < 0., -1., 1.)
802 gamma = num.abs(gamma)
804 z1 = num.cos((a-b)/2.)*num.cos(gamma/2.)
805 n1 = num.cos((a+b)/2.)*num.sin(gamma/2.)
806 z2 = num.sin((a-b)/2.)*num.cos(gamma/2.)
807 n2 = num.sin((a+b)/2.)*num.sin(gamma/2.)
808 t1 = num.arctan2(z1, n1)
809 t2 = num.arctan2(z2, n2)
811 alpha = t1 + t2
813 sin_t1 = num.sin(t1)
814 sin_t2 = num.sin(t2)
815 c = num.where(
816 num.abs(sin_t1) > num.abs(sin_t2),
817 num.arccos(z1/sin_t1)*2.,
818 num.arcsin(z2/sin_t2)*2.)
820 lat = r2d * (math.pi/2. - c)
821 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.)
822 return lat, lon
825def latlon_to_ne(*args):
826 '''
827 Relative cartesian coordinates with respect to a reference location.
829 For two locations, a reference location A and another location B, given in
830 geographical coordinates in degrees, the corresponding cartesian
831 coordinates are calculated.
832 Assisting functions are :py:func:`pyrocko.orthodrome.azimuth` and
833 :py:func:`pyrocko.orthodrome.distance_accurate50m`.
835 .. math::
837 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
838 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
839 \\mathrm{)}\\\\[0.3cm]
841 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180}
842 \\varphi_{\\mathrm{azi},AB} )\\\\
843 e &= D_{AB} \\cdot
844 \\sin( \\frac{\\pi }{180} \\varphi_{\\mathrm{azi},AB})
846 :param refloc: Location reference point.
847 :type refloc: :py:class:`pyrocko.orthodrome.Loc`
848 :param loc: Location of interest.
849 :type loc: :py:class:`pyrocko.orthodrome.Loc`
851 :return: Northing and easting from refloc to location in [m].
852 :rtype: tuple[float, float]
854 '''
856 azi = azimuth(*args)
857 dist = distance_accurate50m(*args)
858 n, e = math.cos(azi*d2r)*dist, math.sin(azi*d2r)*dist
859 return n, e
862def latlon_to_ne_numpy(lat0, lon0, lat, lon):
863 '''
864 Relative cartesian coordinates with respect to a reference location.
866 For two locations, a reference location (``lat0``, ``lon0``) and another
867 location B, given in geographical coordinates in degrees,
868 the corresponding cartesian coordinates are calculated.
869 Assisting functions are :py:func:`azimuth`
870 and :py:func:`distance_accurate50m`.
872 :param lat0: Latitude of the reference location in [deg].
873 :param lon0: Longitude of the reference location in [deg].
874 :param lat: Latitude of the absolute location in [deg].
875 :param lon: Longitude of the absolute location in [deg].
877 :return: ``(n, e)``: relative north and east positions in [m].
878 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
880 Implemented formulations:
882 .. math::
884 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad
885 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B
886 \\mathrm{)}\\\\[0.3cm]
888 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180} \\varphi_{
889 \\mathrm{azi},AB} )\\\\
890 e &= D_{AB} \\cdot \\sin( \\frac{\\pi }{180} \\varphi_{
891 \\mathrm{azi},AB} )
892 '''
894 azi = azimuth_numpy(lat0, lon0, lat, lon)
895 dist = distance_accurate50m_numpy(lat0, lon0, lat, lon)
896 n = num.cos(azi*d2r)*dist
897 e = num.sin(azi*d2r)*dist
898 return n, e
901_wgs84 = None
904def get_wgs84():
905 global _wgs84
906 if _wgs84 is None:
907 from geographiclib.geodesic import Geodesic
908 _wgs84 = Geodesic.WGS84
910 return _wgs84
913def amap(n):
914 def wrap(f):
915 if n == 1:
916 def func(*args):
917 it = num.nditer(args + (None,))
918 for ops in it:
919 ops[-1][...] = f(*ops[:-1])
921 return it.operands[-1]
922 elif n == 2:
923 def func(*args):
924 it = num.nditer(args + (None, None))
925 for ops in it:
926 ops[-2][...], ops[-1][...] = f(*ops[:-2])
928 return it.operands[-2], it.operands[-1]
930 return func
931 return wrap
934@amap(2)
935def ne_to_latlon2(lat0, lon0, north_m, east_m):
936 wgs84 = get_wgs84()
937 az = num.arctan2(east_m, north_m)*r2d
938 dist = num.sqrt(east_m**2 + north_m**2)
939 x = wgs84.Direct(lat0, lon0, az, dist)
940 return x['lat2'], x['lon2']
943@amap(2)
944def latlon_to_ne2(lat0, lon0, lat1, lon1):
945 wgs84 = get_wgs84()
946 x = wgs84.Inverse(lat0, lon0, lat1, lon1)
947 dist = x['s12']
948 az = x['azi1']
949 n = num.cos(az*d2r)*dist
950 e = num.sin(az*d2r)*dist
951 return n, e
954@amap(1)
955def distance_accurate15nm(lat1, lon1, lat2, lon2):
956 wgs84 = get_wgs84()
957 return wgs84.Inverse(lat1, lon1, lat2, lon2)['s12']
960def positive_region(region):
961 '''
962 Normalize parameterization of a rectangular geographical region.
964 :param region: ``(west, east, south, north)`` in [deg].
965 :type region: tuple of float
967 :returns: ``(west, east, south, north)``, where ``west <= east`` and
968 where ``west`` and ``east`` are in the range ``[-180., 180.+360.]``.
969 :rtype: tuple of float
970 '''
971 west, east, south, north = [float(x) for x in region]
973 assert -180. - 360. <= west < 180.
974 assert -180. < east <= 180. + 360.
975 assert -90. <= south < 90.
976 assert -90. < north <= 90.
978 if east < west:
979 east += 360.
981 if west < -180.:
982 west += 360.
983 east += 360.
985 return (west, east, south, north)
988def points_in_region(p, region):
989 '''
990 Check what points are contained in a rectangular geographical region.
992 :param p: ``(lat, lon)`` pairs in [deg].
993 :type p: :py:class:`numpy.ndarray` ``(N, 2)``
994 :param region: ``(west, east, south, north)`` region boundaries in [deg].
995 :type region: tuple of float
997 :returns: Mask, returning ``True`` for each point within the region.
998 :rtype: :py:class:`numpy.ndarray` of bool, shape ``(N)``
999 '''
1001 w, e, s, n = positive_region(region)
1002 return num.logical_and(
1003 num.logical_and(s <= p[:, 0], p[:, 0] <= n),
1004 num.logical_or(
1005 num.logical_and(w <= p[:, 1], p[:, 1] <= e),
1006 num.logical_and(w-360. <= p[:, 1], p[:, 1] <= e-360.)))
1009def point_in_region(p, region):
1010 '''
1011 Check if a point is contained in a rectangular geographical region.
1013 :param p: ``(lat, lon)`` in [deg].
1014 :type p: tuple of float
1015 :param region: ``(west, east, south, north)`` region boundaries in [deg].
1016 :type region: tuple of float
1018 :returns: ``True``, if point is in region, else ``False``.
1019 :rtype: bool
1020 '''
1022 w, e, s, n = positive_region(region)
1023 return num.logical_and(
1024 num.logical_and(s <= p[0], p[0] <= n),
1025 num.logical_or(
1026 num.logical_and(w <= p[1], p[1] <= e),
1027 num.logical_and(w-360. <= p[1], p[1] <= e-360.)))
1030def radius_to_region(lat, lon, radius):
1031 '''
1032 Get a rectangular region which fully contains a given circular region.
1034 :param lat: Latitude of the center point of circular region in [deg].
1035 :type lat: float
1036 :param lon: Longitude of the center point of circular region in [deg].
1037 :type lon: float
1038 :param radius: Radius of circular region in [m].
1039 :type radius: float
1041 :returns: Rectangular region as ``(east, west, south, north)`` in [deg] or
1042 ``None``.
1043 :rtype: tuple[float, float, float, float]
1044 '''
1045 radius_deg = radius * m2d
1046 if radius_deg < 45.:
1047 lat_min = max(-90., lat - radius_deg)
1048 lat_max = min(90., lat + radius_deg)
1049 absmaxlat = max(abs(lat_min), abs(lat_max))
1050 if absmaxlat > 89:
1051 lon_min = -180.
1052 lon_max = 180.
1053 else:
1054 lon_min = max(
1055 -180. - 360.,
1056 lon - radius_deg / math.cos(absmaxlat*d2r))
1057 lon_max = min(
1058 180. + 360.,
1059 lon + radius_deg / math.cos(absmaxlat*d2r))
1061 lon_min, lon_max, lat_min, lat_max = positive_region(
1062 (lon_min, lon_max, lat_min, lat_max))
1064 return lon_min, lon_max, lat_min, lat_max
1066 else:
1067 return None
1070def geographic_midpoint(lats, lons, weights=None):
1071 '''
1072 Calculate geographic midpoints by finding the center of gravity.
1074 This method suffers from instabilities if points are centered around the
1075 poles.
1077 :param lats: Latitudes in [deg].
1078 :param lons: Longitudes in [deg].
1079 :param weights: Weighting factors.
1080 :type lats: :py:class:`numpy.ndarray`, ``(N)``
1081 :type lons: :py:class:`numpy.ndarray`, ``(N)``
1082 :type weights: optional, :py:class:`numpy.ndarray`, ``(N)``
1084 :return: Latitudes and longitudes of the midpoints in [deg].
1085 :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
1086 '''
1087 if not weights:
1088 weights = num.ones(len(lats))
1090 total_weigth = num.sum(weights)
1091 weights /= total_weigth
1092 lats = lats * d2r
1093 lons = lons * d2r
1094 x = num.sum(num.cos(lats) * num.cos(lons) * weights)
1095 y = num.sum(num.cos(lats) * num.sin(lons) * weights)
1096 z = num.sum(num.sin(lats) * weights)
1098 lon = num.arctan2(y, x)
1099 hyp = num.sqrt(x**2 + y**2)
1100 lat = num.arctan2(z, hyp)
1102 return lat/d2r, lon/d2r
1105def geographic_midpoint_locations(locations, weights=None):
1106 coords = num.array([loc.effective_latlon
1107 for loc in locations])
1108 return geographic_midpoint(coords[:, 0], coords[:, 1], weights)
1111def geodetic_to_ecef(lat, lon, alt):
1112 '''
1113 Convert geodetic coordinates to Earth-Centered, Earth-Fixed (ECEF)
1114 Cartesian coordinates. [#1]_ [#2]_
1116 :param lat: Geodetic latitude in [deg].
1117 :param lon: Geodetic longitude in [deg].
1118 :param alt: Geodetic altitude (height) in [m] (positive for points outside
1119 the geoid).
1120 :type lat: float
1121 :type lon: float
1122 :type alt: float
1124 :return: ECEF Cartesian coordinates (X, Y, Z) in [m].
1125 :rtype: tuple[float, float, float]
1127 .. [#1] https://en.wikipedia.org/wiki/ECEF
1128 .. [#2] https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1129 #From_geodetic_to_ECEF_coordinates
1130 '''
1132 f = earth_oblateness
1133 a = earthradius_equator
1134 e2 = 2*f - f**2
1136 lat, lon = num.radians(lat), num.radians(lon)
1137 # Normal (plumb line)
1138 N = a / num.sqrt(1.0 - (e2 * num.sin(lat)**2))
1140 X = (N+alt) * num.cos(lat) * num.cos(lon)
1141 Y = (N+alt) * num.cos(lat) * num.sin(lon)
1142 Z = (N*(1.0-e2) + alt) * num.sin(lat)
1144 return (X, Y, Z)
1147def ecef_to_geodetic(X, Y, Z):
1148 '''
1149 Convert Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates to
1150 geodetic coordinates (Ferrari's solution).
1152 :param X, Y, Z: Cartesian coordinates in ECEF system in [m].
1153 :type X, Y, Z: float
1155 :return: Geodetic coordinates (lat, lon, alt). Latitude and longitude are
1156 in [deg] and altitude is in [m]
1157 (positive for points outside the geoid).
1158 :rtype: tuple, float
1160 .. seealso ::
1161 https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
1162 #The_application_of_Ferrari.27s_solution
1163 '''
1164 f = earth_oblateness
1165 a = earthradius_equator
1166 b = a * (1. - f)
1167 e2 = 2.*f - f**2
1169 # usefull
1170 a2 = a**2
1171 b2 = b**2
1172 e4 = e2**2
1173 X2 = X**2
1174 Y2 = Y**2
1175 Z2 = Z**2
1177 r = num.sqrt(X2 + Y2)
1178 r2 = r**2
1180 e_prime2 = (a2 - b2)/b2
1181 E2 = a2 - b2
1182 F = 54. * b2 * Z2
1183 G = r2 + (1.-e2)*Z2 - (e2*E2)
1184 C = (e4 * F * r2) / (G**3)
1185 S = cbrt(1. + C + num.sqrt(C**2 + 2.*C))
1186 P = F / (3. * (S + 1./S + 1.)**2 * G**2)
1187 Q = num.sqrt(1. + (2.*e4*P))
1189 dum1 = -(P*e2*r) / (1.+Q)
1190 dum2 = 0.5 * a2 * (1. + 1./Q)
1191 dum3 = (P * (1.-e2) * Z2) / (Q * (1.+Q))
1192 dum4 = 0.5 * P * r2
1193 r0 = dum1 + num.sqrt(dum2 - dum3 - dum4)
1195 U = num.sqrt((r - e2*r0)**2 + Z2)
1196 V = num.sqrt((r - e2*r0)**2 + (1.-e2)*Z2)
1197 Z0 = (b2*Z) / (a*V)
1199 alt = U * (1. - (b2 / (a*V)))
1200 lat = num.arctan((Z + e_prime2 * Z0)/r)
1201 lon = num.arctan2(Y, X)
1203 return (lat*r2d, lon*r2d, alt)
1206class Farside(Exception):
1207 pass
1210def latlon_to_xyz(latlons):
1211 if latlons.ndim == 1:
1212 return latlon_to_xyz(latlons[num.newaxis, :])[0]
1214 points = num.zeros((latlons.shape[0], 3))
1215 lats = latlons[:, 0]
1216 lons = latlons[:, 1]
1217 points[:, 0] = num.cos(lats*d2r) * num.cos(lons*d2r)
1218 points[:, 1] = num.cos(lats*d2r) * num.sin(lons*d2r)
1219 points[:, 2] = num.sin(lats*d2r)
1220 return points
1223def xyz_to_latlon(xyz):
1224 if xyz.ndim == 1:
1225 return xyz_to_latlon(xyz[num.newaxis, :])[0]
1227 latlons = num.zeros((xyz.shape[0], 2))
1228 latlons[:, 0] = num.arctan2(
1229 xyz[:, 2], num.sqrt(xyz[:, 0]**2 + xyz[:, 1]**2)) * r2d
1230 latlons[:, 1] = num.arctan2(
1231 xyz[:, 1], xyz[:, 0]) * r2d
1232 return latlons
1235def rot_to_00(lat, lon):
1236 rot0 = euler_to_matrix(0., -90.*d2r, 0.).A
1237 rot1 = euler_to_matrix(-d2r*lat, 0., -d2r*lon).A
1238 return num.dot(rot0.T, num.dot(rot1, rot0)).T
1241def distances3d(a, b):
1242 return num.sqrt(num.sum((a-b)**2, axis=a.ndim-1))
1245def circulation(points2):
1246 return num.sum(
1247 (points2[1:, 0] - points2[:-1, 0])
1248 * (points2[1:, 1] + points2[:-1, 1]))
1251def stereographic(points):
1252 dists = distances3d(points[1:, :], points[:-1, :])
1253 if dists.size > 0:
1254 maxdist = num.max(dists)
1255 cutoff = maxdist**2 / 2.
1256 else:
1257 cutoff = 1.0e-5
1259 points = points.copy()
1260 if num.any(points[:, 0] < -1. + cutoff):
1261 raise Farside()
1263 points_out = points[:, 1:].copy()
1264 factor = 1.0 / (1.0 + points[:, 0])
1265 points_out *= factor[:, num.newaxis]
1267 return points_out
1270def stereographic_poly(points):
1271 dists = distances3d(points[1:, :], points[:-1, :])
1272 if dists.size > 0:
1273 maxdist = num.max(dists)
1274 cutoff = maxdist**2 / 2.
1275 else:
1276 cutoff = 1.0e-5
1278 points = points.copy()
1279 if num.any(points[:, 0] < -1. + cutoff):
1280 raise Farside()
1282 points_out = points[:, 1:].copy()
1283 factor = 1.0 / (1.0 + points[:, 0])
1284 points_out *= factor[:, num.newaxis]
1286 if circulation(points_out) >= 0:
1287 raise Farside()
1289 return points_out
1292def gnomonic_x(points, cutoff=0.01):
1293 points_out = points[:, 1:].copy()
1294 if num.any(points[:, 0] < cutoff):
1295 raise Farside()
1297 factor = 1.0 / points[:, 0]
1298 points_out *= factor[:, num.newaxis]
1299 return points_out
1302def cneg(i, x):
1303 if i == 1:
1304 return x
1305 else:
1306 return num.logical_not(x)
1309def contains_points(polygon, points):
1310 '''
1311 Test which points are inside polygon on a sphere.
1313 The inside of the polygon is defined as the area which is to the left hand
1314 side of an observer walking the polygon line, points in order, on the
1315 sphere. Lines between the polygon points are treated as great circle paths.
1316 The polygon may be arbitrarily complex, as long as it does not have any
1317 crossings or thin parts with zero width. The polygon may contain the poles
1318 and is allowed to wrap around the sphere multiple times.
1320 The algorithm works by consecutive cutting of the polygon into (almost)
1321 hemispheres and subsequent Gnomonic projections to perform the
1322 point-in-polygon tests on a 2D plane.
1324 :param polygon: Point coordinates defining the polygon [deg].
1325 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1326 0=lat, 1=lon
1327 :param points: Coordinates of points to test [deg].
1328 :type points: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1329 0=lat, 1=lon
1331 :returns: Boolean mask array.
1332 :rtype: :py:class:`numpy.ndarray` of shape ``(N,)``.
1333 '''
1335 and_ = num.logical_and
1337 points_xyz = latlon_to_xyz(points)
1338 mask_x = 0. <= points_xyz[:, 0]
1339 mask_y = 0. <= points_xyz[:, 1]
1340 mask_z = 0. <= points_xyz[:, 2]
1342 result = num.zeros(points.shape[0], dtype=int)
1344 for ix in [-1, 1]:
1345 for iy in [-1, 1]:
1346 for iz in [-1, 1]:
1347 mask = and_(
1348 and_(cneg(ix, mask_x), cneg(iy, mask_y)),
1349 cneg(iz, mask_z))
1351 center_xyz = num.array([ix, iy, iz], dtype=float)
1353 lat, lon = xyz_to_latlon(center_xyz)
1354 rot = rot_to_00(lat, lon)
1356 points_rot_xyz = num.dot(rot, points_xyz[mask, :].T).T
1357 points_rot_pro = gnomonic_x(points_rot_xyz)
1359 offset = 0.01
1361 poly_xyz = latlon_to_xyz(polygon)
1362 poly_rot_xyz = num.dot(rot, poly_xyz.T).T
1363 poly_rot_xyz[:, 0] -= offset
1364 groups = spoly_cut([poly_rot_xyz], axis=0)
1366 for poly_rot_group_xyz in groups[1]:
1367 poly_rot_group_xyz[:, 0] += offset
1369 poly_rot_group_pro = gnomonic_x(
1370 poly_rot_group_xyz)
1372 if circulation(poly_rot_group_pro) > 0:
1373 result[mask] += path_contains_points(
1374 poly_rot_group_pro, points_rot_pro)
1375 else:
1376 result[mask] -= path_contains_points(
1377 poly_rot_group_pro, points_rot_pro)
1379 return result.astype(num.bool)
1382def contains_point(polygon, point):
1383 '''
1384 Test if point is inside polygon on a sphere.
1386 Convenience wrapper to :py:func:`contains_points` to test a single point.
1388 :param polygon: Point coordinates defining the polygon [deg].
1389 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index
1390 0=lat, 1=lon
1391 :param point: Coordinates ``(lat, lon)`` of point to test [deg].
1392 :type point: tuple of float
1394 :returns: ``True``, if point is located within polygon, else ``False``.
1395 :rtype: bool
1396 '''
1398 return bool(
1399 contains_points(polygon, num.asarray(point)[num.newaxis, :])[0])