Coverage for pygeodesy/azimuthal.py: 98%

312 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-03-31 10:52 -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 .datums 

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 _spherical_datum, _WGS84, _xinstanceof 

51# from pygeodesy.datums import _spherical_datum, _WGS84, _xinstanceof 

52from pygeodesy.errors import _ValueError, _xdatum, _xkwds 

53from pygeodesy.fmath import euclid, Fsum, hypot 

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

55# from pygeodesy.formy import antipode # from latlonBase 

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

57 _scale_, _SPACE_, _x_, _y_ 

58from pygeodesy.karney import _norm180 

59from pygeodesy.latlonBase import antipode, LatLonBase as _LLB 

60from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS 

61from pygeodesy.named import _NamedBase, _NamedTuple, notOverloaded, _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_, \ 

67 Northing, Scalar, Scalar_ 

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

69 sincos2d, sincos2d_ 

70 

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

72 

73__all__ = _ALL_LAZY.azimuthal 

74__version__ = '23.03.21' 

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 html/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 @Property_RO 

137 def flattening(self): 

138 '''Get the geodesic's flattening (C{float}). 

139 ''' 

140 return self.datum.ellipsoid.f 

141 

142 def forward(self, lat, lon, name=NN): # PYCHOK no cover 

143 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}. 

144 ''' 

145 notOverloaded(self, lat, lon, name=name) 

146 

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

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

149 ''' 

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

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

152 s0, c0 = self._sc0 

153 

154 cb *= ca 

155 k, t = _k_t_2(s0 * sa + c0 * cb) 

156 if t: 

157 r = k * self.radius 

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

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

160 else: # 0 or 180 

161 e = n = z = _0_0 

162 

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

164 name=name or self.name) 

165 return t 

166 

167 def _forwards(self, *lls): 

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

169 ''' 

170 _fwd = self.forward 

171 for ll in lls: 

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

173 

174 @Property_RO 

175 def lat0(self): 

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

177 ''' 

178 return self._latlon0.lat 

179 

180 @property 

181 def latlon0(self): 

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

183 ''' 

184 return self._latlon0 

185 

186 @latlon0.setter # PYCHOK setter! 

187 def latlon0(self, latlon0): 

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

189 

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

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

192 ''' 

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

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

195 if hasattr(latlon0, _datum_): 

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

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

198 

199 @Property_RO 

200 def lon0(self): 

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

202 ''' 

203 return self._latlon0.lon 

204 

205 @deprecated_Property_RO 

206 def majoradius(self): # PYCHOK no cover 

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

208 return self.equatoradius 

209 

210 @Property_RO 

211 def radius(self): 

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

213 ''' 

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

215 

216 def reset(self, lat0, lon0): 

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

218 

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

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

221 

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

223 ''' 

224 _update_all(self) # zap caches 

225 self._reset(lat0, lon0) 

226 

227 def _reset(self, lat0, lon0): 

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

229 ''' 

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

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

232 self._sc0 = sincos2d(self.lat0) 

233 

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

235 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}. 

236 ''' 

237 notOverloaded(self, x, y, name=name, **LatLon_and_kwds) 

238 

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

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

241 ''' 

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

243 

244 c = _c(r / self.radius) 

245 if c is None: 

246 lat, lon = self.latlon0 

247 k, z = _1_0, _0_0 

248 else: 

249 s0, c0 = self._sc0 

250 sc, cc = sincos2(c) 

251 k = c / sc 

252 s = s0 * cc 

253 if r > EPS0: 

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

255 lat = degrees(asin1(s)) 

256 if lea or fabs(c0) > EPS: 

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

258 else: 

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

260 lon = _norm180(self.lon0 + d) 

261 

262 if LatLon is None: 

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

264 name=name or self.name) 

265 else: 

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

267 return t 

268 

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

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

271 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne. 

272 ''' 

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

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

275 return t, d 

276 

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

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

279 ''' 

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

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

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

283 _xinstanceof(B, LatLon=r) 

284 return r 

285 

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

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

288 

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

290 

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

292 (C{str}). 

293 ''' 

294 return _fstrLL0(self, prec, True) 

295 

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

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

298 

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

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

301 

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

303 ''' 

