Coverage for pygeodesy/azimuthal.py: 99%

320 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-06-01 11:43 -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, _umod_360, \ 

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

48from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB, \ 

49 _xinstanceof 

50from pygeodesy.datums import _spherical_datum, _WGS84 

51from pygeodesy.errors import _ValueError, _xdatum, _xkwds 

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

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

54# from pygeodesy.formy import antipode # _MODS 

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

56 _SPACE_, _x_, _y_ 

57from pygeodesy.karney import _norm180 

58from pygeodesy.latlonBase import _MODS, LatLonBase as _LLB 

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

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

61from pygeodesy.namedTuples import LatLon2Tuple, LatLon4Tuple 

62from pygeodesy.props import deprecated_Property_RO, Property_RO, \ 

63 property_doc_, _update_all 

64from pygeodesy.streprs import Fmt, _fstrLL0, unstr 

65from pygeodesy.units import Bearing, Easting, Lat_, Lon_, Northing, \ 

66 Scalar, Scalar_ 

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

68 sincos2d, sincos2d_ 

69 

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

71 

72__all__ = _ALL_LAZY.azimuthal 

73__version__ = '24.05.24' 

74 

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

76_over_horizon_ = 'over horizon' 

77_TRIPS = 21 # numit, 4 sufficient 

78 

79 

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

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

82 ''' 

83 e = Easting( x=x) 

84 n = Northing(y=y) 

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

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

87 

88 

89class _AzimuthalBase(_NamedBase): 

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

91 

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

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

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

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

96 ''' 

97 _datum = _WGS84 # L{Datum} 

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

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

100 

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

102 '''New azimuthal projection. 

103 

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

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

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

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

108 radius (C{meter}). 

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

110 

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

112 

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

114 ''' 

115 if datum not in (None, self._datum): 

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

117 if name: 

118 self.name = name 

119 

120 if lat0 or lon0: # often both 0 

121 self._reset(lat0, lon0) 

122 

123 @Property_RO 

124 def datum(self): 

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

126 ''' 

127 return self._datum 

128 

129 @Property_RO 

130 def equatoradius(self): 

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

132 ''' 

133 return self.datum.ellipsoid.a 

134 

135 a = equatoradius 

136 

137 @Property_RO 

138 def flattening(self): 

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

140 ''' 

141 return self.datum.ellipsoid.f 

142 

143 f = flattening 

144 

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

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

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

148 

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

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

151 ''' 

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

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

154 s0, c0 = self._sc0 

155 

156 cb *= ca 

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

158 if t: 

159 r = k * self.radius 

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

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

162 else: # 0 or 180 

163 e = n = z = _0_0 

164 

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

166 name=self._name__(name)) 

167 return t 

168 

169 def _forwards(self, *lls): 

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

171 ''' 

172 _fwd = self.forward 

173 for ll in lls: 

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

175 

176 @Property_RO 

177 def lat0(self): 

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

179 ''' 

180 return self._latlon0.lat 

181 

182 @property 

183 def latlon0(self): 

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

185 ''' 

186 return self._latlon0 

187 

188 @latlon0.setter # PYCHOK setter! 

189 def latlon0(self, latlon0): 

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

191 

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

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

194 ''' 

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

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

197 if hasattr(latlon0, _datum_): 

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

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

200 

201 @Property_RO 

202 def lon0(self): 

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

204 ''' 

205 return self._latlon0.lon 

206 

207 @deprecated_Property_RO 

208 def majoradius(self): # PYCHOK no cover 

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

210 return self.equatoradius 

211 

212 @Property_RO 

213 def radius(self): 

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

215 ''' 

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

217 

218 def reset(self, lat0, lon0): 

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

220 

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

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

223 

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

225 ''' 

226 _update_all(self) # zap caches 

227 self._reset(lat0, lon0) 

228 

229 def _reset(self, lat0, lon0): 

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

231 ''' 

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

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

234 self._sc0 = sincos2d(self.lat0) 

235 

236 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds): # PYCHOK no cover 

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

238 self._notOverloaded(x, y, LatLon=LatLon, **name_LatLon_kwds) 

239 

240 def _reverse(self, x, y, _c, lea, LatLon, **name_LatLon_kwds): 

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

