Coverage for pygeodesy/azimuthal.py: 99%

316 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-05-06 16:50 -0400

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 _ # PYCHOK semicolon 

44 

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

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

47 _EPStol, _umod_360, _0_0, _0_1, \ 

48 _0_5, _1_0, _N_1_0, _2_0 

49from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB, \ 

50 _xinstanceof 

51from pygeodesy.datums import _spherical_datum, _WGS84 

52from pygeodesy.errors import _ValueError, _xdatum, _xkwds 

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

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

55# from pygeodesy.formy import antipode # _MODS 

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

57 _scale_, _SPACE_, _x_, _y_ 

58from pygeodesy.karney import _norm180 

59from pygeodesy.latlonBase import _MODS, LatLonBase as _LLB 

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

61from pygeodesy.named import _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 Bearing, Easting, Lat_, Lon_, Northing, \ 

67 Scalar, Scalar_ 

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

69 sincos2d, sincos2d_ 

70 

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

72 

73__all__ = _ALL_LAZY.azimuthal 

74__version__ = '24.04.07' 

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 the 

96 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=NN): 

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 name 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 datum not in (None, self._datum): 

117 self._datum = _spherical_datum(datum, name=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=NN): # PYCHOK no cover 

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

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

149 

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

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_2(s0 * sa + c0 * cb) 

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=name or self.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 if hasattr(latlon0, _datum_): 

199 _xdatum(self.datum, latlon0.datum, Error=AzimuthalError) 

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

201 

202 @Property_RO 

203 def lon0(self): 

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

205 ''' 

206 return self._latlon0.lon 

207 

208 @deprecated_Property_RO 

209 def majoradius(self): # PYCHOK no cover 

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

211 return self.equatoradius 

212 

213 @Property_RO 

214 def radius(self): 

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

216 ''' 

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

218 

219 def reset(self, lat0, lon0): 

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

221 

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

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

224 

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

226 ''' 

227 _update_all(self) # zap caches 

228 self._reset(lat0, lon0) 

229 

230 def _reset(self, lat0, lon0): 

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

232 ''' 

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

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

235 self._sc0 = sincos2d(self.lat0) 

236 

237 def reverse(self, x, y, name=NN, **LatLon_and_kwds): # PYCHOK no cover 

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

239 self._notOverloaded(x, y, name=name, **LatLon_and_kwds) 

240 

241 def _reverse(self, x, y, name, _c, lea, LatLon=None, **LatLon_kwds): 

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

243 ''' 

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

245 

246 c = _c(r / self.radius) 

247 if c is None: 

248 lat, lon = self.latlon0 

249 k, z = _1_0, _0_0 

250 else: 

251 s0, c0 = self._sc0 

252 sc, cc = sincos2(c) 

253 k = c / sc 

254 s = s0 * cc 

255 if r > EPS0: 

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

257 lat = degrees(asin1(s)) 

258 if lea or fabs(c0) > EPS: 

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

260 else: 

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

262 lon = _norm180(self.lon0 + d) 

263 

264 if LatLon is None: 

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

266 name=name or self.name) 

267 else: 

268 t = self._toLatLon(lat, lon, LatLon, LatLon_kwds, name) 

269 return t 

270 

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

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

273 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne. 

274 ''' 

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

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

277 return t, d 

278 

279 def _toLatLon(self, lat, lon, LatLon, LatLon_kwds, name): 

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

281 ''' 

282 kwds = _xkwds(LatLon_kwds, datum=self.datum) 

283 r = self._xnamed(LatLon(lat, lon, **kwds), name=name) # handle .classof 

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

285 _xinstanceof(B, LatLon=r) 

286 return r 

287 

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

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

290 

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

292 

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

294 (C{str}). 

295 ''' 

296 return _fstrLL0(self, prec, True) 

297 

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

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

300 

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

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

303 

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

305 ''' 

306 t = _fstrLL0(self, prec, False) 

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

308 

309 

310class AzimuthalError(_ValueError): 

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

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

313 L{Azimuthal7Tuple} issue. 

314 ''' 

315 pass 

316 

317 

318class Azimuthal7Tuple(_NamedTuple): 

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

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

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

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

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

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

325 ''' 

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

327 _Units_ = ( Easting, Northing, Lat_, Lon_, Bearing, Scalar, _Pass) 

328 

329 def antipodal(self, azimuth=None): 

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

331 

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

333 (C{compass degrees360}). 

334 ''' 

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

336 z = self.azimuth if azimuth is None else Bearing(azimuth=azimuth) # PYCHOK named 

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

338 

339 

340class Equidistant(_AzimuthalBase): 

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

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

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

344 

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

346 L{EquidistantGeodSolve} or L{EquidistantKarney} projection 

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

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

