Coverage for pygeodesy/azimuthal.py: 98%

1 

2# -*- coding: utf-8 -*- 

3 

4u'''Equidistant, Equal-Area, and other Azimuthal projections. 

5 

6Classes L{Equidistant}, L{EquidistantExact}, L{EquidistantGeodSolve}, 

7L{EquidistantKarney}, L{Gnomonic}, L{GnomonicExact}, L{GnomonicKarney}, 

8L{LambertEqualArea}, L{Orthographic} and L{Stereographic}, classes 

9L{AzimuthalError}, L{Azimuthal7Tuple} and functions L{equidistant} 

10and L{gnomonic}. 

11 

12L{EquidistantExact} and L{GnomonicExact} are based on exact geodesic classes 

13L{GeodesicExact} and L{GeodesicLineExact}, Python versions of I{Charles Karney}'s 

14C++ original U{GeodesicExact<https://GeographicLib.SourceForge.io/C++/doc/ 

15classGeographicLib_1_1GeodesicExact.html>}, respectively U{GeodesicLineExact 

16<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicLineExact.html>}. 

17 

18Using L{EquidistantGeodSolve} requires I{Karney}'s utility U{GeodSolve 

19<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>} to be 

20executable and set in env variable C{PYGEODESY_GEODSOLVE}, see module 

21L{geodsolve} for more details. 

22 

23L{EquidistantKarney} and L{GnomonicKarney} require I{Karney}'s Python package 

24U{geographiclib<https://PyPI.org/project/geographiclib>} to be installed. 

25 

26Other azimuthal classes implement only (***) U{Snyder's FORMULAS FOR THE SPHERE 

27<https://Pubs.USGS.gov/pp/1395/report.pdf>} and use those for any datum, 

28spherical and ellipsoidal. The radius used for the latter is the ellipsoid's 

29I{mean radius of curvature} at the latitude of the projection center point. For 

30further justification, see the first paragraph under U{Snyder's FORMULAS FOR THE 

31ELLIPSOID, page 197<https://Pubs.USGS.gov/pp/1395/report.pdf>}. 

32 

33Page numbers in C{Snyder} references apply to U{John P. Snyder, "Map Projections 

34-- A Working Manual", 1987<https://Pubs.USGS.gov/pp/1395/report.pdf>}. 

35 

36See also U{here<https://WikiPedia.org/wiki/Azimuthal_equidistant_projection>}, 

37especially the U{Comparison of the Azimuthal equidistant projection and some 

38azimuthal projections centred on 90° N at the same scale, ordered by projection 

39altitude in Earth radii<https://WikiPedia.org/wiki/Azimuthal_equidistant_projection 

40#/media/File:Comparison_azimuthal_projections.svg>}. 

41''' 

42# make sure int/int division yields float quotient, see .basics 

43from __future__ import division as _; del _ # noqa: E702 ; 

44 

45# from pygeodesy.basics import _isin, _xinstanceof # from .ellipsoidalBase 

46from pygeodesy.constants import EPS, EPS0, EPS1, NAN, isnon0, _umod_360, \ 

47 _EPStol, _0_0, _0_1, _0_5, _1_0, _N_1_0, _2_0 

48from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB, \ 

49 _isin, _xinstanceof 

50from pygeodesy.datums import _spherical_datum, _WGS84 

51from pygeodesy.errors import _ValueError, _xattr, _xdatum, _xkwds 

52from pygeodesy.fmath import euclid, hypot as _hypot, Fsum 

53# from pygeodesy.fsums import Fsum # from .fmath 

54# from pygeodesy.formy import antipode # _MODS 

55# from pygeodesy.internals import typename # from .karney 

56from pygeodesy.interns import _azimuth_, _datum_, _lat_, _lon_, _scale_, \ 

57 _SPACE_, _x_, _y_ 

58from pygeodesy.karney import _norm180, typename 

59from pygeodesy.latlonBase import _MODS, LatLonBase as _LLB 

60from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS # ALL_MODS 

61from pygeodesy.named import _name__, _name2__, _NamedBase, _NamedTuple, _Pass 

62from pygeodesy.namedTuples import LatLon2Tuple, LatLon4Tuple 

63from pygeodesy.props import deprecated_Property_RO, Property_RO, \ 

64 property_doc_, _update_all 

65from pygeodesy.streprs import Fmt, _fstrLL0, unstr 

66from pygeodesy.units import Azimuth, Easting, Lat_, Lon_, Northing, \ 

67 Scalar, Scalar_ 

68from pygeodesy.utily import asin1, atan1, atan2, atan2b, atan2d, \ 

69 sincos2, sincos2d, sincos2d_ 

70 

71from math import acos, degrees, fabs, sin, sqrt 

72 

73__all__ = _ALL_LAZY.azimuthal 

74__version__ = '25.11.29' 

75 

76_EPS_K = _EPStol * _0_1 # Karney's eps_ or _EPSmin * _0_1? 

77_over_horizon_ = 'over horizon' 

78_TRIPS = 21 # numit, 4 sufficient 

79 

80 

81def _enzh4(x, y, *h): 

82 '''(INTERNAL) Return 4-tuple (easting, northing, azimuth, hypot). 

83 ''' 