304 t = _fstrLL0(self, prec, False) 

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

306 

307 

308class AzimuthalError(_ValueError): 

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

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

311 L{Azimuthal7Tuple} issue. 

312 ''' 

313 pass 

314 

315 

316class Azimuthal7Tuple(_NamedTuple): 

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

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

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

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

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

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

323 ''' 

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

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

326 

327 def antipodal(self, azimuth=None): 

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

329 

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

331 (C{compass degrees360}). 

332 ''' 

333 a = antipode(self.lat, self.lon) # PYCHOK named 

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

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

336 

337 

338class Equidistant(_AzimuthalBase): 

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

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

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

342 

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

344 L{EquidistantGeodSolve} or L{EquidistantKarney} projection 

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

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

347 ''' 

348 if _FOR_DOCS: 

349 __init__ = _AzimuthalBase.__init__ 

350 

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

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

353 

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

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

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

357 

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

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

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

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

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

363 

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

365 

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

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

368 ''' 

369 def _k_t(c): 

370 k = _N_1_0 if c < 0 else _1_0 

371 t = fabs(c) < EPS1 

372 if t: 

373 a = acos(c) 

374 s = sin(a) 

375 if s: 

376 k = a / s 

377 return k, t 

378 

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

380 

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

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

383 

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

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

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

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

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

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

390 

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

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

393 

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

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

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

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

398 perpendicular to this. 

399 ''' 

400 def _c(c): 

401 return c if c > EPS else None 

402 

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

404 

405 

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

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

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

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

410 

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

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

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

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

415 radius (C{meter}). 

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

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

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

419 

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

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

422 

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

424 

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

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

427 

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

429 ''' 

430 

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

432 if E is Equidistant: 

433 try: 

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

435 except ImportError: 

436 pass 

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

438 

439 

440class _AzimuthalGeodesic(_AzimuthalBase): 

441 '''(INTERNAL) Base class for azimuthal projections using U{Karney Geodesic 

442 <https://GeographicLib.SourceForge.io/C++/doc/python/code.html>} or 

443 exact geodesic classes L{GeodesicExact} and L{GeodesicLineExact}. 

444 ''' 

445 _mask = 0 

446 

447 @Property_RO 

448 def geodesic(self): # PYCHOK no cover 

449 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}. 

450 ''' 

451 notOverloaded(self) 

452 

453 def _7Tuple(self, e, n, r, M=None, name=NN): 

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

455 ''' 

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

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

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

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

460 name=name or self.name) 

461 

462 

463class _EquidistantBase(_AzimuthalGeodesic): 

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

465 and L{EquidistantKarney}. 

466 ''' 

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

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

469 L{EquidistantKarney} projection. 

470 ''' 

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

472 

473 g = self.geodesic 

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

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

476 

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

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

479 

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

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

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

483 

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

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

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

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

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

489 

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

491 

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

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

494 ''' 

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

496 Lat_(lat), Lon_(lon), self._mask) 

497 x, y = sincos2d(r.azi1) 

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

499 

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

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

502 

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

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

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

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

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

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

509 

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

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

512 

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

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

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

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

517 to this. 

518 ''' 

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

520 

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

522 return self._7Tuple(e, n, r, name=name) if LatLon is None else \ 

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

524 

525 

526class EquidistantExact(_EquidistantBase): 

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

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

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

530 

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

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

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

534 clockwise from true North. 

535 

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

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

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

539 geographic and projected coordinates. 

540 ''' 

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

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

543 

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

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

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

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

548 radius (C{meter}). 

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

550 

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

552 ''' 

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

554 

555 if _FOR_DOCS: 

556 forward = _EquidistantBase.forward 

557 reverse = _EquidistantBase.reverse 

558 

559 @Property_RO 

560 def geodesic(self): 

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

562 ''' 

563 return self.datum.ellipsoid.geodesicx 

564 

565 

566class EquidistantGeodSolve(_EquidistantBase): 

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

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

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

570 I{for testing purposes only}. 

571 

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

573 ''' 

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

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

576 

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

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

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

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

581 radius (C{meter}). 

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

583 

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

585 ''' 

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

587 

588 if _FOR_DOCS: 

589 forward = _EquidistantBase.forward 

