1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5from __future__ import division, absolute_import 

6 

7import math 

8import numpy as num 

9 

10from .moment_tensor import euler_to_matrix 

11from .config import config 

12from .plot.beachball import spoly_cut 

13 

14from matplotlib.path import Path 

15 

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 

23 

24_testpath = Path([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)], closed=True) 

25 

26raise_if_slow_path_contains_points = False 

27 

28 

29class Slow(Exception): 

30 pass 

31 

32 

33if hasattr(_testpath, 'contains_points') and num.all( 

34 _testpath.contains_points([(0.5, 0.5), (1.5, 0.5)]) == [True, False]): 

35 

36 def path_contains_points(verts, points): 

37 p = Path(verts, closed=True) 

38 return p.contains_points(points).astype(bool) 

39 

40else: 

41 # work around missing contains_points and bug in matplotlib ~ v1.2.0 

42 

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() 

47 

48 p = Path(verts, closed=True) 

49 result = num.zeros(points.shape[0], dtype=bool) 

50 for i in range(result.size): 

51 result[i] = p.contains_point(points[i, :]) 

52 

53 return result 

54 

55 

56try: 

57 cbrt = num.cbrt 

58except AttributeError: 

59 def cbrt(x): 

60 return x**(1./3.) 

61 

62 

63def float_array_broadcast(*args): 

64 return num.broadcast_arrays(*[ 

65 num.asarray(x, dtype=float) for x in args]) 

66 

67 

68class Loc(object): 

69 ''' 

70 Simple location representation. 

71 

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 

78 

79 

80def clip(x, mi, ma): 

81 ''' 

82 Clip values of an array. 

83 

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 

90 

91 :return: Clipped data. 

92 :rtype: :py:class:`numpy.ndarray` 

93 ''' 

94 return num.minimum(num.maximum(mi, x), ma) 

95 

96 

97def wrap(x, mi, ma): 

98 ''' 

99 Wrapping continuous data to fundamental phase values. 

100 

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}}]. 

108 

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 

115 

116 :return: Wrapped data. 

117 :rtype: :py:class:`numpy.ndarray` 

118 ''' 

119 return x - num.floor((x-mi)/(ma-mi)) * (ma-mi) 

120 

121 

122def _latlon_pair(args): 

123 if len(args) == 2: 

124 a, b = args 

125 return a.lat, a.lon, b.lat, b.lon 

126 

127 elif len(args) == 4: 

128 return args 

129 

130 

131def cosdelta(*args): 

132 ''' 

133 Cosine of the angular distance between two points ``a`` and ``b`` on a 

134 sphere. 

135 

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. 

140 

141 .. math:: 

142 

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] 

147 

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'}} ) 

152 

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` 

157 

158 :return: Cosdelta. 

159 :rtype: float 

160 ''' 

161 

162 alat, alon, blat, blon = _latlon_pair(args) 

163 

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))) 

169 

170 

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. 

175 

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. 

181 

182 Please find the details of the implementation in the documentation of 

183 the function :py:func:`pyrocko.orthodrome.cosdelta` above. 

184 

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` 

193 

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))) 

202 

203 

204def azimuth(*args): 

205 ''' 

206 Azimuth calculation. 

207 

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. 

211 

212 .. math:: 

213 

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}\\\\ 

218 

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] 

225 

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` 

230 

231 :return: Azimuth in degree 

232 ''' 

233 

234 alat, alon, blat, blon = _latlon_pair(args) 

235 

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)) 

241 

242 

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. 

246 

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. 

250 

251 Please find the details of the implementation in the documentation of the 

252 function :py:func:`pyrocko.orthodrome.azimuth`. 

253 

254 

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)`` 

263 

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) 

269 

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) 

274 

275 

276def azibazi(*args, **kwargs): 

277 ''' 

278 Azimuth and backazimuth from location A towards B and back. 

279 

280 :returns: Azimuth in [deg] from A to B, back azimuth in [deg] from B to A. 

281 :rtype: tuple[float, float] 

282 ''' 