84 e = Easting( x=x) 

85 n = Northing(y=y) 

86 z = atan2b(e, n) # (x, y) for azimuth from true North 

87 return e, n, z, (h[0] if h else _hypot(e, n)) 

88 

89 

90class _AzimuthalBase(_NamedBase): 

91 '''(INTERNAL) Base class for azimuthal projections. 

92 

93 @see: I{Karney}'s C++ class U{AzimuthalEquidistant<https://GeographicLib.SourceForge.io/ 

94 C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>} and U{Gnomonic 

95 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>} or 

96 the C{PyGeodesy} versions thereof L{EquidistantKarney} respectively L{GnomonicKarney}. 

97 ''' 

98 _datum = _WGS84 # L{Datum} 

99 _latlon0 = LatLon2Tuple(_0_0, _0_0) # lat0, lon0 (L{LatLon2Tuple}) 

100 _sc0 = _0_0, _1_0 # 2-Tuple C{sincos2d(lat0)} 

101 

102 def __init__(self, lat0, lon0, datum=None, **name): 

103 '''New azimuthal projection. 

104 

105 @arg lat0: Latitude of the center point (C{degrees90}). 

106 @arg lon0: Longitude of the center point (C{degrees180}). 

107 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid}, 

108 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth 

109 radius (C{meter}). 

110 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}). 

111 

112 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or (spherical) B{C{datum}}. 

113 

114 @raise TypeError: Invalid B{C{datum}}. 

115 ''' 

116 if not _isin(datum, None, self._datum): 

117 self._datum = _spherical_datum(datum, **name) 

118 if name: 

119 self.name = name 

120 

121 if lat0 or lon0: # often both 0 

122 self._reset(lat0, lon0) 

123 

124 @Property_RO 

125 def datum(self): 

126 '''Get the datum (L{Datum}). 

127 ''' 

128 return self._datum 

129 

130 @Property_RO 

131 def equatoradius(self): 

132 '''Get the geodesic's equatorial radius, semi-axis (C{meter}). 

133 ''' 

134 return self.datum.ellipsoid.a 

135 

136 a = equatoradius 

137 

138 @Property_RO 

139 def flattening(self): 

140 '''Get the geodesic's flattening (C{scalar}). 

141 ''' 

142 return self.datum.ellipsoid.f 

143 

144 f = flattening 

145 

146 def forward(self, lat, lon, **name): # PYCHOK no cover 

147 '''I{Must be overloaded}.''' 

148 self._notOverloaded(lat, lon, **name) 

149 

150 def _forward(self, lat, lon, name, _k_t): 

151 '''(INTERNAL) Azimuthal (spherical) forward C{lat, lon} to C{x, y}. 

152 ''' 

153 lat, lon = Lat_(lat), Lon_(lon) 

154 sa, ca, sb, cb = sincos2d_(lat, lon - self.lon0) 

155 s0, c0 = self._sc0 

156 

157 cb *= ca 

158 k, t = _k_t(c0 * cb + s0 * sa) 

159 if t: 

160 r = k * self.radius 

161 y = r * (c0 * sa - s0 * cb) 

162 e, n, z, _ = _enzh4(r * sb * ca, y, None) 

163 else: # 0 or 180 

164 e = n = z = _0_0 

165 

166 t = Azimuthal7Tuple(e, n, lat, lon, z, k, self.datum, 

167 name=self._name__(name)) 

168 return t 

169 

170 def _forwards(self, *lls): 

171 '''(INTERNAL) One or more C{.forward} calls, see .ellipsoidalBaseDI. 

172 ''' 

173 _fwd = self.forward 

174 for ll in lls: 

175 yield _fwd(ll.lat, ll.lon) 

176 

177 @Property_RO 

178 def lat0(self): 

179 '''Get the center latitude (C{degrees90}). 

180 ''' 

181 return self._latlon0.lat 

182 

183 @property 

184 def latlon0(self): 

185 '''Get the center lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}) in (C{degrees90}, C{degrees180}). 

186 ''' 

187 return self._latlon0 

188 

189 @latlon0.setter # PYCHOK setter! 

190 def latlon0(self, latlon0): 

191 '''Set the center lat- and longitude (C{LatLon}, L{LatLon2Tuple} or L{LatLon4Tuple}). 

192 

193 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or ellipsoidal mismatch 

194 of B{C{latlon0}} and this projection. 

195 ''' 

196 B = _LLEB if self.datum.isEllipsoidal else _LLB 

197 _xinstanceof(B, LatLon2Tuple, LatLon4Tuple, latlon0=latlon0) 

198 D = _xattr(latlon0, datum=B) 

199 if D is not B: 

200 _xdatum(self.datum, D, Error=AzimuthalError) 

201 self.reset(latlon0.lat, latlon0.lon) 

202 

203 @Property_RO 

204 def lon0(self): 

205 '''Get the center longitude (C{degrees180}). 

206 ''' 

207 return self._latlon0.lon 

208 

209 @deprecated_Property_RO 

210 def majoradius(self): # PYCHOK no cover 

211 '''DEPRECATED, use property C{equatoradius}.''' 

212 return self.equatoradius 

213 

214 @Property_RO 