349 ''' 

350 if _FOR_DOCS: 

351 __init__ = _AzimuthalBase.__init__ 

352 

353 def forward(self, lat, lon, name=NN): 

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

355 

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

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

358 @kwarg name: Optional name for the location (C{str}). 

359 

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

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

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

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

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

365 

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

367 

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

369 the projection center C{(lat0, lon0)}. 

370 ''' 

371 def _k_t(c): 

372 k = _N_1_0 if c < 0 else _1_0 

373 t = fabs(c) < EPS1 

374 if t: 

375 a = acos(c) 

376 s = sin(a) 

377 if s: 

378 k = a / s 

379 return k, t 

380 

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

382 

383 def reverse(self, x, y, name=NN, **LatLon_and_kwds): 

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

385 

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

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

388 @kwarg name: Optional name for the location (C{str}). 

389 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and 

390 additional B{C{LatLon}} keyword arguments, 

391 ignored if C{B{LatLon} is None} or not given. 

392 

393 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{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, name, _c, False, **LatLon_and_kwds) 

406 

407 

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

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 name 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=name) # PYCHOK types 

437 except ImportError: 

438 pass 

439 return E(lat0, lon0, datum=datum, name=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, M=None, name=NN): 

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

457 ''' 

458 s = M if M is not None else ( # reciprocal, azimuthal scale 

459 (r.m12 / r.s12) if r.a12 > _EPS_K else _1_0) 

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

461 return Azimuthal7Tuple(e, n, r.lat2, r.lon2, z, s, self.datum, 

462 name=name or self.name) 

463 

464 

465class _EquidistantBase(_AzimuthalGeodesic): 

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

467 and L{EquidistantKarney}. 

468 ''' 

469 def __init__(self, lat0, lon0, datum=_WGS84, name=NN): 

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

471 L{EquidistantKarney} projection. 

472 ''' 

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

474 

475 g = self.geodesic 

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

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

478 

479 def forward(self, lat, lon, name=NN): 

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

481 

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

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

484 @kwarg name: Optional name for the location (C{str}). 

485 

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

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

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

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

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

491 

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

493 

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

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

496 ''' 

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

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

499 x, y = sincos2d(r.azi1) 

500 return self._7Tuple(x * r.s12, y * r.s12, r, name=name) 

501 

502 def reverse(self, x, y, name=NN, LatLon=None, **LatLon_kwds): # PYCHOK signature 

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

504 

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

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

507 @kwarg name: Optional name for the location (C{str}). 

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

509 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword 

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

511 

512 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{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=name) if LatLon is None else \ 

525 self._toLatLon(r.lat2, r.lon2, LatLon, LatLon_kwds, name) 

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=NN): 

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 name 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=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=NN): 

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 name 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=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=NN): 

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 name 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=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=NN): 

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 name 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=NN, **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: Optional name for the location (C{str}). 

679 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and 

680 additional B{C{LatLon}} keyword arguments, 

681 ignored if C{B{LatLon} is None} or not given. 

682 

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

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

685 

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

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

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

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

690 perpendicular to this. 

691 ''' 

692 def _c(c): 

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

694 

695 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds) 

696 

697 

698def gnomonic(lat0, lon0, datum=_WGS84, exact=False, geodsolve=False, name=NN): 

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

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

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

702 

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

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

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

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

707 radius (C{meter}). 

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

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

710 @kwarg name: Optional name for the projection (C{str}). 

711 

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

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

714 

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

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

717 

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

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

720 

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

722 ''' 

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

724 if G is Gnomonic: 

725 try: 

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

727 except ImportError: 

728 pass 

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

730 

731 

732class _GnomonicBase(_AzimuthalGeodesic): 

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

734 and L{GnomonicKarney}. 

735 ''' 

736 def __init__(self, lat0, lon0, datum=_WGS84, name=NN): 

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

738 ''' 

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

740 

741 g = self.geodesic 

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

743 

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

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

746 and northing. 

747 

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

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

750 @kwarg name: Optional name for the location (C{str}). 

751 @kwarg raiser: Do or don't throw an error (C{bool}) if 

752 the location lies over the horizon. 

753 

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

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

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

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

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

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

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

761 

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

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

764 ''' 

765 self._iteration = 0 

766 

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

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

769 M = r.M21 

770 if M > EPS0: 

771 q = r.m12 / M # .M12 

772 e, n = sincos2d(r.azi1) 

773 e *= q 

774 n *= q 

775 elif raiser: # PYCHOK no cover 

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

777 else: # PYCHOK no cover 

778 e = n = NAN 

779 

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

781 t._iteraton = 0 

782 return t 

783 

784 def reverse(self, x, y, name=NN, LatLon=None, **LatLon_kwds): # PYCHOK signature 

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

786 

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

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

789 @kwarg name: Optional name for the location (C{str}). 

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

791 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword 

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

793 

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

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

796 

797 @raise AzimuthalError: No convergence. 

798 

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

800 in the range C{[-180..180] degrees}. The C{azimuth} is clockwise 

801 from true North. The scale is C{1 / reciprocal**2} in C{radial} 

802 direction and C{1 / reciprocal} in the direction perpendicular 

803 to this. 

804 ''' 

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

