Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/orthodrome.py: 76%

418 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-23 12:35 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Some basic geodetic functions. 

8''' 

9 

10import math 

11import numpy as num 

12 

13from .moment_tensor import euler_to_matrix 

14from .config import config 

15from .plot.beachball import spoly_cut 

16 

17from matplotlib.path import Path 

18 

19d2r = math.pi/180. 

20r2d = 1./d2r 

21earth_oblateness = 1./298.257223563 

22earthradius_equator = 6378.14 * 1000. 

23earthradius = config().earthradius 

24d2m = earthradius_equator*math.pi/180. 

25m2d = 1./d2m 

26 

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

28 

29raise_if_slow_path_contains_points = False 

30 

31 

32class Slow(Exception): 

33 pass 

34 

35 

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

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

38 

39 def path_contains_points(verts, points): 

40 p = Path(verts, closed=True) 

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

42 

43else: 

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

45 

46 def path_contains_points(verts, points): 

47 if raise_if_slow_path_contains_points: 

48 # used by unit test to skip slow gshhg_example.py 

49 raise Slow() 

50 

51 p = Path(verts, closed=True) 

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

53 for i in range(result.size): 

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

55 

56 return result 

57 

58 

59try: 

60 cbrt = num.cbrt 

61except AttributeError: 

62 def cbrt(x): 

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

64 

65 

66def float_array_broadcast(*args): 

67 return num.broadcast_arrays(*[ 

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

69 

70 

71class Loc(object): 

72 ''' 

73 Simple location representation. 

74 

75 :attrib lat: Latitude in [deg]. 

76 :attrib lon: Longitude in [deg]. 

77 ''' 

78 def __init__(self, lat, lon): 

79 self.lat = lat 

80 self.lon = lon 

81 

82 

83def clip(x, mi, ma): 

84 ''' 

85 Clip values of an array. 

86 

87 :param x: Continunous data to be clipped. 

88 :param mi: Clip minimum. 

89 :param ma: Clip maximum. 

90 :type x: :py:class:`numpy.ndarray` 

91 :type mi: float 

92 :type ma: float 

93 

94 :return: Clipped data. 

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

96 ''' 

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

98 

99 

100def wrap(x, mi, ma): 

101 ''' 

102 Wrapping continuous data to fundamental phase values. 

103 

104 .. math:: 

105 x_{\\mathrm{wrapped}} = x_{\\mathrm{cont},i} - 

106 \\frac{ x_{\\mathrm{cont},i} - r_{\\mathrm{min}} } 

107 { r_{\\mathrm{max}} - r_{\\mathrm{min}}} 

108 \\cdot ( r_{\\mathrm{max}} - r_{\\mathrm{min}}),\\quad 

109 x_{\\mathrm{wrapped}}\\; \\in 

110 \\;[ r_{\\mathrm{min}},\\, r_{\\mathrm{max}}]. 

111 

112 :param x: Continunous data to be wrapped. 

113 :param mi: Minimum value of wrapped data. 

114 :param ma: Maximum value of wrapped data. 

115 :type x: :py:class:`numpy.ndarray` 

116 :type mi: float 

117 :type ma: float 

118 

119 :return: Wrapped data. 

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

121 ''' 

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

123 

124 

125def _latlon_pair(args): 

126 if len(args) == 2: 

127 a, b = args 

128 return a.lat, a.lon, b.lat, b.lon 

129 

130 elif len(args) == 4: 

131 return args 

132 

133 

134def cosdelta(*args): 

135 ''' 

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

137 sphere. 

138 

139 This function (find implementation below) returns the cosine of the 

140 distance angle 'delta' between two points ``a`` and ``b``, coordinates of 

141 which are expected to be given in geographical coordinates and in degrees. 

142 For numerical stability a maximum of 1.0 is enforced. 

143 

144 .. math:: 

145 

146 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad 

147 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad 

148 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad 

149 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\[0.5cm] 

150 