215 def radius(self): 

216 '''Get this projection's mean radius of curvature (C{meter}). 

217 ''' 

218 return self.datum.ellipsoid.rocMean(self.lat0) 

219 

220 def reset(self, lat0, lon0): 

221 '''Set or reset the center point of this azimuthal projection. 

222 

223 @arg lat0: Center point latitude (C{degrees90}). 

224 @arg lon0: Center point longitude (C{degrees180}). 

225 

226 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}. 

227 ''' 

228 _update_all(self) # zap caches 

229 self._reset(lat0, lon0) 

230 

231 def _reset(self, lat0, lon0): 

232 '''(INTERNAL) Update the center point. 

233 ''' 

234 self._latlon0 = LatLon2Tuple(Lat_(lat0=lat0, Error=AzimuthalError), 

235 Lon_(lon0=lon0, Error=AzimuthalError)) 

236 self._sc0 = sincos2d(self.lat0) 

237 

238 def reverse(self, x, y, **name_LatLon_and_kwds): 

239 '''I{Must be overloaded}.''' 

240 self._notOverloaded(x, y, **name_LatLon_and_kwds) # PYCHOK no cover 

241 

242 def _reverse(self, x, y, _c, lea, LatLon=None, **name_LatLon_kwds): 

243 '''(INTERNAL) Azimuthal (spherical) reverse C{x, y} to C{lat, lon}. 

244 ''' 

245 e, n, z, r = _enzh4(x, y) 

246 

247 c = _c(r / self.radius) 

248 if c is None: 

249 lat, lon = self.latlon0 

250 k, z = _1_0, _0_0 

251 else: 

252 s0, c0 = self._sc0 

253 sc, cc = sincos2(c) 

254 k = c / sc 

255 s = s0 * cc 

256 if r > EPS0: 

257 s += c0 * sc * (n / r) 

258 lat = degrees(asin1(s)) 

259 if lea or fabs(c0) > EPS: 

260 d = atan2d(e * sc, r * c0 * cc - n * s0 * sc) 

261 else: 

262 d = atan2d(e, (n if s0 < 0 else -n)) 

263 lon = _norm180(self.lon0 + d) 

264 

265 if LatLon is None: 

266 t, _ = _name2__(name_LatLon_kwds, _or_nameof=self) 

267 t = Azimuthal7Tuple(e, n, lat, lon, z, k, self.datum, name=t) 

268 else: 

269 t = self._toLatLon(lat, lon, LatLon, name_LatLon_kwds) 

270 return t 

271 

272 def _reverse2(self, x_t, *y): 

273 '''(INTERNAL) See iterating functions .ellipsoidalBaseDI._intersect3, 

274 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne. 

275 ''' 

276 t = self.reverse(x_t, *y) if y else self.reverse(x_t.x, x_t.y) # LatLon=None 

277 d = euclid(t.lat - self.lat0, t.lon - self.lon0) # degrees 

278 return t, d 

279 

280 def _toLatLon(self, lat, lon, LatLon, name_LatLon_kwds): 

281 '''(INTERNAL) Check B{C{LatLon}} and return an instance. 

282 ''' 

283 kwds = _xkwds(name_LatLon_kwds, datum=self.datum, _or_nameof=self) 

284 r = LatLon(lat, lon, **kwds) # handle .classof 

285 B = _LLEB if self.datum.isEllipsoidal else _LLB 

286 _xinstanceof(B, LatLon=r) 

287 return r 

288 

289 def toRepr(self, prec=6, **unused): # PYCHOK expected 

290 '''Return a string representation of this projection. 

291 

292 @kwarg prec: Number of (decimal) digits, unstripped (C{int}). 

293 

294 @return: This projection as C{"<classname>(lat0, lon0, ...)"} 

295 (C{str}). 

296 ''' 

297 return _fstrLL0(self, prec, True) 

298 

299 def toStr(self, prec=6, sep=_SPACE_, **unused): # PYCHOK expected 

300 '''Return a string representation of this projection. 

301 

302 @kwarg prec: Number of (decimal) digits, unstripped (C{int}). 

303 @kwarg sep: Separator to join (C{str}). 

304 

305 @return: This projection as C{"lat0 lon0"} (C{str}). 

306 ''' 

307 t = _fstrLL0(self, prec, False) 

308 return t if sep is None else sep.join(t) 

309 

310 

311class AzimuthalError(_ValueError): 

312 '''An azimuthal L{Equidistant}, L{EquidistantKarney}, L{Gnomonic}, 

313 L{LambertEqualArea}, L{Orthographic}, L{Stereographic} or 

314 L{Azimuthal7Tuple} issue. 

315 ''' 

316 pass 

317 

318 

319class Azimuthal7Tuple(_NamedTuple): 

320 '''7-Tuple C{(x, y, lat, lon, azimuth, scale, datum)}, in C{meter}, C{meter}, 

321 C{degrees90}, C{degrees180}, compass C{degrees}, C{scalar} and C{Datum} 

322 where C{(x, y)} is the easting and northing of a projected point, C{(lat, 

323 lon)} the geodetic location, C{azimuth} the azimuth, clockwise from true 

324 North and C{scale} is the projection scale, either C{1 / reciprocal} or 

325 C{1} or C{-1} in the L{Equidistant} case. 

326 ''' 