283 

284 alat, alon, blat, blon = _latlon_pair(args) 

285 if alat == blat and alon == blon: 

286 return 0., 180. 

287 

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) 

293 

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) 

303 

304 return azi, bazi 

305 

306 

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. 

310 

311 Arguments are given as :py:class:`numpy.ndarray`. 

312 

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` 

321 

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 ''' 

326 

327 a_lats, a_lons, b_lats, b_lons = float_array_broadcast( 

328 a_lats, a_lons, b_lats, b_lons) 

329 

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) 

334 

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) 

338 

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 

344 

345 

346def azidist_numpy(*args): 

347 ''' 

348 Calculation of the azimuth (*track angle*) and the distance from locations 

349 A towards B on a sphere. 

350 

351 The assisting functions used are :py:func:`pyrocko.orthodrome.cosdelta` and 

352 :py:func:`pyrocko.orthodrome.azimuth` 

353 

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)`` 

362 

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) 

369 

370 

371def distance_accurate50m(*args, **kwargs): 

372 ''' 

373 Accurate distance calculation based on a spheroid of rotation. 

374 

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`. 

381 

382 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on: 

383 

384 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell, 

385 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1`` 

386 

387 .. math:: 

388 

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) 

400 

401 .. math:: 

402 

403 w = \\arctan \\left( \\sqrt{ \\frac{S}{C}} \\right) , \\quad 

404 r = \\sqrt{\\frac{S}{C} } 

405 

406 The spherical-earth distance D between A and B, can be given with: 

407 

408 .. math:: 

409 

410 D_{sphere} = 2w \\cdot R_{equator} 

411 

412 The oblateness of the Earth requires some correction with 

413 correction factors h1 and h2: 

414 

415 .. math:: 

416 

417 h_1 = \\frac{3r - 1}{2C}, \\quad 

418 h_2 = \\frac{3r +1 }{2S}\\\\[0.5cm] 

419 

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)] 

424 

425 

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` 

430 

431 :return: Distance in [m]. 

432 :rtype: float 

433 ''' 

434 

435 alat, alon, blat, blon = _latlon_pair(args) 

436 

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) 

442 

443 f = (alat + blat)*d2r / 2. 

444 g = (alat - blat)*d2r / 2. 

445 h = (alon - blon)*d2r / 2. 

446 

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 

449 

450 w = math.atan(math.sqrt(s/c)) 

451 

452 if w == 0.0: 

453 return 0.0 

454 

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) 

459 

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) 

463 

464 

465def distance_accurate50m_numpy( 

466 a_lats, a_lons, b_lats, b_lons, implementation='c'): 

467 

468 ''' 

469 Accurate distance calculation based on a spheroid of rotation. 

470 

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`. 

478 

479 From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on: 

480 

481 ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell, 

482 Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1`` 

483 

484 .. math:: 

485 

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) 

496 

497 .. math:: 

498 

499 w_i = \\arctan \\left( \\sqrt{\\frac{S_i}{C_i}} \\right), \\quad 

500 r_i = \\sqrt{\\frac{S_i}{C_i} } 

501 

502 The spherical-earth distance ``D`` between ``a`` and ``b``, 

503 can be given with: 

504 

505 .. math:: 

506 

507 D_{\\mathrm{sphere},i} = 2w_i \\cdot R_{\\mathrm{equator}} 

508 

509 The oblateness of the Earth requires some correction with 

510 correction factors ``h1`` and ``h2``: 

511 

512 .. math:: 

513 

514 h_{1.i} = \\frac{3r - 1}{2C_i}, \\quad 

515 h_{2,i} = \\frac{3r +1 }{2S_i}\\\\[0.5cm] 

516 

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)] 

522 

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)`` 

531 

532 :return: Distances in [m]. 

533 :rtype: :py:class:`numpy.ndarray`, ``(N)`` 