151 \\cos(\\Delta) = \\min( 1.0, \\quad \\sin( A_{\\mathrm{lat'}}) 

152 \\sin( B_{\\mathrm{lat'}} ) + 

153 \\cos(A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} ) 

154 \\cos( B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} ) 

155 

156 :param a: Location point A. 

157 :type a: :py:class:`pyrocko.orthodrome.Loc` 

158 :param b: Location point B. 

159 :type b: :py:class:`pyrocko.orthodrome.Loc` 

160 

161 :return: Cosdelta. 

162 :rtype: float 

163 ''' 

164 

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

166 

167 return min( 

168 1.0, 

169 math.sin(alat*d2r) * math.sin(blat*d2r) + 

170 math.cos(alat*d2r) * math.cos(blat*d2r) * 

171 math.cos(d2r*(blon-alon))) 

172 

173 

174def cosdelta_numpy(a_lats, a_lons, b_lats, b_lons): 

175 ''' 

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

177 sphere. 

178 

179 This function returns the cosines of the distance 

180 angles *delta* between two points ``a`` and ``b`` given as 

181 :py:class:`numpy.ndarray`. 

182 The coordinates are expected to be given in geographical coordinates 

183 and in degrees. For numerical stability a maximum of ``1.0`` is enforced. 

184 

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

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

187 

188 :param a_lats: Latitudes in [deg] point A. 

189 :param a_lons: Longitudes in [deg] point A. 

190 :param b_lats: Latitudes in [deg] point B. 

191 :param b_lons: Longitudes in [deg] point B. 

192 :type a_lats: :py:class:`numpy.ndarray` 

193 :type a_lons: :py:class:`numpy.ndarray` 

194 :type b_lats: :py:class:`numpy.ndarray` 

195 :type b_lons: :py:class:`numpy.ndarray` 

196 

197 :return: Cosdelta. 

198 :type b_lons: :py:class:`numpy.ndarray`, ``(N)`` 

199 ''' 

200 return num.minimum( 

201 1.0, 

202 num.sin(a_lats*d2r) * num.sin(b_lats*d2r) + 

203 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) * 

204 num.cos(d2r*(b_lons-a_lons))) 

205 

206 

207def azimuth(*args): 

208 ''' 

209 Azimuth calculation. 

210 

211 This function (find implementation below) returns azimuth ... 

212 between points ``a`` and ``b``, coordinates of 

213 which are expected to be given in geographical coordinates and in degrees. 

214 

215 .. math:: 

216 

217 A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad 

218 A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad 

219 B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad 

220 B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\ 

221 

222 \\varphi_{\\mathrm{azi},AB} = \\frac{180}{\\pi} \\arctan \\left[ 

223 \\frac{ 

224 \\cos( A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} ) 

225 \\sin(B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )} 

226 {\\sin ( B_{\\mathrm{lat'}} ) - \\sin( A_{\\mathrm{lat'}} 

227 cosdelta) } \\right] 

228 

229 :param a: Location point A. 

230 :type a: :py:class:`pyrocko.orthodrome.Loc` 

231 :param b: Location point B. 

232 :type b: :py:class:`pyrocko.orthodrome.Loc` 

233 

234 :return: Azimuth in degree 

235 ''' 

236 

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

238 

239 return r2d*math.atan2( 

240 math.cos(alat*d2r) * math.cos(blat*d2r) * 

241 math.sin(d2r*(blon-alon)), 

242 math.sin(d2r*blat) - math.sin(d2r*alat) * cosdelta( 

243 alat, alon, blat, blon)) 

244 

245 

246def azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta=None): 

247 ''' 

248 Calculation of the azimuth (*track angle*) from a location A towards B. 

249 

250 This function returns azimuths (*track angles*) from locations A towards B 

251 given in :py:class:`numpy.ndarray`. Coordinates are expected to be given in 

252 geographical coordinates and in degrees. 

253 

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

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

256 

257 

258 :param a_lats: Latitudes in [deg] point A. 

259 :param a_lons: Longitudes in [deg] point A. 

260 :param b_lats: Latitudes in [deg] point B. 

261 :param b_lons: Longitudes in [deg] point B. 

262 :type a_lats: :py:class:`numpy.ndarray`, ``(N)`` 

263 :type a_lons: :py:class:`numpy.ndarray`, ``(N)`` 

264 :type b_lats: :py:class:`numpy.ndarray`, ``(N)`` 

265 :type b_lons: :py:class:`numpy.ndarray`, ``(N)`` 

266 

267 :return: Azimuths in [deg]. 

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

269 ''' 

270 if _cosdelta is None: 

271 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons) 

272 

273 return r2d*num.arctan2( 

274 num.cos(a_lats*d2r) * num.cos(b_lats*d2r) * 

275 num.sin(d2r*(b_lons-a_lons)), 

276 num.sin(d2r*b_lats) - num.sin(d2r*a_lats) * _cosdelta) 

277 

278 

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

280 ''' 

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

282 

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

284 :rtype: tuple[float, float] 

285 ''' 

286 

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

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

289 return 0., 180. 

290 

291 implementation = kwargs.get('implementation', 'c') 

292 assert implementation in ('c', 'python') 

293 if implementation == 'c': 

294 from pyrocko import orthodrome_ext 

295 return orthodrome_ext.azibazi(alat, alon, blat, blon) 

296 

297 cd = cosdelta(alat, alon, blat, blon) 

298 azi = r2d*math.atan2( 

299 math.cos(alat*d2r) * math.cos(blat*d2r) * 

300 math.sin(d2r*(blon-alon)), 

301 math.sin(d2r*blat) - math.sin(d2r*alat) * cd) 

302 bazi = r2d*math.atan2( 

303 math.cos(blat*d2r) * math.cos(alat*d2r) * 

304 math.sin(d2r*(alon-blon)), 

305 math.sin(d2r*alat) - math.sin(d2r*blat) * cd) 

306 

307 return azi, bazi 

308 

309 

310def azibazi_numpy(a_lats, a_lons, b_lats, b_lons, implementation='c'): 

311 ''' 

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

313 

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

315 

316 :param a_lats: Latitude(s) in [deg] of point A. 

317 :type a_lats: :py:class:`numpy.ndarray` 

318 :param a_lons: Longitude(s) in [deg] of point A. 

319 :type a_lons: :py:class:`numpy.ndarray` 

320 :param b_lats: Latitude(s) in [deg] of point B. 

321 :type b_lats: :py:class:`numpy.ndarray` 

322 :param b_lons: Longitude(s) in [deg] of point B. 

323 :type b_lons: :py:class:`numpy.ndarray` 

324 

325 :returns: Azimuth(s) in [deg] from A to B, 

326 back azimuth(s) in [deg] from B to A. 

327 :rtype: :py:class:`numpy.ndarray`, :py:class:`numpy.ndarray` 

328 ''' 

329 

330 a_lats, a_lons, b_lats, b_lons = float_array_broadcast( 

331 a_lats, a_lons, b_lats, b_lons) 

332 

333 assert implementation in ('c', 'python') 

334 if implementation == 'c': 

335 from pyrocko import orthodrome_ext 

336 return orthodrome_ext.azibazi_numpy(a_lats, a_lons, b_lats, b_lons) 

337 

338 _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons) 

339 azis = azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta) 

340 bazis = azimuth_numpy(b_lats, b_lons, a_lats, a_lons, _cosdelta) 

341 

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

343 ii_eq = num.where(eq)[0] 

344 azis[ii_eq] = 0.0 

345 bazis[ii_eq] = 180.0 

346 return azis, bazis 

347 

348 

349def azidist_numpy(*args): 

350 ''' 

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

352 A towards B on a sphere. 

353 

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

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

356 

357 :param a_lats: Latitudes in [deg] point A. 

358 :param a_lons: Longitudes in [deg] point A. 

359 :param b_lats: Latitudes in [deg] point B. 

360 :param b_lons: Longitudes in [deg] point B. 

361 :type a_lats: :py:class:`numpy.ndarray`, ``(N)`` 

362 :type a_lons: :py:class:`numpy.ndarray`, ``(N)`` 

363 :type b_lats: :py:class:`numpy.ndarray`, ``(N)`` 

364 :type b_lons: :py:class:`numpy.ndarray`, ``(N)`` 

365 

366 :return: Azimuths in [deg], distances in [deg]. 

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

368 ''' 

369 _cosdelta = cosdelta_numpy(*args) 

370 _azimuths = azimuth_numpy(_cosdelta=_cosdelta, *args) 

371 return _azimuths, r2d*num.arccos(_cosdelta) 

372 

373 

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

375 ''' 

376 Accurate distance calculation based on a spheroid of rotation. 

377 

378 Function returns distance in meter between points A and B, coordinates of 

379 which must be given in geographical coordinates and in degrees. 

380 The returned distance should be accurate to 50 m using WGS84. 

381 Values for the Earth's equator radius and the Earth's oblateness 

382 (``f_oblate``) are defined in the pyrocko configuration file 

383 :py:class:`pyrocko.config`. 

384 

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

386 

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

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

389 

390 .. math:: 

391 

392 F = \\frac{\\pi}{180} 

393 \\frac{(A_{lat} + B_{lat})}{2}, \\quad 

394 G = \\frac{\\pi}{180} 

395 \\frac{(A_{lat} - B_{lat})}{2}, \\quad 

396 l = \\frac{\\pi}{180} 

397 \\frac{(A_{lon} - B_{lon})}{2} \\quad 

398 \\\\[0.5cm] 

399 S = \\sin^2(G) \\cdot \\cos^2(l) + 

400 \\cos^2(F) \\cdot \\sin^2(l), \\quad \\quad 

401 C = \\cos^2(G) \\cdot \\cos^2(l) + 

402 \\sin^2(F) \\cdot \\sin^2(l) 

403 

404 .. math:: 

405 

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

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

408 

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

410 

411 .. math:: 

412 

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

414 

415 The oblateness of the Earth requires some correction with 

416 correction factors h1 and h2: 

417 

418 .. math:: 

419 

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

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

422 

423 D = D_{\\mathrm{sphere}} \\cdot [ 1 + h_1 \\,f_{\\mathrm{oblate}} 

424 \\cdot \\sin^2(F) 

425 \\cos^2(G) - h_2\\, f_{\\mathrm{oblate}} 

426 \\cdot \\cos^2(F) \\sin^2(G)] 

427 

428 

429 :param a: Location point A. 

430 :type a: :py:class:`pyrocko.orthodrome.Loc` 

431 :param b: Location point B. 

432 :type b: :py:class:`pyrocko.orthodrome.Loc` 

433 

434 :return: Distance in [m]. 

435 :rtype: float 

436 ''' 

437 

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

439 

440 implementation = kwargs.get('implementation', 'c') 

441 assert implementation in ('c', 'python') 

442 if implementation == 'c': 

443 from pyrocko import orthodrome_ext 

444 return orthodrome_ext.distance_accurate50m(alat, alon, blat, blon) 

445 

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

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

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

449 

450 s = math.sin(g)**2 * math.cos(h)**2 + math.cos(f)**2 * math.sin(h)**2 

451 c = math.cos(g)**2 * math.cos(h)**2 + math.sin(f)**2 * math.sin(h)**2 

452 

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

454 

455 if w == 0.0: 

456 return 0.0 

457 

458 r = math.sqrt(s*c)/w 

459 d = 2.*w*earthradius_equator 

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

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

462 

463 return d * (1. + 

464 earth_oblateness * h1 * math.sin(f)**2 * math.cos(g)**2 - 

465 earth_oblateness * h2 * math.cos(f)**2 * math.sin(g)**2) 

466 

467 

468def distance_accurate50m_numpy( 

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

470 

471 ''' 

472 Accurate distance calculation based on a spheroid of rotation. 

473 

474 Function returns distance in meter between points ``a`` and ``b``, 

475 coordinates of which must be given in geographical coordinates and in 

476 degrees. 

477 The returned distance should be accurate to 50 m using WGS84. 

478 Values for the Earth's equator radius and the Earth's oblateness 

479 (``f_oblate``) are defined in the pyrocko configuration file 

480 :py:class:`pyrocko.config`. 

481 

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

483 

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

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

486 

487 .. math:: 

488 

489 F_i = \\frac{\\pi}{180} 

490 \\frac{(a_{lat,i} + a_{lat,i})}{2}, \\quad 

491 G_i = \\frac{\\pi}{180} 

492 \\frac{(a_{lat,i} - b_{lat,i})}{2}, \\quad 

493 l_i= \\frac{\\pi}{180} 

494 \\frac{(a_{lon,i} - b_{lon,i})}{2} \\\\[0.5cm] 

495 S_i = \\sin^2(G_i) \\cdot \\cos^2(l_i) + 

496 \\cos^2(F_i) \\cdot \\sin^2(l_i), \\quad \\quad 

497 C_i = \\cos^2(G_i) \\cdot \\cos^2(l_i) + 

498 \\sin^2(F_i) \\cdot \\sin^2(l_i) 

499 

500 .. math:: 

501 

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

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

504 

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

506 can be given with: 

507 

508 .. math:: 

509 

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

511 

512 The oblateness of the Earth requires some correction with 

513 correction factors ``h1`` and ``h2``: 

514 

515 .. math:: 

516 

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

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

519 

520 D_{AB,i} = D_{\\mathrm{sphere},i} \\cdot [1 + h_{1,i} 

521 \\,f_{\\mathrm{oblate}} 

522 \\cdot \\sin^2(F_i) 

523 \\cos^2(G_i) - h_{2,i}\\, f_{\\mathrm{oblate}} 

524 \\cdot \\cos^2(F_i) \\sin^2(G_i)] 

525 

526 :param a_lats: Latitudes in [deg] point A. 

527 :param a_lons: Longitudes in [deg] point A. 

528 :param b_lats: Latitudes in [deg] point B. 

529 :param b_lons: Longitudes in [deg] point B. 

530 :type a_lats: :py:class:`numpy.ndarray`, ``(N)`` 

531 :type a_lons: :py:class:`numpy.ndarray`, ``(N)`` 

532 :type b_lats: :py:class:`numpy.ndarray`, ``(N)`` 

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

534 

535 :return: Distances in [m]. 

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

537 ''' 

538 

539 a_lats, a_lons, b_lats, b_lons = float_array_broadcast( 

540 a_lats, a_lons, b_lats, b_lons) 

541 

542 assert implementation in ('c', 'python') 

543 if implementation == 'c': 

544 from pyrocko import orthodrome_ext 

545 return orthodrome_ext.distance_accurate50m_numpy( 

546 a_lats, a_lons, b_lats, b_lons) 

547 

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

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

550 

551 if num.all(eq): 

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

553 

554 def extr(x): 

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

556 return x[ii_neq] 

557 else: 

558 return x 

559 

560 a_lats = extr(a_lats) 

561 a_lons = extr(a_lons) 

562 b_lats = extr(b_lats) 

563 b_lons = extr(b_lons) 

564 

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

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

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

568 

569 s = num.sin(g)**2 * num.cos(h)**2 + num.cos(f)**2 * num.sin(h)**2 

570 c = num.cos(g)**2 * num.cos(h)**2 + num.sin(f)**2 * num.sin(h)**2 

571 

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

573 

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

575 

576 d = 2.*w*earthradius_equator 

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

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

579 

580 dists = num.zeros(eq.size, dtype=float) 

581 dists[ii_neq] = d * ( 

582 1. + 

583 earth_oblateness * h1 * num.sin(f)**2 * num.cos(g)**2 - 

584 earth_oblateness * h2 * num.cos(f)**2 * num.sin(g)**2) 

585 

586 return dists 

587 

588 

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

590 ''' 

591 Transform local cartesian coordinates to latitude and longitude. 

592 

593 From east and north coordinates (``x`` and ``y`` coordinate 

594 :py:class:`numpy.ndarray`) relative to a reference differences in 

595 longitude and latitude are calculated, which are effectively changes in 

596 azimuth and distance, respectively: 

597 

598 .. math:: 

599 

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

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

602 

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

604 / \\bf{y}). 

605 

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

607 

608 :param lat0: Latitude origin of the cartesian coordinate system in [deg]. 

609 :param lon0: Longitude origin of the cartesian coordinate system in [deg]. 

610 :param north_m: Northing distances from origin in [m]. 

611 :param east_m: Easting distances from origin in [m]. 

612 :type north_m: :py:class:`numpy.ndarray`, ``(N)`` 

613 :type east_m: :py:class:`numpy.ndarray`, ``(N)`` 

614 :type lat0: float 

615 :type lon0: float 

616 

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

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

619 

620 ''' 

621 

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

623 gamma = num.arctan2(east_m, north_m) 

624 

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

626 

627 

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

629 ''' 

630 Absolute latitudes and longitudes are calculated from relative changes. 

631 

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

633 distance given in degrees. 

634 

635 :param lat0: Latitude origin of the cartesian coordinate system in [deg]. 

636 :type lat0: float 

637 :param lon0: Longitude origin of the cartesian coordinate system in [deg]. 

638 :type lon0: float 

639 :param azimuth_deg: Azimuth from origin in [deg]. 

640 :type azimuth_deg: :py:class:`numpy.ndarray`, ``(N)`` 

641 :param distance_deg: Distances from origin in [deg]. 

642 :type distance_deg: :py:class:`numpy.ndarray`, ``(N)`` 

643 

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

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

646 ''' 

647 

648 return azidist_to_latlon_rad( 

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

650 

651 

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

653 ''' 

654 Absolute latitudes and longitudes are calculated from relative changes. 

655 

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

657 enforced for ``c`` and ``alpha``. 

658 

659 .. math:: 

660 

661 \\Delta {\\bf a}_i \\; \\text{and} \\; \\Delta \\gamma_i \\; 

662 \\text{are relative distances and azimuths from lat0 and lon0 for 

663 \\textit{i} source points of a finite source.} 

664 

665 .. math:: 

666 

667 \\mathrm{b} &= \\frac{\\pi}{2} -\\frac{\\pi}{180}\\;\\mathrm{lat_0}\\\\ 

668 {\\bf c}_i &=\\arccos[\\; \\cos(\\Delta {\\bf{a}}_i) 

669 \\cos(\\mathrm{b}) + |\\Delta \\gamma_i| \\, 

670 \\sin(\\Delta {\\bf a}_i) 

671 \\sin(\\mathrm{b})\\; ] \\\\ 

672 \\mathrm{lat}_i &= \\frac{180}{\\pi} 

673 \\left(\\frac{\\pi}{2} - {\\bf c}_i \\right) 

674 

675 .. math:: 

676 

677 \\alpha_i &= \\arcsin \\left[ \\; \\frac{ \\sin(\\Delta {\\bf a}_i ) 

678 \\sin(|\\Delta \\gamma_i|)}{\\sin({\\bf c}_i)}\\; 

679 \\right] \\\\ 

680 \\alpha_i &= \\begin{cases} 

681 \\alpha_i, &\\text{if} \\; \\cos(\\Delta {\\bf a}_i) - 

682 \\cos(\\mathrm{b}) \\cos({\\bf{c}}_i) > 0, \\; 

683 \\text{else} \\\\ 

684 \\pi - \\alpha_i, & \\text{if} \\; \\alpha_i > 0,\\; 

685 \\text{else}\\\\ 

686 -\\pi - \\alpha_i, & \\text{if} \\; \\alpha_i < 0. 

687 \\end{cases} \\\\ 

688 \\mathrm{lon}_i &= \\mathrm{lon_0} + 

689 \\frac{180}{\\pi} \\, 

690 \\frac{\\Delta \\gamma_i }{|\\Delta \\gamma_i|} 

691 \\cdot \\alpha_i 

692 \\text{, with $\\alpha_i \\in [-\\pi,\\pi]$} 

693 

694 :param lat0: Latitude origin of the cartesian coordinate system in [deg]. 

695 :param lon0: Longitude origin of the cartesian coordinate system in [deg]. 

696 :param distance_rad: Distances from origin in [rad]. 

697 :param azimuth_rad: Azimuth from origin in [rad]. 

698 :type distance_rad: :py:class:`numpy.ndarray`, ``(N)`` 

699 :type azimuth_rad: :py:class:`numpy.ndarray`, ``(N)`` 

700 :type lat0: float 

701 :type lon0: float 

702 

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

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

705 ''' 

706 

707 a = distance_rad 

708 gamma = azimuth_rad 

709 

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

711 

712 alphasign = 1. 

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

714 gamma = num.abs(gamma) 

715 

716 c = num.arccos(clip( 

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

718 

719 alpha = num.arcsin(clip( 

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

721 

722 alpha = num.where( 

723 num.cos(a)-num.cos(b)*num.cos(c) < 0, 

724 num.where(alpha > 0, math.pi-alpha, -math.pi-alpha), 

725 alpha) 

726 

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

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

729 

730 return lat, lon 

731 

732 

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

734 ''' 

735 Transform local cartesian coordinates to latitude and longitude. 

736 

737 Like :py:func:`pyrocko.orthodrome.ne_to_latlon`, 

738 but this method (implementation below), although it should be numerically 

739 more stable, suffers problems at points which are *across the pole* 

740 as seen from the cartesian origin. 

741 

742 .. math:: 

743 

744 \\text{distance change:}\\; \\Delta {{\\bf a}_i} &= 

745 \\sqrt{{\\bf{y}}^2_i + {\\bf{x}}^2_i }/ \\mathrm{R_E},\\\\ 

746 \\text{azimuth change:}\\; \\Delta {\\bf \\gamma}_i &= 

747 \\arctan( {\\bf x}_i {\\bf y}_i). \\\\ 

748 \\mathrm{b} &= 

749 \\frac{\\pi}{2} -\\frac{\\pi}{180} \\;\\mathrm{lat_0}\\\\ 

750 

751 .. math:: 

752 

753 {{\\bf z}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i - 

754 \\mathrm{b}}{2} \\right)} 

755 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\ 

756 {{\\bf n}_1}_i &= \\cos{\\left( \\frac{\\Delta {\\bf a}_i + 

757 \\mathrm{b}}{2} \\right)} 

758 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\ 

759 {{\\bf z}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i - 

760 \\mathrm{b}}{2} \\right)} 