806 

807 d = a = self.equatoradius 

808 s = a * atan1(q, a) 

809 if q > a: # PYCHOK no cover 

810 def _d(r, q): 

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

812 

813 q = _1_0 / q 

814 else: # little == True 

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

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

817 

818 a *= _EPS_K 

819 m = self._mask 

820 g = self.geodesic 

821 

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

823 _S2 = Fsum(s).fsum2f_ 

824 _abs = fabs 

825 for i in range(1, _TRIPS): 

826 r = _P(s, outmask=m) 

827 if _abs(d) < a: 

828 break 

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

830 else: # PYCHOK no cover 

831 self._iteration = _TRIPS 

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

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

834 

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

836 self._toLatLon(r.lat2, r.lon2, LatLon, LatLon_kwds, name) 

837 t._iteration = self._iteration = i 

838 return t 

839 

840 

841class GnomonicExact(_GnomonicBase): 

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

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

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

845 

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

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

848 ''' 

849 def __init__(self, lat0, lon0, datum=_WGS84, name=NN): 

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

851 

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

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

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

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

856 radius (C{meter}). 

857 @kwarg name: Optional name for the projection (C{str}). 

858 

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

860 ''' 

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

862 

863 if _FOR_DOCS: 

864 forward = _GnomonicBase.forward 

865 reverse = _GnomonicBase.reverse 

866 

867 @Property_RO 

868 def geodesic(self): 

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

870 ''' 

871 return self.datum.ellipsoid.geodesicx 

872 

873 

874class GnomonicGeodSolve(_GnomonicBase): 

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

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

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

878 intended I{for testing purposes only}. 

879 

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

881 ''' 

882 def __init__(self, lat0, lon0, datum=_WGS84, name=NN): 

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

884 

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

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

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

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

889 radius (C{meter}). 

890 @kwarg name: Optional name for the projection (C{str}). 

891 

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

893 ''' 

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

895 

896 if _FOR_DOCS: 

897 forward = _GnomonicBase.forward 

898 reverse = _GnomonicBase.reverse 

899 

900 @Property_RO 

901 def geodesic(self): 

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

903 ''' 

904 return self.datum.ellipsoid.geodsolve 

905 

906 

907class GnomonicKarney(_GnomonicBase): 

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

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

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

911 

912 @see: L{GnomonicExact}. 

913 ''' 

914 def __init__(self, lat0, lon0, datum=_WGS84, name=NN): 

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

916 

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

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

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

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

921 radius (C{meter}). 

922 @kwarg name: Optional name for the projection (C{str}). 

923 

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

925 not installed or not found. 

926 

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

928 ''' 

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

930 

931 if _FOR_DOCS: 

932 forward = _GnomonicBase.forward 

933 reverse = _GnomonicBase.reverse 

934 

935 @Property_RO 

936 def geodesic(self): 

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

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

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

940 package is installed. 