242 ''' 

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

244 

245 c = _c(r / self.radius) 

246 if c is None: 

247 lat, lon = self.latlon0 

248 k, z = _1_0, _0_0 

249 else: 

250 s0, c0 = self._sc0 

251 sc, cc = sincos2(c) 

252 k = c / sc 

253 s = s0 * cc 

254 if r > EPS0: 

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

256 lat = degrees(asin1(s)) 

257 if lea or fabs(c0) > EPS: 

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

259 else: 

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

261 lon = _norm180(self.lon0 + d) 

262 

263 if LatLon is None: 

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

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

266 else: 

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

268 return t 

269 

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

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

272 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne. 

273 ''' 

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

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

276 return t, d 

277 

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

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

280 ''' 

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

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

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

284 _xinstanceof(B, LatLon=r) 

285 return r 

286 

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

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

289 

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

291 

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

293 (C{str}). 

294 ''' 

295 return _fstrLL0(self, prec, True) 

296 

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

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

299 

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

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

302 

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

304 ''' 

305 t = _fstrLL0(self, prec, False) 

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

307 

308 

309class AzimuthalError(_ValueError): 

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

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

312 L{Azimuthal7Tuple} issue. 

313 ''' 

314 pass 

315 

316 

317class Azimuthal7Tuple(_NamedTuple): 

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

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

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

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

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

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

324 ''' 

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

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

327 

328 def antipodal(self, azimuth=None): 

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

330 

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

332 (C{compass degrees360}). 

333 ''' 

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

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

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

337 

338 

339class Equidistant(_AzimuthalBase): 

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

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

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

343 

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

345 L{EquidistantGeodSolve} or L{EquidistantKarney} projection 

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

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

348 ''' 

349 if _FOR_DOCS: 

350 __init__ = _AzimuthalBase.__init__ 

351 

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

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

354 

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

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

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

358 

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

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

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

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

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

364 

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

366 

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

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

369 ''' 

370 def _k_t(c): 

371 k = _N_1_0 if c < 0 else _1_0 

372 t = fabs(c) < EPS1 

373 if t: 

374 a = acos(c) 

375 s = sin(a) 

376 if s: 

377 k = a / s 

378 return k, t 

379 

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

381 

382 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds): 

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

384 

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

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

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

388 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} forthe location and 

389 optional, additional B{C{LatLon}} keyword arguments, 

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

391 

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

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

394 

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

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

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

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

399 perpendicular to this. 

400 ''' 

401 def _c(c): 

402 return c if c > EPS else None 

403 

404 return self._reverse(x, y, _c, False, LatLon, **name_LatLon_kwds) 

405 

406 

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

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

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

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

411 

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

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

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

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

416 radius (C{meter}). 

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

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

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

420 

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

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

423 

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

425 

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

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

428 

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

430 ''' 

431 

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

433 if E is Equidistant: 

434 try: 

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

436 except ImportError: 

437 pass 

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

439 

440 

441class _AzimuthalGeodesic(_AzimuthalBase): 

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

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

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

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

446 ''' 

447 _mask = 0 

448 

449 @Property_RO 

450 def geodesic(self): # PYCHOK no cover 

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

452 self._notOverloaded() 

453 

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

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

456 ''' 

457 s = M 

458 if s is None: # reciprocal, azimuthal scale 