761 \\cos {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\ 

762 {{\\bf n}_2}_i &= \\sin{\\left( \\frac{\\Delta {\\bf a}_i + 

763 \\mathrm{b}}{2} \\right)} 

764 \\sin {\\left( \\frac{|\\gamma_i|}{2} \\right) }\\\\ 

765 {{\\bf t}_1}_i &= \\arctan{\\left( \\frac{{{\\bf z}_1}_i} 

766 {{{\\bf n}_1}_i} \\right) }\\\\ 

767 {{\\bf t}_2}_i &= \\arctan{\\left( \\frac{{{\\bf z}_2}_i} 

768 {{{\\bf n}_2}_i} \\right) } \\\\[0.5cm] 

769 c &= \\begin{cases} 

770 2 \\cdot \\arccos \\left( {{\\bf z}_1}_i / \\sin({{\\bf t}_1}_i) 

771 \\right),\\; \\text{if } 

772 |\\sin({{\\bf t}_1}_i)| > 

773 |\\sin({{\\bf t}_2}_i)|,\\; \\text{else} \\\\ 

774 2 \\cdot \\arcsin{\\left( {{\\bf z}_2}_i / 

775 \\sin({{\\bf t}_2}_i) \\right)}. 

776 \\end{cases}\\\\ 

777 

778 .. math:: 

779 

780 {\\bf {lat}}_i &= \\frac{180}{ \\pi } \\left( \\frac{\\pi}{2} 

781 - {\\bf {c}}_i \\right) \\\\ 

782 {\\bf {lon}}_i &= {\\bf {lon}}_0 + \\frac{180}{ \\pi } 

783 \\frac{\\gamma_i}{|\\gamma_i|}, 

784 \\text{ with}\\; \\gamma \\in [-\\pi,\\pi] 

785 

786 :param lat0: Latitude origin of the cartesian coordinate system in [deg]. 

787 :param lon0: Longitude origin of the cartesian coordinate system in [deg]. 

788 :param north_m: Northing distances from origin in [m]. 

789 :param east_m: Easting distances from origin in [m]. 

790 :type north_m: :py:class:`numpy.ndarray`, ``(N)`` 

791 :type east_m: :py:class:`numpy.ndarray`, ``(N)`` 

792 :type lat0: float 

793 :type lon0: float 

794 

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

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

797 ''' 

798 

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

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

801 

802 gamma = num.arctan2(east_m, north_m) 

803 alphasign = 1. 

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

805 gamma = num.abs(gamma) 

806 

807 z1 = num.cos((a-b)/2.)*num.cos(gamma/2.) 

808 n1 = num.cos((a+b)/2.)*num.sin(gamma/2.) 

809 z2 = num.sin((a-b)/2.)*num.cos(gamma/2.) 

810 n2 = num.sin((a+b)/2.)*num.sin(gamma/2.) 

811 t1 = num.arctan2(z1, n1) 

812 t2 = num.arctan2(z2, n2) 

813 

814 alpha = t1 + t2 

815 

816 sin_t1 = num.sin(t1) 

817 sin_t2 = num.sin(t2) 

818 c = num.where( 

819 num.abs(sin_t1) > num.abs(sin_t2), 

820 num.arccos(z1/sin_t1)*2., 

821 num.arcsin(z2/sin_t2)*2.) 

822 

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

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

825 return lat, lon 

826 

827 

828def latlon_to_ne(*args): 

829 ''' 

830 Relative cartesian coordinates with respect to a reference location. 

831 

832 For two locations, a reference location A and another location B, given in 

833 geographical coordinates in degrees, the corresponding cartesian 

834 coordinates are calculated. 

835 Assisting functions are :py:func:`pyrocko.orthodrome.azimuth` and 

836 :py:func:`pyrocko.orthodrome.distance_accurate50m`. 

837 

838 .. math:: 

839 

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

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

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

843 

844 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180} 