534 ''' 

535 

536 a_lats, a_lons, b_lats, b_lons = float_array_broadcast( 

537 a_lats, a_lons, b_lats, b_lons) 

538 

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) 

544 

545 eq = num.logical_and(a_lats == b_lats, a_lons == b_lons) 

546 ii_neq = num.where(num.logical_not(eq))[0] 

547 

548 if num.all(eq): 

549 return num.zeros_like(eq, dtype=float) 

550 

551 def extr(x): 

552 if isinstance(x, num.ndarray) and x.size > 1: 

553 return x[ii_neq] 

554 else: 

555 return x 

556 

557 a_lats = extr(a_lats) 

558 a_lons = extr(a_lons) 

559 b_lats = extr(b_lats) 

560 b_lons = extr(b_lons) 

561 

562 f = (a_lats + b_lats)*d2r / 2. 

563 g = (a_lats - b_lats)*d2r / 2. 

564 h = (a_lons - b_lons)*d2r / 2. 

565 

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 

568 

569 w = num.arctan(num.sqrt(s/c)) 

570 

571 r = num.sqrt(s*c)/w 

572 

573 d = 2.*w*earthradius_equator 

574 h1 = (3.*r-1.)/(2.*c) 

575 h2 = (3.*r+1.)/(2.*s) 

576 

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) 

582 

583 return dists 

584 

585 

586def ne_to_latlon(lat0, lon0, north_m, east_m): 

587 ''' 

588 Transform local cartesian coordinates to latitude and longitude. 

589 

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: 

594 

595 .. math:: 

596 

597 \\text{distance change:}\\; \\Delta {\\bf{a}} &= \\sqrt{{\\bf{y}}^2 + 

598 {\\bf{x}}^2 }/ \\mathrm{R_E}, 

599 

600 \\text{azimuth change:}\\; \\Delta \\bf{\\gamma} &= \\arctan( \\bf{x} 

601 / \\bf{y}). 

602 

603 The projection used preserves the azimuths of the input points. 

604 

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 

613 

614 :return: Array with latitudes and longitudes in [deg]. 

615 :rtype: :py:class:`numpy.ndarray`, ``(2xN)`` 

616 

617 ''' 

618 

619 a = num.sqrt(north_m**2+east_m**2)/earthradius 

620 gamma = num.arctan2(east_m, north_m) 

621 

622 return azidist_to_latlon_rad(lat0, lon0, gamma, a) 

623 

624 

625def azidist_to_latlon(lat0, lon0, azimuth_deg, distance_deg): 

626 ''' 

627 Absolute latitudes and longitudes are calculated from relative changes. 

628 

629 Convenience wrapper to :py:func:`azidist_to_latlon_rad` with azimuth and 

630 distance given in degrees. 

631 

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)`` 

640 

641 :return: Array with latitudes and longitudes in [deg]. 

642 :rtype: :py:class:`numpy.ndarray`, ``(2xN)`` 

643 ''' 

644 

645 return azidist_to_latlon_rad( 

646 lat0, lon0, azimuth_deg/180.*num.pi, distance_deg/180.*num.pi) 

647 

648 

649def azidist_to_latlon_rad(lat0, lon0, azimuth_rad, distance_rad): 

650 ''' 

651 Absolute latitudes and longitudes are calculated from relative changes. 

652 

653 For numerical stability a range between of ``-1.0`` and ``1.0`` is 

654 enforced for ``c`` and ``alpha``. 

655 

656 .. math:: 

657 

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.} 

661 

662 .. math:: 

663 

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) 

671 

672 .. math:: 

673 

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]$} 

690 

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 

699 

700 :return: Array with latitudes and longitudes in [deg]. 

701 :rtype: :py:class:`numpy.ndarray`, ``(2xN)`` 