590 reverse = _EquidistantBase.reverse 

591 

592 @Property_RO 

593 def geodesic(self): 

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

595 ''' 

596 return self.datum.ellipsoid.geodsolve 

597 

598 

599class EquidistantKarney(_EquidistantBase): 

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

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

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

603 

604 @see: L{EquidistantExact}. 

605 ''' 

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

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

608 

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

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

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

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

613 radius (C{meter}). 

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

615 

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

617 not installed or not found. 

618 

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

620 ''' 

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

622 

623 if _FOR_DOCS: 

624 forward = _EquidistantBase.forward 

625 reverse = _EquidistantBase.reverse 

626 

627 @Property_RO 

628 def geodesic(self): 

629 '''Get this projection's I{wrapped} U{Karney Geodesic 

630 <https://GeographicLib.SourceForge.io/C++/doc/python/code.html>}, 

631 provided package U{geographiclib 

632 <https://PyPI.org/project/geographiclib>} is installed. 

633 ''' 

634 return self.datum.ellipsoid.geodesic 

635 

636 

637_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve, 

638 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI 

639 

640 

641class Gnomonic(_AzimuthalBase): 

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

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

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

645 ''' 

646 if _FOR_DOCS: 

647 __init__ = _AzimuthalBase.__init__ 

648 

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

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

651 

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

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

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

655 

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

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

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

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

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

661 

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

663 ''' 

664 def _k_t(c): 

665 t = c > EPS 

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

667 return k, t 

668 

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

670 

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

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

673 

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

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

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

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

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

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

680 

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

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

683 

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

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

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

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

688 perpendicular to this. 

689 ''' 

690 def _c(c): 

691 return atan(c) if c > EPS else None 

692 

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

694 

695 

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

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

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

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

700 

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

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

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

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

705 radius (C{meter}). 

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

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

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

709 

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

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

712 

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

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

715 

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

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

718 

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

720 ''' 

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

722 if G is Gnomonic: 

723 try: 

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

725 except ImportError: 

726 pass 

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

728 

729 

730class _GnomonicBase(_AzimuthalGeodesic): 

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

732 and L{GnomonicKarney}. 

733 ''' 

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

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

736 ''' 

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

738 

739 g = self.geodesic 

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

741 

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

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

744 and northing. 

745 

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

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

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

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

750 the location lies over the horizon. 

751 

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

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

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

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

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

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

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

759 

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

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

762 ''' 

763 self._iteration = 0 

764 

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

766 Lat_(lat), Lon_(lon), self._mask) 

767 M = r.M21 

768 if M > EPS0: 

769 q = r.m12 / M # .M12 

770 e, n = sincos2d(r.azi1) 

771 e *= q 

772 n *= q 

773 elif raiser: # PYCHOK no cover 

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

775 else: # PYCHOK no cover 

776 e = n = NAN 

777 

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

779 t._iteraton = 0 

780 return t 

781 

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

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

784 

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

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

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

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

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

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

791 

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

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

794 

795 @raise AzimuthalError: No convergence. 

796 

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

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

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

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

801 to this. 

802 ''' 

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

804 

805 d = a = self.equatoradius 

806 s = a * atan(q / a) 

807 if q > a: # PYCHOK no cover 

808 def _d(r, q): 

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

810 

811 q = _1_0 / q 

812 else: # little == True 

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

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

815 

816 a *= _EPS_K 

817 m = self._mask 

818 g = self.geodesic 

819 

820 P = g.Line(self.lat0, self.lon0, z, m | g.LINE_OFF).Position 

821 _S2 = Fsum(s).fsum2_ 

822 for i in range(1, _TRIPS): 

823 r = P(s, m) 

824 if fabs(d) < a: 

825 break 

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

827 else: # PYCHOK no cover 

828 self._iteration = _TRIPS 

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

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

831 

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

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

834 t._iteration = self._iteration = i 

835 return t 

836 

837 

838class GnomonicExact(_GnomonicBase): 

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

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

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

842 

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

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

845 ''' 

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

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

848 

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

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

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

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

853 radius (C{meter}). 

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

855 

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