845 \\varphi_{\\mathrm{azi},AB} )\\\\ 

846 e &= D_{AB} \\cdot 

847 \\sin( \\frac{\\pi }{180} \\varphi_{\\mathrm{azi},AB}) 

848 

849 :param refloc: Location reference point. 

850 :type refloc: :py:class:`pyrocko.orthodrome.Loc` 

851 :param loc: Location of interest. 

852 :type loc: :py:class:`pyrocko.orthodrome.Loc` 

853 

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

855 :rtype: tuple[float, float] 

856 

857 ''' 

858 

859 azi = azimuth(*args) 

860 dist = distance_accurate50m(*args) 

861 n, e = math.cos(azi*d2r)*dist, math.sin(azi*d2r)*dist 

862 return n, e 

863 

864 

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

866 ''' 

867 Relative cartesian coordinates with respect to a reference location. 

868 

869 For two locations, a reference location (``lat0``, ``lon0``) and another 

870 location B, given in geographical coordinates in degrees, 

871 the corresponding cartesian coordinates are calculated. 

872 Assisting functions are :py:func:`azimuth` 

873 and :py:func:`distance_accurate50m`. 

874 

875 :param lat0: Latitude of the reference location in [deg]. 

876 :param lon0: Longitude of the reference location in [deg]. 