327 _Names_ = (_x_, _y_, _lat_, _lon_, _azimuth_, _scale_, _datum_) 

328 _Units_ = ( Easting, Northing, Lat_, Lon_, Azimuth, Scalar, _Pass) 

329 

330 def antipodal(self, azimuth=None): 

331 '''Return this tuple with the antipodal C{lat} and C{lon}. 

332 

333 @kwarg azimuth: Optional azimuth, overriding the current azimuth 

334 (C{compass degrees360}). 

335 ''' 

336 a = _MODS.formy.antipode(self.lat, self.lon) # PYCHOK named 

337 z = self.azimuth if azimuth is None else Azimuth(azimuth) # PYCHOK named 

338 return _NamedTuple.dup(self, lat=a.lat, lon=a.lon, azimuth=z) 

339 

340 

341class Equidistant(_AzimuthalBase): 

342 '''Azimuthal equidistant projection for the sphere***, see U{Snyder, pp 195-197 

343 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram 

344 <https://MathWorld.Wolfram.com/AzimuthalEquidistantProjection.html>}. 

345 

346 @note: Results from this L{Equidistant} and an L{EquidistantExact}, 

347 L{EquidistantGeodSolve} or L{EquidistantKarney} projection 

348 C{may differ} by 10% or more. For an example, see method 

349 C{testDiscrepancies} in module C{testAzimuthal.py}. 

350 ''' 

351 if _FOR_DOCS: 

352 __init__ = _AzimuthalBase.__init__ 

353 

354 def forward(self, lat, lon, **name): 

355 '''Convert a geodetic location to azimuthal equidistant east- and northing. 

356 

357 @arg lat: Latitude of the location (C{degrees90}). 

358 @arg lon: Longitude of the location (C{degrees180}). 

359 @kwarg name: Optional C{B{name}=NN} for the location (C{str}). 

360 

361 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)} 

362 with easting C{x} and northing C{y} of point in C{meter} and C{lat} 

363 and C{lon} in C{degrees} and C{azimuth} clockwise from true North. 

364 The C{scale} of the projection is C{1} in I{radial} direction and 

365 is C{1 / reciprocal} in the direction perpendicular to this. 

366 

367 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}. 

368 

369 @note: The C{scale} will be C{-1} if B{C{(lat, lon)}} is antipodal to the 

370 projection center C{(lat0, lon0)}. 

371 ''' 

372 def _k_t(c): 

373 k = _N_1_0 if c < 0 else _1_0 

374 t = fabs(c) < EPS1 

375 if t: 

376 a = acos(c) 

377 s = sin(a) 

378 if s: 

379 k = a / s 

380 return k, t 

381 

382 return self._forward(lat, lon, name, _k_t) 

383 

384 def reverse(self, x, y, **name_LatLon_and_kwds): 

385 '''Convert an azimuthal equidistant location to geodetic lat- and longitude. 

386 

387 @arg x: Easting of the location (C{meter}). 

388 @arg y: Northing of the location (C{meter}). 

389 @kwarg name_LatLon_and_kwds: Optional C{B{name}=NN} and class C{B{LatLon}=None} 

390 to use and optionally, additional B{C{LatLon}} keyword arguments, 

391 ignored if C{B{LatLon} is None}. 

392 

393 @return: The geodetic (C{LatLon}) or if C{B{LatLon} is None} an 

394 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}. 

395 

396 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon} 

397 in the range C{[-180..180] degrees}. The C{scale} of the 

398 projection is C{1} in I{radial} direction, C{azimuth} clockwise 

399 from true North and is C{1 / reciprocal} in the direction 

400 perpendicular to this. 

401 ''' 

402 def _c(c): 

403 return c if c > EPS else None 

404 

405 return self._reverse(x, y, _c, False, **name_LatLon_and_kwds) 

406 

407 

408def equidistant(lat0, lon0, datum=_WGS84, exact=False, geodsolve=False, **name): 

409 '''Return an L{EquidistantExact}, L{EquidistantGeodSolve} or (if I{Karney}'s 

410 U{geographiclib<https://PyPI.org/project/geographiclib>} package is 

411 installed) an L{EquidistantKarney}, otherwise an L{Equidistant} instance. 

412 

413 @arg lat0: Latitude of center point (C{degrees90}). 

414 @arg lon0: Longitude of center point (C{degrees180}). 

415 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid}, 

416 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth 

417 radius (C{meter}). 

418 @kwarg exact: Return an L{EquidistantExact} instance. 

419 @kwarg geodsolve: Return an L{EquidistantGeodSolve} instance. 

420 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}). 

421 

422 @return: An L{EquidistantExact}, L{EquidistantGeodSolve}, 

423 L{EquidistantKarney} or L{Equidistant} instance. 

424 

425 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or (spherical) B{C{datum}}. 

426 

427 @raise GeodesicError: Issue with L{GeodesicExact}, L{GeodesicSolve} 

428 or I{Karney}'s wrapped C{Geodesic}. 

429 

430 @raise TypeError: Invalid B{C{datum}}. 

431 ''' 

432 