702 ''' 

703 

704 a = distance_rad 

705 gamma = azimuth_rad 

706 

707 b = math.pi/2.-lat0*d2r 

708 

709 alphasign = 1. 

710 alphasign = num.where(gamma < 0, -1., 1.) 

711 gamma = num.abs(gamma) 

712 

713 c = num.arccos(clip( 

714 num.cos(a)*num.cos(b)+num.sin(a)*num.sin(b)*num.cos(gamma), -1., 1.)) 

715 

716 alpha = num.arcsin(clip( 

717 num.sin(a)*num.sin(gamma)/num.sin(c), -1., 1.)) 

718 

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) 

723 

724 lat = r2d * (math.pi/2. - c) 

725 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.) 

726 

727 return lat, lon 

728 

729 

730def ne_to_latlon_alternative_method(lat0, lon0, north_m, east_m): 

731 ''' 

732 Transform local cartesian coordinates to latitude and longitude. 

733 

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. 

738 

739 .. math:: 

740 

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}\\\\ 

747 

748 .. math:: 

749 

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}\\\\ 

774 

775 .. math:: 

776 

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] 

782 

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 

791 

792 :return: Array with latitudes and longitudes in [deg]. 

793 :rtype: :py:class:`numpy.ndarray`, ``(2xN)`` 

794 ''' 

795 

796 b = math.pi/2.-lat0*d2r 

797 a = num.sqrt(north_m**2+east_m**2)/earthradius 

798 

799 gamma = num.arctan2(east_m, north_m) 

800 alphasign = 1. 

801 alphasign = num.where(gamma < 0., -1., 1.) 

802 gamma = num.abs(gamma) 

803 

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) 

810 

811 alpha = t1 + t2 

812 

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.) 

819 

820 lat = r2d * (math.pi/2. - c) 

821 lon = wrap(lon0 + r2d*alpha*alphasign, -180., 180.) 

822 return lat, lon 

823 

824 

825def latlon_to_ne(*args): 

826 ''' 

827 Relative cartesian coordinates with respect to a reference location. 

828 

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`. 

834 

835 .. math:: 

836 

837 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad 

838 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B 

839 \\mathrm{)}\\\\[0.3cm] 

840 

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}) 

845 

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` 

850 

851 :return: Northing and easting from refloc to location in [m]. 

852 :rtype: tuple[float, float] 

853 

854 ''' 

855 

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 

860 

861 

862def latlon_to_ne_numpy(lat0, lon0, lat, lon): 

863 ''' 

864 Relative cartesian coordinates with respect to a reference location. 