877 :param lat: Latitude of the absolute location in [deg]. 

878 :param lon: Longitude of the absolute location in [deg]. 

879 

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

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

882 

883 Implemented formulations: 

884 

885 .. math:: 

886 

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

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

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

890 

891 n &= D_{AB} \\cdot \\cos( \\frac{\\pi }{180} \\varphi_{ 

892 \\mathrm{azi},AB} )\\\\ 

893 e &= D_{AB} \\cdot \\sin( \\frac{\\pi }{180} \\varphi_{ 

894 \\mathrm{azi},AB} ) 

895 ''' 

896 

897 azi = azimuth_numpy(lat0, lon0, lat, lon) 

898 dist = distance_accurate50m_numpy(lat0, lon0, lat, lon) 

899 n = num.cos(azi*d2r)*dist 

900 e = num.sin(azi*d2r)*dist 

901 return n, e 

902 

903 

904_wgs84 = None 

905 

906 

907def get_wgs84(): 

908 global _wgs84 

909 if _wgs84 is None: 

910 from geographiclib.geodesic import Geodesic 

911 _wgs84 = Geodesic.WGS84 

912 

913 return _wgs84 

914 

915 

916def amap(n): 

917 def wrap(f): 

918 if n == 1: 

919 def func(*args): 

920 it = num.nditer(args + (None,)) 

921 for ops in it: 

922 ops[-1][...] = f(*ops[:-1]) 

923 

924 return it.operands[-1] 

925 elif n == 2: 

926 def func(*args): 

927 it = num.nditer(args + (None, None)) 

928 for ops in it: 

929 ops[-2][...], ops[-1][...] = f(*ops[:-2]) 

930 

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

932 

933 return func 

934 return wrap 

935 

936 

937@amap(2) 

938def ne_to_latlon2(lat0, lon0, north_m, east_m): 

939 wgs84 = get_wgs84() 

940 az = num.arctan2(east_m, north_m)*r2d 

941 dist = num.sqrt(east_m**2 + north_m**2) 

942 x = wgs84.Direct(lat0, lon0, az, dist) 

943 return x['lat2'], x['lon2'] 

944 

945 

946@amap(2) 

947def latlon_to_ne2(lat0, lon0, lat1, lon1): 

948 wgs84 = get_wgs84() 

949 x = wgs84.Inverse(lat0, lon0, lat1, lon1) 

950 dist = x['s12'] 

951 az = x['azi1'] 

952 n = num.cos(az*d2r)*dist 

953 e = num.sin(az*d2r)*dist 

954 return n, e 

955 

956 

957@amap(1) 

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

959 wgs84 = get_wgs84() 

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

961 

962 

963def positive_region(region): 

964 ''' 