459 s = (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 t, _ = _name2__(name_LatLon_kwds, _or_nameof=self) 

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

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

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

471 L{EquidistantKarney} projection. 

472 ''' 

473 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, **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): 

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 C{B{name}=NN} 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, LatLon=None, **name_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 LatLon: Class to use (C{LatLon}) or C{None}. 

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

509 optional, additional B{C{LatLon}} keyword arguments, 

510 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_LatLon_kwds) if LatLon is None else \ 

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

526 

527 

528class EquidistantExact(_EquidistantBase): 

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

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

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

532 

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

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

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

536 clockwise from true North. 

537 

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

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

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

541 geographic and projected coordinates. 

542 ''' 

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

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

545 

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

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

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

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

550 radius (C{meter}). 

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

552 

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

554 ''' 

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

556 

557 if _FOR_DOCS: 

558 forward = _EquidistantBase.forward 

559 reverse = _EquidistantBase.reverse 

560 

561 @Property_RO 

562 def geodesic(self): 

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

564 ''' 

565 return self.datum.ellipsoid.geodesicx 

566 

567 

568class EquidistantGeodSolve(_EquidistantBase): 

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

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

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

572 I{for testing purposes only}. 

573 

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

575 ''' 

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

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

578 

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

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

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

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

583 radius (C{meter}). 

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

585 

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

587 ''' 

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

589 

590 if _FOR_DOCS: 

591 forward = _EquidistantBase.forward 

592 reverse = _EquidistantBase.reverse 

593 

594 @Property_RO 

595 def geodesic(self): 

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

597 ''' 

598 return self.datum.ellipsoid.geodsolve 

599 

600 

601class EquidistantKarney(_EquidistantBase): 

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

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

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

605 

606 @see: L{EquidistantExact}. 

607 ''' 

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

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

610 

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

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

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

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

615 radius (C{meter}). 

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

617 

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

619 not installed or not found. 

620 

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

622 ''' 

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

624 

625 if _FOR_DOCS: 

626 forward = _EquidistantBase.forward 

627 reverse = _EquidistantBase.reverse 

628 

629 @Property_RO 

630 def geodesic(self): 

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

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

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

634 package is installed. 

635 ''' 

636 return self.datum.ellipsoid.geodesic 

637 

638 

639_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve, 

640 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI 

641 

642 

643class Gnomonic(_AzimuthalBase): 

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

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

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

647 ''' 

648 if _FOR_DOCS: 

649 __init__ = _AzimuthalBase.__init__ 

650 

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

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

653 

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

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

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

657 

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

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

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

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

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

663 

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

665 ''' 

666 def _k_t(c): 

667 t = c > EPS 

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

669 return k, t 

670 

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

672 

673 def reverse(self, x, y, LatLon=None, **name_LatLon_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 LatLon: Class to use (C{LatLon}) or C{None}. 

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

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

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

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, _c, False, LatLon, **name_LatLon_kwds) 

696 

697 

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

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 C{B{name}=NN} 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) # PYCHOK types 

727 except ImportError: 

728 pass 

729 return G(lat0, lon0, datum=datum, **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): 

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

738 ''' 

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

740 

741 g = self.geodesic 

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

743 

744 def forward(self, lat, lon, raiser=True, **name): # 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 raiser: Do or don't throw an error (C{bool}) if 

751 the location lies over the horizon. 

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

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, _name__(name), M=M) 

781 t._iteraton = 0 

782 return t 

783 

784 def reverse(self, x, y, LatLon=None, **name_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 LatLon: Class to use (C{LatLon}) or C{None}. 

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

791 optional, additional B{C{LatLon}} keyword arguments, 

792 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, name_LatLon_kwds, M=r.M12) if LatLon is None else \ 

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

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

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 C{B{name}=NN} 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) 

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

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 C{B{name}=NN} 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) 

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

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 C{B{name}=NN} 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) 

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

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 C{B{name}=NN} 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, LatLon=None, **name_LatLon_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 LatLon: Class to use (C{LatLon}) or C{None}. 

983 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} for the location 

984 and optional, additional B{C{LatLon}} keyword 

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

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, _c, True, LatLon, **name_LatLon_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): 

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 C{B{name}=NN} 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, LatLon=None, **name_LatLon_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 LatLon: Class to use (C{LatLon}) or C{None}. 

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

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

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

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, _c, False, LatLon, **name_LatLon_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): 

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 C{B{name}=NN} 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__ # 'k0', name__=Stereographic.k0.fget 

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, LatLon=None, **name_LatLon_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 LatLon: Class to use (C{LatLon}) or C{None}. 

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

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

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

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 range C{[-90..90] degrees}, C{lon} in range 

1118 C{[-180..180] degrees} and C{azimuth} clockwise from true North. 

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

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

1121 ''' 

1122 def _c(c): 

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

1124 

1125 return self._reverse(x, y, _c, False, LatLon, **name_LatLon_kwds) 

1126 

1127 

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

1129 

1130# **) MIT License 

1131# 

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

1133# 

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

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

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

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

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

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

1140# 

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

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

1143# 

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

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

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

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

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

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

1150# OTHER DEALINGS IN THE SOFTWARE.