865 

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`. 

871 

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]. 

876 

877 :return: ``(n, e)``: relative north and east positions in [m]. 

878 :rtype: :py:class:`numpy.ndarray`, ``(2xN)`` 

879 

880 Implemented formulations: 

881 

882 .. math:: 

883 

884 D_{AB} &= \\mathrm{distance\\_accurate50m(}A, B \\mathrm{)}, \\quad 

885 \\varphi_{\\mathrm{azi},AB} = \\mathrm{azimuth(}A,B 

886 \\mathrm{)}\\\\[0.3cm] 

887 

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 ''' 

893 

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 

899 

900 

901_wgs84 = None 

902 

903 

904def get_wgs84(): 

905 global _wgs84 

906 if _wgs84 is None: 

907 from geographiclib.geodesic import Geodesic 

908 _wgs84 = Geodesic.WGS84 

909 

910 return _wgs84 

911 

912 

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]) 

920 

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]) 

927 

928 return it.operands[-2], it.operands[-1] 

929 

930 return func 

931 return wrap 

932 

933 

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'] 

941 

942 

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 

952 

953 

954@amap(1) 

955def distance_accurate15nm(lat1, lon1, lat2, lon2): 

956 wgs84 = get_wgs84() 

957 return wgs84.Inverse(lat1, lon1, lat2, lon2)['s12'] 

958 

959 

960def positive_region(region): 

961 ''' 

962 Normalize parameterization of a rectangular geographical region. 

963 

964 :param region: ``(west, east, south, north)`` in [deg]. 

965 :type region: tuple of float 

966 

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] 

972 

973 assert -180. - 360. <= west < 180. 

974 assert -180. < east <= 180. + 360. 

975 assert -90. <= south < 90. 

976 assert -90. < north <= 90. 

977 

978 if east < west: 

979 east += 360. 

980 

981 if west < -180.: 

982 west += 360. 

983 east += 360. 

984 

985 return (west, east, south, north) 

986 

987 

988def points_in_region(p, region): 

989 ''' 

990 Check what points are contained in a rectangular geographical region. 

991 

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 

996 

997 :returns: Mask, returning ``True`` for each point within the region. 

998 :rtype: :py:class:`numpy.ndarray` of bool, shape ``(N)`` 

999 ''' 

1000 

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.))) 

1007 

1008 

1009def point_in_region(p, region): 

1010 ''' 

1011 Check if a point is contained in a rectangular geographical region. 

1012 

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 

1017 

1018 :returns: ``True``, if point is in region, else ``False``. 

1019 :rtype: bool 

1020 ''' 

1021 

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.))) 

1028 

1029 

1030def radius_to_region(lat, lon, radius): 

1031 ''' 

1032 Get a rectangular region which fully contains a given circular region. 

1033 

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 

1040 

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)) 

1060 

1061 lon_min, lon_max, lat_min, lat_max = positive_region( 

1062 (lon_min, lon_max, lat_min, lat_max)) 

1063 

1064 return lon_min, lon_max, lat_min, lat_max 

1065 

1066 else: 

1067 return None 

1068 

1069 

1070def geographic_midpoint(lats, lons, weights=None): 

1071 ''' 

1072 Calculate geographic midpoints by finding the center of gravity. 

1073 

1074 This method suffers from instabilities if points are centered around the 

1075 poles. 

1076 

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)`` 

1083 

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)) 

1089 

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) 

1097 

1098 lon = num.arctan2(y, x) 

1099 hyp = num.sqrt(x**2 + y**2) 

1100 lat = num.arctan2(z, hyp) 

1101 

1102 return lat/d2r, lon/d2r 

1103 

1104 

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) 

1109 

1110 

1111def geodetic_to_ecef(lat, lon, alt): 

1112 ''' 

1113 Convert geodetic coordinates to Earth-Centered, Earth-Fixed (ECEF) 

1114 Cartesian coordinates. [#1]_ [#2]_ 

1115 

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 

1123 

1124 :return: ECEF Cartesian coordinates (X, Y, Z) in [m]. 

1125 :rtype: tuple[float, float, float] 

1126 

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 ''' 

1131 

1132 f = earth_oblateness 

1133 a = earthradius_equator 

1134 e2 = 2*f - f**2 

1135 

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)) 

1139 

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) 

1143 

1144 return (X, Y, Z) 

1145 

1146 

1147def ecef_to_geodetic(X, Y, Z): 

1148 ''' 

1149 Convert Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates to 

1150 geodetic coordinates (Ferrari's solution). 

1151 

1152 :param X, Y, Z: Cartesian coordinates in ECEF system in [m]. 

1153 :type X, Y, Z: float 

1154 

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 

1159 

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 

1168 

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 

1176 

1177 r = num.sqrt(X2 + Y2) 

1178 r2 = r**2 

1179 

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)) 

1188 

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) 

1194 

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) 

1198 

1199 alt = U * (1. - (b2 / (a*V))) 

1200 lat = num.arctan((Z + e_prime2 * Z0)/r) 

1201 lon = num.arctan2(Y, X) 

1202 

1203 return (lat*r2d, lon*r2d, alt) 

1204 

1205 

1206class Farside(Exception): 

1207 pass 

1208 

1209 

1210def latlon_to_xyz(latlons): 

1211 if latlons.ndim == 1: 

1212 return latlon_to_xyz(latlons[num.newaxis, :])[0] 

1213 

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 

1221 

1222 

1223def xyz_to_latlon(xyz): 

1224 if xyz.ndim == 1: 

1225 return xyz_to_latlon(xyz[num.newaxis, :])[0] 

1226 

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 

1233 

1234 

1235def rot_to_00(lat, lon): 

1236 rot0 = euler_to_matrix(0., -90.*d2r, 0.) 

1237 rot1 = euler_to_matrix(-d2r*lat, 0., -d2r*lon) 

1238 return num.dot(rot0.T, num.dot(rot1, rot0)).T 

1239 

1240 

1241def distances3d(a, b): 

1242 return num.sqrt(num.sum((a-b)**2, axis=a.ndim-1)) 

1243 

1244 

1245def circulation(points2): 

1246 return num.sum( 

1247 (points2[1:, 0] - points2[:-1, 0]) 

1248 * (points2[1:, 1] + points2[:-1, 1])) 

1249 

1250 

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 

1258 

1259 points = points.copy() 

1260 if num.any(points[:, 0] < -1. + cutoff): 

1261 raise Farside() 

1262 

1263 points_out = points[:, 1:].copy() 

1264 factor = 1.0 / (1.0 + points[:, 0]) 

1265 points_out *= factor[:, num.newaxis] 

1266 

1267 return points_out 

1268 

1269 

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 

1277 

1278 points = points.copy() 

1279 if num.any(points[:, 0] < -1. + cutoff): 

1280 raise Farside() 

1281 

1282 points_out = points[:, 1:].copy() 

1283 factor = 1.0 / (1.0 + points[:, 0]) 

1284 points_out *= factor[:, num.newaxis] 

1285 

1286 if circulation(points_out) >= 0: 

1287 raise Farside() 

1288 

1289 return points_out 

1290 

1291 

1292def gnomonic_x(points, cutoff=0.01): 

1293 points_out = points[:, 1:].copy() 

1294 if num.any(points[:, 0] < cutoff): 

1295 raise Farside() 

1296 

1297 factor = 1.0 / points[:, 0] 

1298 points_out *= factor[:, num.newaxis] 

1299 return points_out 

1300 

1301 

1302def cneg(i, x): 

1303 if i == 1: 

1304 return x 

1305 else: 

1306 return num.logical_not(x) 

1307 

1308 

1309def contains_points(polygon, points): 

1310 ''' 

1311 Test which points are inside polygon on a sphere. 

1312 

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. 

1319 

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. 

1323 

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 

1330 

1331 :returns: Boolean mask array. 

1332 :rtype: :py:class:`numpy.ndarray` of shape ``(N,)``. 

1333 ''' 

1334 

1335 and_ = num.logical_and 

1336 

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] 

1341 

1342 result = num.zeros(points.shape[0], dtype=int) 

1343 

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)) 

1350 

1351 center_xyz = num.array([ix, iy, iz], dtype=float) 

1352 

1353 lat, lon = xyz_to_latlon(center_xyz) 

1354 rot = rot_to_00(lat, lon) 

1355 

1356 points_rot_xyz = num.dot(rot, points_xyz[mask, :].T).T 

1357 points_rot_pro = gnomonic_x(points_rot_xyz) 

1358 

1359 offset = 0.01 

1360 

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) 

1365 

1366 for poly_rot_group_xyz in groups[1]: 

1367 poly_rot_group_xyz[:, 0] += offset 

1368 

1369 poly_rot_group_pro = gnomonic_x( 

1370 poly_rot_group_xyz) 

1371 

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) 

1378 

1379 return result.astype(bool) 

1380 

1381 

1382def contains_point(polygon, point): 

1383 ''' 

1384 Test if point is inside polygon on a sphere. 

1385 

1386 Convenience wrapper to :py:func:`contains_points` to test a single point. 

1387 

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 

1393 

1394 :returns: ``True``, if point is located within polygon, else ``False``. 

1395 :rtype: bool 

1396 ''' 

1397 

1398 return bool( 

1399 contains_points(polygon, num.asarray(point)[num.newaxis, :])[0])