965 Normalize parameterization of a rectangular geographical region. 

966 

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

968 :type region: :py:class:`tuple` of :py:class:`float` 

969 

970 :returns: ``(west, east, south, north)``, where ``west <= east`` and 

971 where ``west`` and ``east`` are in the range ``[-180., 180.+360.]``. 

972 :rtype: :py:class:`tuple` of :py:class:`float` 

973 ''' 

974 west, east, south, north = [float(x) for x in region] 

975 

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

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

978 assert -90. <= south < 90. 

979 assert -90. < north <= 90. 

980 

981 if east < west: 

982 east += 360. 

983 

984 if west < -180.: 

985 west += 360. 

986 east += 360. 

987 

988 return (west, east, south, north) 

989 

990 

991def points_in_region(p, region): 

992 ''' 

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

994 

995 :param p: ``(lat, lon)`` pairs in [deg]. 

996 :type p: :py:class:`numpy.ndarray` ``(N, 2)`` 

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

998 :type region: :py:class:`tuple` of :py:class:`float` 

999 

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

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

1002 ''' 

1003 

1004 w, e, s, n = positive_region(region) 

1005 return num.logical_and( 

1006 num.logical_and(s <= p[:, 0], p[:, 0] <= n), 

1007 num.logical_or( 

1008 num.logical_and(w <= p[:, 1], p[:, 1] <= e), 

1009 num.logical_and(w-360. <= p[:, 1], p[:, 1] <= e-360.))) 

1010 

1011 

1012def point_in_region(p, region): 

1013 ''' 

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

1015 

1016 :param p: ``(lat, lon)`` in [deg]. 

1017 :type p: :py:class:`tuple` of :py:class:`float` 

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

1019 :type region: :py:class:`tuple` of :py:class:`float` 

1020 

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

1022 :rtype: bool 

1023 ''' 

1024 

1025 w, e, s, n = positive_region(region) 

1026 return num.logical_and( 

1027 num.logical_and(s <= p[0], p[0] <= n), 

1028 num.logical_or( 

1029 num.logical_and(w <= p[1], p[1] <= e), 

1030 num.logical_and(w-360. <= p[1], p[1] <= e-360.))) 

1031 

1032 

1033def radius_to_region(lat, lon, radius): 

1034 ''' 

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

1036 

1037 :param lat: Latitude of the center point of circular region in [deg]. 

1038 :type lat: float 

1039 :param lon: Longitude of the center point of circular region in [deg]. 

1040 :type lon: float 

1041 :param radius: Radius of circular region in [m]. 

1042 :type radius: float 

1043 

1044 :returns: Rectangular region as ``(east, west, south, north)`` in [deg] or 

1045 ``None``. 

1046 :rtype: tuple[float, float, float, float] 