433 E = EquidistantExact if exact else (EquidistantGeodSolve if geodsolve else Equidistant) 

434 if E is Equidistant: 

435 try: 

436 return EquidistantKarney(lat0, lon0, datum=datum, **name) # PYCHOK types 

437 except ImportError: 

438 pass 

439 return E(lat0, lon0, datum=datum, **name) # PYCHOK types 

440 

441 

442class _AzimuthalGeodesic(_AzimuthalBase): 

443 '''(INTERNAL) Base class for azimuthal projections using the 

444 I{wrapped} U{geodesic.Geodesic and geodesicline.GeodesicLine 

445 <https://GeographicLib.SourceForge.io/Python/doc/code.html>} or the 

446 I{exact} geodesic classes L{GeodesicExact} and L{GeodesicLineExact}. 

447 ''' 

448 _mask = 0 

449 

450 @Property_RO 

451 def geodesic(self): # PYCHOK no cover 

452 '''I{Must be overloaded}.''' 

453 self._notOverloaded() 

454 

455 def _7Tuple(self, e, n, r, name_LatLon_kwds, M=None): 

456 '''(INTERNAL) Return an C{Azimuthal7Tuple}. 

457 ''' 

458 s = M 

459 if s is None: # reciprocal, azimuthal scale 

460 s = (r.m12 / r.s12) if r.a12 > _EPS_K else _1_0 

461 z = _umod_360(r.azi2) # -180 <= r.azi2 < 180 ... 0 <= z < 360 

462 t, _ = _name2__(name_LatLon_kwds, _or_nameof=self) 

463 return Azimuthal7Tuple(e, n, r.lat2, r.lon2, z, s, self.datum, name=t) 

464 

465 

466class _EquidistantBase(_AzimuthalGeodesic): 

467 '''(INTERNAL) Base for classes L{EquidistantExact}, L{EquidistantGeodSolve} 

468 and L{EquidistantKarney}. 

469 ''' 

470 def __init__(self, lat0, lon0, datum=_WGS84, **name): 

471 '''New azimuthal L{EquidistantExact}, L{EquidistantGeodSolve} or 

472 L{EquidistantKarney} projection. 

473 ''' 

474 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, **name) 

475 

476 g = self.geodesic 

477 # g.STANDARD = g.AZIMUTH | g.DISTANCE | g.LATITUDE | g.LONGITUDE 

478 self._mask = g.REDUCEDLENGTH | g.STANDARD # | g.LONG_UNROLL 

479 

480 def forward(self, lat, lon, **name): 

481 '''Convert an (ellipsoidal) geodetic location to azimuthal equidistant east- and northing. 

482 

483 @arg lat: Latitude of the location (C{degrees90}). 

484 @arg lon: Longitude of the location (C{degrees180}). 

485 @kwarg name: Optional C{B{name}=NN} for the location (C{str}). 

486 

487 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)} 

488 with easting C{x} and northing C{y} of point in C{meter} and C{lat} 

489 and C{lon} in C{degrees} and C{azimuth} clockwise from true North. 

490 The C{scale} of the projection is C{1} in I{radial} direction and 

491 is C{1 / reciprocal} in the direction perpendicular to this. 

492 

493 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}. 

494 

495 @note: A call to C{.forward} followed by a call to C{.reverse} will return 

496 the original C{lat, lon} to within roundoff. 

497 ''' 

498 r = self.geodesic.Inverse(self.lat0, self.lon0, 

499 Lat_(lat), Lon_(lon), outmask=self._mask) 

500 x, y = sincos2d(r.azi1) 

501 return self._7Tuple(x * r.s12, y * r.s12, r, _name__(name)) 

502 

503 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds): # PYCHOK signature 

504 '''Convert an azimuthal equidistant location to (ellipsoidal) geodetic lat- and longitude. 

505 

506 @arg x: Easting of the location (C{meter}). 

507 @arg y: Northing of the location (C{meter}). 

508 @kwarg LatLon: Class to use (C{LatLon}) or C{None}. 

509 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} and optionally, additional 

510 B{C{LatLon}} keyword arguments, ignored if C{B{LatLon} is None}. 

511 

512 @return: The geodetic (C{LatLon}) or if C{B{LatLon} is None} an 

513 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}. 

514 

515 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon} 

516 in the range C{[-180..180] degrees}. The scale of the projection 

517 is C{1} in I{radial} direction, C{azimuth} clockwise from true 

518 North and is C{1 / reciprocal} in the direction perpendicular 

519 to this. 

520 ''' 

521 e, n, z, s = _enzh4(x, y) 

522 

523 r = self.geodesic.Direct(self.lat0, self.lon0, z, s, outmask=self._mask) 

524 return self._7Tuple(e, n, r, name_LatLon_kwds) if LatLon is None else \ 

525 self._toLatLon(r.lat2, r.lon2, LatLon, name_LatLon_kwds) 

526 

527 

528class EquidistantExact(_EquidistantBase): 

529 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant 

530 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>}, 

531 based on exact geodesic classes L{GeodesicExact} and L{GeodesicLineExact}. 

532 

533 An azimuthal equidistant projection is centered at an arbitrary position on the ellipsoid. 