857 ''' 

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

859 

860 if _FOR_DOCS: 

861 forward = _GnomonicBase.forward 

862 reverse = _GnomonicBase.reverse 

863 

864 @Property_RO 

865 def geodesic(self): 

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

867 ''' 

868 return self.datum.ellipsoid.geodesicx 

869 

870 

871class GnomonicGeodSolve(_GnomonicBase): 

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

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

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

875 intended I{for testing purposes only}. 

876 

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

878 ''' 

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

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

881 

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

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

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

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

886 radius (C{meter}). 

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

888 

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

890 ''' 

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

892 

893 if _FOR_DOCS: 

894 forward = _GnomonicBase.forward 

895 reverse = _GnomonicBase.reverse 

896 

897 @Property_RO 

898 def geodesic(self): 

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

900 ''' 

901 return self.datum.ellipsoid.geodsolve 

902 

903 

904class GnomonicKarney(_GnomonicBase): 

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

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

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

908 

909 @see: L{GnomonicExact}. 

910 ''' 

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

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

913 

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

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

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

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

918 radius (C{meter}). 

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

920 

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

922 not installed or not found. 

923 

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

925 ''' 

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

927 

928 if _FOR_DOCS: 

929 forward = _GnomonicBase.forward 

930 reverse = _GnomonicBase.reverse 

931 

932 @Property_RO 

933 def geodesic(self): 

934 '''Get this projection's I{wrapped} U{Karney Geodesic 

935 <https://GeographicLib.SourceForge.io/C++/doc/python/code.html>}, provided package 

936 U{geographiclib<https://PyPI.org/project/geographiclib>} is installed. 

937 ''' 

938 return self.datum.ellipsoid.geodesic 

939 

940 

941class LambertEqualArea(_AzimuthalBase): 

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

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

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

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

946 ''' 

947 if _FOR_DOCS: 

948 __init__ = _AzimuthalBase.__init__ 

949 

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

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

952 

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

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

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

956 

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

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

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

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

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

962 

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

964 ''' 

965 def _k_t(c): 

966 c += _1_0 

967 t = c > EPS0 

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

969 return k, t 

970 

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

972 

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

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

975 

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

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

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

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

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

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

982 

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

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

985 

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

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

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

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

990 perpendicular to this. 

991 ''' 

992 def _c(c): 

993 c *= _0_5 

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

995 

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

997 

998 

999class Orthographic(_AzimuthalBase): 

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

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

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

1003 ''' 

1004 if _FOR_DOCS: 

1005 __init__ = _AzimuthalBase.__init__ 

1006 

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

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

1009 

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

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

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

1013 

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

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

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

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

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

1019 

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

1021 ''' 

1022 def _k_t(c): 

1023 return _1_0, (c >= 0) 

1024 

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

1026 

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

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

1029 

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

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

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

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

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

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

1036 

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

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

1039 

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

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

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

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

1044 perpendicular to this. 

1045 ''' 

1046 def _c(c): 

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

1048 

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

1050 

1051 

1052class Stereographic(_AzimuthalBase): 

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

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

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

1056 ''' 

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

1058 _k02 = _2_0 # double ._k0 

1059 

1060 if _FOR_DOCS: 

1061 __init__ = _AzimuthalBase.__init__ 

1062 

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

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

1065 

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

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

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

1069 

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

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

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

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

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

1075 

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

1077 ''' 

1078 def _k_t(c): 

1079 c += _1_0 

1080 t = isnon0(c) 

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

1082 return k, t 

1083 

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

1085 

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

1087 def k0(self): 

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

1089 ''' 

1090 return self._k0 

1091 

1092 @k0.setter # PYCHOK setter! 

1093 def k0(self, factor): 

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

1095 ''' 

1096 n = Stereographic.k0.fget.__name__ 

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

1098 self._k02 = self._k0 * _2_0 

1099 

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

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

1102 

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

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

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

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

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

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

1109 

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

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

1112 

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

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

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

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

1117 perpendicular to this. 

1118 ''' 

1119 def _c(c): 

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

1121 

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

1123 

1124 

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

1126 

1127# **) MIT License 

1128# 

1129# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

1130# 

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

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

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

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

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

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

1137# 

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

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

1140# 

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

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

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

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

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

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

1147# OTHER DEALINGS IN THE SOFTWARE.