1047 ''' 

1048 radius_deg = radius * m2d 

1049 if radius_deg < 45.: 

1050 lat_min = max(-90., lat - radius_deg) 

1051 lat_max = min(90., lat + radius_deg) 

1052 absmaxlat = max(abs(lat_min), abs(lat_max)) 

1053 if absmaxlat > 89: 

1054 lon_min = -180. 

1055 lon_max = 180. 

1056 else: 

1057 lon_min = max( 

1058 -180. - 360., 

1059 lon - radius_deg / math.cos(absmaxlat*d2r)) 

1060 lon_max = min( 

1061 180. + 360., 

1062 lon + radius_deg / math.cos(absmaxlat*d2r)) 

1063 

1064 lon_min, lon_max, lat_min, lat_max = positive_region( 

1065 (lon_min, lon_max, lat_min, lat_max)) 

1066 

1067 return lon_min, lon_max, lat_min, lat_max 

1068 

1069 else: 

1070 return None 

1071 

1072 

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

1074 ''' 

1075 Calculate geographic midpoints by finding the center of gravity. 

1076 

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

1078 poles. 

1079 

1080 :param lats: Latitudes in [deg]. 

1081 :param lons: Longitudes in [deg]. 

1082 :param weights: Weighting factors. 

1083 :type lats: :py:class:`numpy.ndarray`, ``(N)`` 

1084 :type lons: :py:class:`numpy.ndarray`, ``(N)`` 

1085 :type weights: optional, :py:class:`numpy.ndarray`, ``(N)`` 

1086 

1087 :return: Latitudes and longitudes of the midpoints in [deg]. 

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

1089 ''' 

1090 if not weights: 

1091 weights = num.ones(len(lats)) 

1092 

1093 total_weigth = num.sum(weights) 

1094 weights /= total_weigth 

1095 lats = lats * d2r 

1096 lons = lons * d2r 

1097 x = num.sum(num.cos(lats) * num.cos(lons) * weights) 

1098 y = num.sum(num.cos(lats) * num.sin(lons) * weights) 

1099 z = num.sum(num.sin(lats) * weights) 

1100 

1101 lon = num.arctan2(y, x) 

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

1103 lat = num.arctan2(z, hyp) 

1104 

1105 return lat/d2r, lon/d2r 

1106 

1107 

1108def geographic_midpoint_locations(locations, weights=None): 

1109 coords = num.array([loc.effective_latlon 

1110 for loc in locations]) 

1111 return geographic_midpoint(coords[:, 0], coords[:, 1], weights) 

1112 

1113 

1114def geodetic_to_ecef(lat, lon, alt): 

1115 ''' 

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

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

1118 

1119 :param lat: Geodetic latitude in [deg]. 

1120 :param lon: Geodetic longitude in [deg]. 

1121 :param alt: Geodetic altitude (height) in [m] (positive for points outside 

1122 the geoid). 

1123 :type lat: float 

1124 :type lon: float 

1125 :type alt: float 

1126 

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

1128 :rtype: tuple[float, float, float] 

1129 