534 For a point in projected space C{(x, y)}, the geodesic distance from the center position 

535 is C{hypot(x, y)} and the C{azimuth} of the geodesic from the center point is C{atan2(x, y)}, 

536 clockwise from true North. 

537 

538 The C{.forward} and C{.reverse} methods also return the C{azimuth} of the geodesic at C{(x, 

539 y)} and the C{scale} in the azimuthal direction which, together with the basic properties 

540 of the projection, serve to specify completely the local affine transformation between 

541 geographic and projected coordinates. 

542 ''' 

543 def __init__(self, lat0, lon0, datum=_WGS84, **name): 

544 '''New azimuthal L{EquidistantExact} projection. 

545 

546 @arg lat0: Latitude of center point (C{degrees90}). 

547 @arg lon0: Longitude of center point (C{degrees180}). 

548 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid}, 

549 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth 

550 radius (C{meter}). 

551 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}). 

552 

553 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}. 

554 ''' 

555 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, **name) 

556 

557 if _FOR_DOCS: 

558 forward = _EquidistantBase.forward 

559 reverse = _EquidistantBase.reverse 

560 

561 @Property_RO 

562 def geodesic(self): 

563 '''Get this projection's exact geodesic (L{GeodesicExact}). 

564 ''' 

565 return self.datum.ellipsoid.geodesicx 

566 

567 

568class EquidistantGeodSolve(_EquidistantBase): 

569 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant 

570 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>}, 

571 based on (exact) geodesic I{wrappers} L{GeodesicSolve} and L{GeodesicLineSolve} and intended 

572 I{for testing purposes only}. 

573 

574 @see: L{EquidistantExact} and module L{geodsolve}. 

575 ''' 

576 def __init__(self, lat0, lon0, datum=_WGS84, **name): 

577 '''New azimuthal L{EquidistantGeodSolve} projection. 

578 

579 @arg lat0: Latitude of center point (C{degrees90}). 

580 @arg lon0: Longitude of center point (C{degrees180}). 

581 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid}, 

582 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth 

583 radius (C{meter}). 

584 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}). 

585 

586 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}. 

587 ''' 

588 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, **name) 

589 

590 if _FOR_DOCS: 

591 forward = _EquidistantBase.forward 

592 reverse = _EquidistantBase.reverse 

593 

594 @Property_RO 

595 def geodesic(self): 

596 '''Get this projection's (exact) geodesic (L{GeodesicSolve}). 

597 ''' 

598 return self.datum.ellipsoid.geodsolve 

599 

600 

601class EquidistantKarney(_EquidistantBase): 

602 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant 

603 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>}, 

604 requiring package U{geographiclib<https://PyPI.org/project/geographiclib>} to be installed. 

605 

606 @see: L{EquidistantExact}. 

607 ''' 

608 def __init__(self, lat0, lon0, datum=_WGS84, **name): 

609 '''New azimuthal L{EquidistantKarney} projection. 

610 

611 @arg lat0: Latitude of center point (C{degrees90}). 

612 @arg lon0: Longitude of center point (C{degrees180}). 

613 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid}, 

614 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth 

615 radius (C{meter}). 

616 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}). 

617 

618 @raise ImportError: Package U{geographiclib<https://PyPI.org/project/geographiclib>} 

619 not installed or not found. 

620 

621 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}. 

622 ''' 

623 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, **name) 

624 

625 if _FOR_DOCS: 

626 forward = _EquidistantBase.forward 

627 reverse = _EquidistantBase.reverse 

628 

629 @Property_RO 

630 def geodesic(self): 

631 '''Get this projection's I{wrapped} U{geodesic.Geodesic 

632 <https://GeographicLib.SourceForge.io/Python/doc/code.html>}, provided 

633 I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>} 

634 package is installed. 

635 ''' 

636 return self.datum.ellipsoid.geodesic 

637 

638 

639_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve, 

640 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI 

641 

642 

643class Gnomonic(_AzimuthalBase): 

644 '''Azimuthal gnomonic projection for the sphere***, see U{Snyder, pp 164-168 

645 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram 

646 <https://MathWorld.Wolfram.com/GnomonicProjection.html>}. 

647 ''' 

648 if _FOR_DOCS: 

649 __init__ = _AzimuthalBase.__init__ 

650 

651 def forward(self, lat, lon, **name): 

652 '''Convert a geodetic location to azimuthal equidistant east- and northing. 

653 

654 @arg lat: Latitude of the location (C{degrees90}). 

655 @arg lon: Longitude of the location (C{degrees180}). 

656 @kwarg name: Optional C{B{name}=NN} for the location (C{str}). 

657 

658 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)} 

659 with easting C{x} and northing C{y} of point in C{meter} and C{lat} 

660 and C{lon} in C{degrees} and C{azimuth} clockwise from true North. 

661 The C{scale} of the projection is C{1} in I{radial} direction and 

662 is C{1 / reciprocal} in the direction perpendicular to this. 

663 

664 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}. 

665 ''' 

666 def _k_t(c): 

667 t = c > EPS 

668 k = (_1_0 / c) if t else _1_0 

669 return k, t 

670 

671 return self._forward(lat, lon, name, _k_t) 

672 

673 def reverse(self, x, y, **name_LatLon_and_kwds): 

674 '''Convert an azimuthal equidistant location to geodetic lat- and longitude. 