941 ''' 

942 return self.datum.ellipsoid.geodesic 

943 

944 

945class LambertEqualArea(_AzimuthalBase): 

946 '''Lambert-equal-area projection for the sphere*** (aka U{Lambert zenithal equal-area 

947 projection<https://WikiPedia.org/wiki/Lambert_azimuthal_equal-area_projection>}, see 

948 U{Snyder, pp 185-187<https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram 

949 <https://MathWorld.Wolfram.com/LambertAzimuthalEqual-AreaProjection.html>}. 

950 ''' 

951 if _FOR_DOCS: 

952 __init__ = _AzimuthalBase.__init__ 

953 

954 def forward(self, lat, lon, name=NN): 

955 '''Convert a geodetic location to azimuthal Lambert-equal-area east- and northing. 

956 

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

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

959 @kwarg name: Optional name for the location (C{str}). 

960 

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

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

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

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

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

966 

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

968 ''' 

969 def _k_t(c): 

970 c += _1_0 

971 t = c > EPS0 

972 k = sqrt(_2_0 / c) if t else _1_0 

973 return k, t 

974 

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

976 

977 def reverse(self, x, y, name=NN, **LatLon_and_kwds): 

978 '''Convert an azimuthal Lambert-equal-area location to geodetic lat- and longitude. 

979 

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

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

982 @kwarg name: Optional name for the location (C{str}). 

983 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and 

984 additional B{C{LatLon}} keyword arguments, 

985 ignored if C{B{LatLon} is None} or not given. 

986 

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

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

989 

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

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

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

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

994 perpendicular to this. 

995 ''' 

996 def _c(c): 

997 c *= _0_5 

998 return (asin1(c) * _2_0) if c > EPS else None 

999 

1000 return self._reverse(x, y, name, _c, True, **LatLon_and_kwds) 

1001 

1002 

1003class Orthographic(_AzimuthalBase): 

1004 '''Orthographic projection for the sphere***, see U{Snyder, pp 148-153 

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

1006 <https://MathWorld.Wolfram.com/OrthographicProjection.html>}. 

1007 ''' 

1008 if _FOR_DOCS: 

1009 __init__ = _AzimuthalBase.__init__ 

1010 

1011 def forward(self, lat, lon, name=NN): 

1012 '''Convert a geodetic location to azimuthal orthographic east- and northing. 

1013 

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

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

1016 @kwarg name: Optional name for the location (C{str}). 

1017 

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

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

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

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

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

1023 

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

1025 ''' 

1026 def _k_t(c): 

1027 return _1_0, (c >= 0) 

1028 

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

1030 

1031 def reverse(self, x, y, name=NN, **LatLon_and_kwds): 

1032 '''Convert an azimuthal orthographic location to geodetic lat- and longitude. 

1033 

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

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

1036 @kwarg name: Optional name for the location (C{str}). 

1037 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and 

1038 additional B{C{LatLon}} keyword arguments, 

1039 ignored if C{B{LatLon} is None} or not given. 

1040 

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

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

1043 

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

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

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

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

1048 perpendicular to this. 

1049 ''' 

1050 def _c(c): 

1051 return asin1(c) if c > EPS else None 

1052 

1053 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds) 

1054 

1055 

1056class Stereographic(_AzimuthalBase): 

1057 '''Stereographic projection for the sphere***, see U{Snyder, pp 157-160 

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

1059 <https://MathWorld.Wolfram.com/StereographicProjection.html>}. 

1060 ''' 

1061 _k0 = _1_0 # central scale factor (C{scalar}) 

1062 _k02 = _2_0 # double ._k0 

1063 

1064 if _FOR_DOCS: 

1065 __init__ = _AzimuthalBase.__init__ 

1066 

1067 def forward(self, lat, lon, name=NN): 

1068 '''Convert a geodetic location to azimuthal stereographic east- and northing. 

1069 

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

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

1072 @kwarg name: Optional name for the location (C{str}). 

1073 

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

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

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

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

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

1079 

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

1081 ''' 

1082 def _k_t(c): 

1083 c += _1_0 

1084 t = isnon0(c) 

1085 k = (self._k02 / c) if t else _1_0 

1086 return k, t 

1087 

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

1089 

1090 @property_doc_(''' optional, central scale factor (C{scalar}).''') 

1091 def k0(self): 

1092 '''Get the central scale factor (C{scalar}). 

1093 ''' 

1094 return self._k0 

1095 

1096 @k0.setter # PYCHOK setter! 

1097 def k0(self, factor): 

1098 '''Set the central scale factor (C{scalar}). 

1099 ''' 

1100 n = Stereographic.k0.fget.__name__ 

1101 self._k0 = Scalar_(factor, name=n, low=EPS, high=2) # XXX high=1, 2, other? 

1102 self._k02 = self._k0 * _2_0 

1103 

1104 def reverse(self, x, y, name=NN, **LatLon_and_kwds): 

1105 '''Convert an azimuthal stereographic location to geodetic lat- and longitude. 

1106 

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

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

1109 @kwarg name: Optional name for the location (C{str}). 

1110 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and 

1111 additional B{C{LatLon}} keyword arguments, 

1112 ignored if C{B{LatLon} is None} or not given. 

1113 

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

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

1116 

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

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

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

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

1121 perpendicular to this. 

1122 ''' 

1123 def _c(c): 

1124 return (atan2(c, self._k02) * _2_0) if c > EPS else None 

1125 

1126 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds) 

1127 

1128 

1129__all__ += _ALL_DOCS(_AzimuthalBase, _AzimuthalGeodesic, _EquidistantBase, _GnomonicBase) 

1130 

1131# **) MIT License 

1132# 

1133# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved. 

1134# 

1135# Permission is hereby granted, free of charge, to any person obtaining a 

1136# copy of this software and associated documentation files (the "Software"), 

1137# to deal in the Software without restriction, including without limitation 

1138# the rights to use, copy, modify, merge, publish, distribute, sublicense, 

1139# and/or sell copies of the Software, and to permit persons to whom the 

1140# Software is furnished to do so, subject to the following conditions: 

1141# 

1142# The above copyright notice and this permission notice shall be included 

1143# in all copies or substantial portions of the Software. 

1144# 

1145# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 

1146# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 

1147# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 

1148# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 

1149# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 

1150# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 

1151# OTHER DEALINGS IN THE SOFTWARE.