1130 .. [#1] https://en.wikipedia.org/wiki/ECEF 

1131 .. [#2] https://en.wikipedia.org/wiki/Geographic_coordinate_conversion 

1132 #From_geodetic_to_ECEF_coordinates 

1133 ''' 

1134 

1135 f = earth_oblateness 

1136 a = earthradius_equator 

1137 e2 = 2*f - f**2 

1138 

1139 lat, lon = num.radians(lat), num.radians(lon) 

1140 # Normal (plumb line) 

1141 N = a / num.sqrt(1.0 - (e2 * num.sin(lat)**2)) 

1142 

1143 X = (N+alt) * num.cos(lat) * num.cos(lon) 

1144 Y = (N+alt) * num.cos(lat) * num.sin(lon) 

1145 Z = (N*(1.0-e2) + alt) * num.sin(lat) 

1146 

1147 return (X, Y, Z) 

1148 

1149 

1150def ecef_to_geodetic(X, Y, Z): 

1151 ''' 

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

1153 geodetic coordinates (Ferrari's solution). 

1154 

1155 :param X: ECEF X coordinate [m]. 

1156 :type X: float 

1157 :param Y: ECEF Y coordinate [m]. 

1158 :type Y: float 

1159 :param Z: ECEF Z coordinate [m]. 

1160 :type Z: float 

1161 

1162 :return: Geodetic coordinates (lat, lon, alt). Latitude and longitude are 

1163 in [deg] and altitude is in [m] 

1164 (positive for points outside the geoid). 

1165 :rtype: :py:class:`tuple` of :py:class:`float` 

1166 

1167 .. seealso :: 

1168 https://en.wikipedia.org/wiki/Geographic_coordinate_conversion 

1169 #The_application_of_Ferrari.27s_solution 

1170 ''' 

1171 f = earth_oblateness 

1172 a = earthradius_equator 

1173 b = a * (1. - f) 

1174 e2 = 2.*f - f**2 

1175 

1176 # usefull 

1177 a2 = a**2 

1178 b2 = b**2 

1179 e4 = e2**2 

1180 X2 = X**2 

1181 Y2 = Y**2 

1182 Z2 = Z**2 

1183 

1184 r = num.sqrt(X2 + Y2) 

1185 r2 = r**2 

1186 

1187 e_prime2 = (a2 - b2)/b2 

1188 E2 = a2 - b2 

1189 F = 54. * b2 * Z2 

1190 G = r2 + (1.-e2)*Z2 - (e2*E2) 

1191 C = (e4 * F * r2) / (G**3) 

1192 S = cbrt(1. + C + num.sqrt(C**2 + 2.*C)) 

1193 P = F / (3. * (S + 1./S + 1.)**2 * G**2) 

1194 Q = num.sqrt(1. + (2.*e4*P)) 

1195 

1196 dum1 = -(P*e2*r) / (1.+Q) 

1197 dum2 = 0.5 * a2 * (1. + 1./Q) 

1198 dum3 = (P * (1.-e2) * Z2) / (Q * (1.+Q)) 

1199 dum4 = 0.5 * P * r2 

1200 r0 = dum1 + num.sqrt(dum2 - dum3 - dum4) 

1201 

1202 U = num.sqrt((r - e2*r0)**2 + Z2) 

1203 V = num.sqrt((r - e2*r0)**2 + (1.-e2)*Z2) 

1204 Z0 = (b2*Z) / (a*V) 

1205 

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

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

1208 lon = num.arctan2(Y, X) 

1209 

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

1211 

1212 

1213class Farside(Exception): 

1214 pass 

1215 

1216 

1217def latlon_to_xyz(latlons): 

1218 if latlons.ndim == 1: 

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

1220 

1221 points = num.zeros((latlons.shape[0], 3)) 

1222 lats = latlons[:, 0] 

1223 lons = latlons[:, 1] 

1224 points[:, 0] = num.cos(lats*d2r) * num.cos(lons*d2r) 

1225 points[:, 1] = num.cos(lats*d2r) * num.sin(lons*d2r) 

1226 points[:, 2] = num.sin(lats*d2r) 

1227 return points 

1228 

1229 

1230def xyz_to_latlon(xyz): 

1231 if xyz.ndim == 1: 

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

1233 

1234 latlons = num.zeros((xyz.shape[0], 2)) 

1235 latlons[:, 0] = num.arctan2( 

1236 xyz[:, 2], num.sqrt(xyz[:, 0]**2 + xyz[:, 1]**2)) * r2d 

1237 latlons[:, 1] = num.arctan2( 

1238 xyz[:, 1], xyz[:, 0]) * r2d 

1239 return latlons 

1240 

1241 

1242def rot_to_00(lat, lon): 

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

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

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

1246 

1247 

1248def distances3d(a, b): 

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

1250 

1251 

1252def circulation(points2): 

1253 return num.sum( 

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

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

1256 

1257 

1258def stereographic(points): 

1259 dists = distances3d(points[1:, :], points[:-1, :]) 

1260 if dists.size > 0: 

1261 maxdist = num.max(dists) 

1262 cutoff = maxdist**2 / 2. 

1263 else: 

1264 cutoff = 1.0e-5 

1265 

1266 points = points.copy() 

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

1268 raise Farside() 

1269 

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

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

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

1273 

1274 return points_out 

1275 

1276 

1277def stereographic_poly(points): 

1278 dists = distances3d(points[1:, :], points[:-1, :]) 

1279 if dists.size > 0: 

1280 maxdist = num.max(dists) 

1281 cutoff = maxdist**2 / 2. 

1282 else: 

1283 cutoff = 1.0e-5 

1284 

1285 points = points.copy() 

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

1287 raise Farside() 

1288 

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

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

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

1292 

1293 if circulation(points_out) >= 0: 

1294 raise Farside() 

1295 

1296 return points_out 

1297 

1298 

1299def gnomonic_x(points, cutoff=0.01): 

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

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

1302 raise Farside() 

1303 

1304 factor = 1.0 / points[:, 0] 

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

1306 return points_out 

1307 

1308 

1309def cneg(i, x): 

1310 if i == 1: 

1311 return x 

1312 else: 

1313 return num.logical_not(x) 

1314 

1315 

1316def contains_points(polygon, points): 

1317 ''' 

1318 Test which points are inside polygon on a sphere. 

1319 

1320 The inside of the polygon is defined as the area which is to the left hand 

1321 side of an observer walking the polygon line, points in order, on the 

1322 sphere. Lines between the polygon points are treated as great circle paths. 

1323 The polygon may be arbitrarily complex, as long as it does not have any 

1324 crossings or thin parts with zero width. The polygon may contain the poles 

1325 and is allowed to wrap around the sphere multiple times. 

1326 

1327 The algorithm works by consecutive cutting of the polygon into (almost) 

1328 hemispheres and subsequent Gnomonic projections to perform the 

1329 point-in-polygon tests on a 2D plane. 

1330 

1331 :param polygon: Point coordinates defining the polygon [deg]. 

1332 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index 

1333 0=lat, 1=lon 

1334 :param points: Coordinates of points to test [deg]. 

1335 :type points: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index 

1336 0=lat, 1=lon 

1337 

1338 :returns: Boolean mask array. 

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

1340 ''' 

1341 

1342 and_ = num.logical_and 

1343 

1344 points_xyz = latlon_to_xyz(points) 

1345 mask_x = 0. <= points_xyz[:, 0] 

1346 mask_y = 0. <= points_xyz[:, 1] 

1347 mask_z = 0. <= points_xyz[:, 2] 

1348 

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

1350 

1351 for ix in [-1, 1]: 

1352 for iy in [-1, 1]: 

1353 for iz in [-1, 1]: 

1354 mask = and_( 

1355 and_(cneg(ix, mask_x), cneg(iy, mask_y)), 

1356 cneg(iz, mask_z)) 

1357 

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

1359 

1360 lat, lon = xyz_to_latlon(center_xyz) 

1361 rot = rot_to_00(lat, lon) 

1362 

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

1364 points_rot_pro = gnomonic_x(points_rot_xyz) 

1365 

1366 offset = 0.01 

1367 

1368 poly_xyz = latlon_to_xyz(polygon) 

1369 poly_rot_xyz = num.dot(rot, poly_xyz.T).T 

1370 poly_rot_xyz[:, 0] -= offset 

1371 groups = spoly_cut([poly_rot_xyz], axis=0) 

1372 

1373 for poly_rot_group_xyz in groups[1]: 

1374 poly_rot_group_xyz[:, 0] += offset 

1375 

1376 poly_rot_group_pro = gnomonic_x( 

1377 poly_rot_group_xyz) 

1378 

1379 if circulation(poly_rot_group_pro) > 0: 

1380 result[mask] += path_contains_points( 

1381 poly_rot_group_pro, points_rot_pro) 

1382 else: 

1383 result[mask] -= path_contains_points( 

1384 poly_rot_group_pro, points_rot_pro) 

1385 

1386 return result.astype(bool) 

1387 

1388 

1389def contains_point(polygon, point): 

1390 ''' 

1391 Test if point is inside polygon on a sphere. 

1392 

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

1394 

1395 :param polygon: Point coordinates defining the polygon [deg]. 

1396 :type polygon: :py:class:`numpy.ndarray` of shape ``(N, 2)``, second index 

1397 0=lat, 1=lon 

1398 :param point: Coordinates ``(lat, lon)`` of point to test [deg]. 

1399 :type point: :py:class:`tuple` of :py:class:`float` 

1400 

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

1402 :rtype: bool 

1403 ''' 

1404 

1405 return bool( 

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