675 

676 @arg x: Easting of the location (C{meter}). 

677 @arg y: Northing of the location (C{meter}). 

678 @kwarg name_LatLon_and_kwds: Optional C{B{name}=NN} and class C{B{LatLon}=None} 

679 for the location and optionally, additional B{C{LatLon}} keyword 

680 arguments, ignored if C{B{LatLon} is None}. 

681 

682 @return: The geodetic (C{LatLon}) or if C{B{LatLon} is None} an 

683 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}. 

684 

685 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon} 

686 in the range C{[-180..180] degrees}. The C{scale} of the 

687 projection is C{1} in I{radial} direction, C{azimuth} clockwise 

688 from true North and C{1 / reciprocal} in the direction 

689 perpendicular to this. 

690 ''' 

691 def _c(c): 

692 return atan1(c) if c > EPS else None 

693 

694 return self._reverse(x, y, _c, False, **name_LatLon_and_kwds) 

695 

696 

697def gnomonic(lat0, lon0, datum=_WGS84, exact=False, geodsolve=False, **name): 

698 '''Return a L{GnomonicExact} or (if I{Karney}'s U{geographiclib 

699 <https://PyPI.org/project/geographiclib>} package is installed) 

700 a L{GnomonicKarney}, otherwise a L{Gnomonic} instance. 

701 

702 @arg lat0: Latitude of center point (C{degrees90}). 

703 @arg lon0: Longitude of center point (C{degrees180}). 

704 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid}, 

705 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth 

706 radius (C{meter}). 

707 @kwarg exact: Return a L{GnomonicExact} instance. 

708 @kwarg geodsolve: Return a L{GnomonicGeodSolve} instance. 

709 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}). 

710 

711 @return: A L{GnomonicExact}, L{GnomonicGeodSolve}, 

712 L{GnomonicKarney} or L{Gnomonic} instance. 

713 

714 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or 

715 (spherical) B{C{datum}}. 

716 

717 @raise GeodesicError: Issue with L{GeodesicExact}, L{GeodesicSolve} 

718 or I{Karney}'s wrapped C{Geodesic}. 

719 

720 @raise TypeError: Invalid B{C{datum}}. 

721 ''' 

722 G = GnomonicExact if exact else (GnomonicGeodSolve if geodsolve else Gnomonic) 

723 if G is Gnomonic: 

724 try: 

725 return GnomonicKarney(lat0, lon0, datum=datum, **name) # PYCHOK types 

726 except ImportError: 

727 pass 

728 return G(lat0, lon0, datum=datum, **name) # PYCHOK types 

729 

730 

731class _GnomonicBase(_AzimuthalGeodesic): 

732 '''(INTERNAL) Base for classes L{GnomonicExact}, L{GnomonicGeodSolve} 

733 and L{GnomonicKarney}. 

734 ''' 

735 def __init__(self, lat0, lon0, datum=_WGS84, **name): 

736 '''New azimuthal L{GnomonicExact} or L{GnomonicKarney} projection. 

737 ''' 

738 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, **name) 

739 

740 g = self.geodesic 

741 self._mask = g.ALL # | g.LONG_UNROLL 

742 

743 def forward(self, lat, lon, raiser=True, **name): # PYCHOK signature 

744 '''Convert an (ellipsoidal) geodetic location to azimuthal gnomonic east- 

745 and northing. 

746 

747 @arg lat: Latitude of the location (C{degrees90}). 

748 @arg lon: Longitude of the location (C{degrees180}). 

749 @kwarg raiser: If C{False}, do not throw an error if the location lies 

750 over the horizon (C{bool}). 

751 @kwarg name: Optional C{B{name}=NN} for the location (C{str}). 

752 

753 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)} 

754 with easting C{x} and northing C{y} in C{meter} and C{lat} and 

755 C{lon} in C{degrees} and C{azimuth} clockwise from true North. 

756 The C{scale} of the projection is C{1 / reciprocal**2} in I{radial} 

757 direction and C{1 / reciprocal} in the direction perpendicular to 

758 this. Both C{x} and C{y} will be C{NAN} if the (geodetic) location 

759 lies over the horizon and C{B{raiser}=False}. 

760 

761 @raise AzimuthalError: Invalid B{C{lat}}, B{C{lon}} or the location lies 

762 over the horizon and C{B{raiser}=True}. 

763 ''' 

764 r = self.geodesic.Inverse(self.lat0, self.lon0, 

765 Lat_(lat), Lon_(lon), outmask=self._mask) 

766 M = r.M21 

767 if M > EPS0: 

768 q = r.m12 / M # .M12 

769 e, n = sincos2d(r.azi1) 

770 e *= q 

771 n *= q 

772 elif raiser: # PYCHOK no cover 

773 raise AzimuthalError(lat=lat, lon=lon, txt=_over_horizon_) 

774 else: # PYCHOK no cover 

775 e = n = NAN 

776 

777 t = self._7Tuple(e, n, r, _name__(name), M=M) 

778 t._iteration = self._iteration = 0 

779 return t 

780 

781 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds): # PYCHOK signature 

782 '''Convert an azimuthal gnomonic location to (ellipsoidal) geodetic lat- and longitude. 

783 

784 @arg x: Easting of the location (C{meter}). 

785 @arg y: Northing of the location (C{meter}). 

786 @kwarg LatLon: Class to use (C{LatLon}) or C{None}. 

787 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} for the location and optionally, 

788 additional B{C{LatLon}} keyword arguments, ignored if C{B{LatLon} is None}. 

789 

790 @return: The geodetic (C{LatLon}) or if C{B{LatLon} is None} an 

791 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}. 

792 

793 @raise AzimuthalError: No convergence. 

794 

795 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon} in the range 

796 C{[-180..180] degrees}. The C{azimuth} is clockwise from true North. The 

797 scale is C{1 / reciprocal**2} in C{radial} direction and C{1 / reciprocal} 

798 in the direction perpendicular to this. 

799 ''' 

800 e, n, z, q = _enzh4(x, y) 

801 

802 d = a = self.equatoradius 

803 s = a * atan1(q, a) 

804 if q > a: # PYCHOK no cover 

805 def _d(r, q): 

806 return (r.M12 - q * r.m12) * r.m12 # negated 

807 

808 q = _1_0 / q 

809 else: # little == True 

810 def _d(r, q): # PYCHOK _d 

811 return (q * r.M12 - r.m12) * r.M12 # negated 

812 

813 a *= _EPS_K 

814 m = self._mask 

815 g = self.geodesic 

816 

817 _P = g.Line(self.lat0, self.lon0, z, caps=m | g.LINE_OFF).Position 

818 _S2 = Fsum(s).fsum2f_ 

819 _abs = fabs 

820 for i in range(1, _TRIPS): 

821 r = _P(s, outmask=m) 

822 if _abs(d) < a: 

823 break 

824 s, d = _S2(_d(r, q)) 

825 else: # PYCHOK no cover 

826 self._iteration = _TRIPS 

827 raise AzimuthalError(Fmt.no_convergence(d, a), 

828 txt=unstr(self.reverse, x, y)) 

829 

830 t = self._7Tuple(e, n, r, name_LatLon_kwds, M=r.M12) if LatLon is None else \ 

831 self._toLatLon(r.lat2, r.lon2, LatLon, name_LatLon_kwds) 

832 t._iteration = self._iteration = i 

833 return t 

834 

835 

836class GnomonicExact(_GnomonicBase): 

837 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic 

838 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>}, 

839 based on exact geodesic classes L{GeodesicExact} and L{GeodesicLineExact}. 

840 

841 @see: I{Karney}'s U{Detailed Description<https://GeographicLib.SourceForge.io/C++/doc/ 

842 classGeographicLib_1_1Gnomonic.html>}, especially the B{Warning}. 

843 ''' 

844 def __init__(self, lat0, lon0, datum=_WGS84, **name): 

845 '''New azimuthal L{GnomonicExact} projection. 

846 

847 @arg lat0: Latitude of center point (C{degrees90}). 

848 @arg lon0: Longitude of center point (C{degrees180}). 

849 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid}, 

850 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth 

851 radius (C{meter}). 

852 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}). 

853 

854 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}. 

855 ''' 

856 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, **name) 

857 

858 if _FOR_DOCS: 

859 forward = _GnomonicBase.forward 

860 reverse = _GnomonicBase.reverse 

861 

862 @Property_RO 

863 def geodesic(self): 

864 '''Get this projection's exact geodesic (L{GeodesicExact}). 

865 ''' 

866 return self.datum.ellipsoid.geodesicx 

867 

868 

869class GnomonicGeodSolve(_GnomonicBase): 

870 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic 

871 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>}, 

872 based on (exact) geodesic I{wrappers} L{GeodesicSolve} and L{GeodesicLineSolve} and 

873 intended I{for testing purposes only}. 

874 

875 @see: L{GnomonicExact} and module L{geodsolve}. 

876 ''' 

877 def __init__(self, lat0, lon0, datum=_WGS84, **name): 

878 '''New azimuthal L{GnomonicGeodSolve} projection. 

879 

880 @arg lat0: Latitude of center point (C{degrees90}). 

881 @arg lon0: Longitude of center point (C{degrees180}). 

882 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid}, 

883 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth 

884 radius (C{meter}). 

885 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}). 

886 

887 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}. 

888 ''' 

889 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, **name) 

890 

891 if _FOR_DOCS: 

892 forward = _GnomonicBase.forward 

893 reverse = _GnomonicBase.reverse 

894 

895 @Property_RO 

896 def geodesic(self): 

897 '''Get this projection's (exact) geodesic (L{GeodesicSolve}). 

898 ''' 

899 return self.datum.ellipsoid.geodsolve 

900 

901 

902class GnomonicKarney(_GnomonicBase): 

903 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic 

904 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>}, 

905 requiring package U{geographiclib<https://PyPI.org/project/geographiclib>} to be installed. 

906 

907 @see: L{GnomonicExact}. 

908 ''' 

909 def __init__(self, lat0, lon0, datum=_WGS84, **name): 

910 '''New azimuthal L{GnomonicKarney} projection.