Coverage for pygeodesy/azimuthal.py: 98%

314 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-09-20 13: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 .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.09.07' 

75 

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

77_over_horizon_ = 'over horizon' 

78_TRIPS = 21 # numit, 4 sufficient 

79 

80 

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

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

83 ''' 

84 e = Easting( x=x) 

85 n = Northing(y=y) 

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

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

88 

89 

90class _AzimuthalBase(_NamedBase): 

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

92 

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

94 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 a = equatoradius 

137 

138 @Property_RO 

139 def flattening(self): 

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

141 ''' 

142 return self.datum.ellipsoid.f 

143 

144 f = flattening 

145 

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

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

148 ''' 

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

150 

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

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

153 ''' 

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

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

156 s0, c0 = self._sc0 

157 

158 cb *= ca 

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

160 if t: 

161 r = k * self.radius 

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

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

164 else: # 0 or 180 

165 e = n = z = _0_0 

166 

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

168 name=name or self.name) 

169 return t 

170 

171 def _forwards(self, *lls): 

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

173 ''' 

174 _fwd = self.forward 

175 for ll in lls: 

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

177 

178 @Property_RO 

179 def lat0(self): 

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

181 ''' 

182 return self._latlon0.lat 

183 

184 @property 

185 def latlon0(self): 

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

187 ''' 

188 return self._latlon0 

189 

190 @latlon0.setter # PYCHOK setter! 

191 def latlon0(self, latlon0): 

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

193 

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

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

196 ''' 

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

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

199 if hasattr(latlon0, _datum_): 

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

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

202 

203 @Property_RO 

204 def lon0(self): 

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

206 ''' 

207 return self._latlon0.lon 

208 

209 @deprecated_Property_RO 

210 def majoradius(self): # PYCHOK no cover 

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

212 return self.equatoradius 

213 

214 @Property_RO 

215 def radius(self): 

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

217 ''' 

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

219 

220 def reset(self, lat0, lon0): 

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

222 

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

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

225 

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

227 ''' 

228 _update_all(self) # zap caches 

229 self._reset(lat0, lon0) 

230 

231 def _reset(self, lat0, lon0): 

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

233 ''' 

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

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

236 self._sc0 = sincos2d(self.lat0) 

237 

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

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

240 ''' 

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

242 

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

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

245 ''' 

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

247 

248 c = _c(r / self.radius) 

249 if c is None: 

250 lat, lon = self.latlon0 

251 k, z = _1_0, _0_0 

252 else: 

253 s0, c0 = self._sc0 

254 sc, cc = sincos2(c) 

255 k = c / sc 

256 s = s0 * cc 

257 if r > EPS0: 

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

259 lat = degrees(asin1(s)) 

260 if lea or fabs(c0) > EPS: 

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

262 else: 

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

264 lon = _norm180(self.lon0 + d) 

265 

266 if LatLon is None: 

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

268 name=name or self.name) 

269 else: 

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

271 return t 

272 

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

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

275 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne. 

276 ''' 

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

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

279 return t, d 

280 

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

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

283 ''' 

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

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

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

287 _xinstanceof(B, LatLon=r) 

288 return r 

289 

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

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

292 

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

294 

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

296 (C{str}). 

297 ''' 

298 return _fstrLL0(self, prec, True) 

299 

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

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

302 

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

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

305 

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

307 ''' 

308 t = _fstrLL0(self, prec, False) 

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

310 

311 

312class AzimuthalError(_ValueError): 

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

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

315 L{Azimuthal7Tuple} issue. 

316 ''' 

317 pass 

318 

319 

320class Azimuthal7Tuple(_NamedTuple): 

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

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

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

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

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

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

327 ''' 

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

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

330 

331 def antipodal(self, azimuth=None): 

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

333 

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

335 (C{compass degrees360}). 

336 ''' 

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

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

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

340 

341 

342class Equidistant(_AzimuthalBase): 

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

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

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

346 

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

348 L{EquidistantGeodSolve} or L{EquidistantKarney} projection 

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

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

351 ''' 

352 if _FOR_DOCS: 

353 __init__ = _AzimuthalBase.__init__ 

354 

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

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

357 

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

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

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

361 

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

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

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

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

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

367 

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

369 

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

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

372 ''' 

373 def _k_t(c): 

374 k = _N_1_0 if c < 0 else _1_0 

375 t = fabs(c) < EPS1 

376 if t: 

377 a = acos(c) 

378 s = sin(a) 

379 if s: 

380 k = a / s 

381 return k, t 

382 

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

384 

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

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

387 

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

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

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

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

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

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

394 

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

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

397 

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

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

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

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

402 perpendicular to this. 

403 ''' 

404 def _c(c): 

405 return c if c > EPS else None 

406 

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

408 

409 

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

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

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

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

414 

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

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

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

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

419 radius (C{meter}). 

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

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

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

423 

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

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

426 

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

428 

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

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

431 

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

433 ''' 

434 

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

436 if E is Equidistant: 

437 try: 

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

439 except ImportError: 

440 pass 

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

442 

443 

444class _AzimuthalGeodesic(_AzimuthalBase): 

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

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

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

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

449 ''' 

450 _mask = 0 

451 

452 @Property_RO 

453 def geodesic(self): # PYCHOK no cover 

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

455 ''' 

456 notOverloaded(self) 

457 

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

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

460 ''' 

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

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

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

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

465 name=name or self.name) 

466 

467 

468class _EquidistantBase(_AzimuthalGeodesic): 

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

470 and L{EquidistantKarney}. 

471 ''' 

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

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

474 L{EquidistantKarney} projection. 

475 ''' 

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

477 

478 g = self.geodesic 

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

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

481 

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

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

484 

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

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

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

488 

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

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

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

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

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

494 

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

496 

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

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

499 ''' 

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

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

502 x, y = sincos2d(r.azi1) 

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

504 

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

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

507 

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

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

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

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

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

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

514 

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

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

517 

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

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

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

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

522 to this. 

523 ''' 

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

525 

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

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

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

529 

530 

531class EquidistantExact(_EquidistantBase): 

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

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

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

535 

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

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

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

539 clockwise from true North. 

540 

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

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

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

544 geographic and projected coordinates. 

545 ''' 

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

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

548 

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

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

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

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

553 radius (C{meter}). 

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

555 

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

557 ''' 

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

559 

560 if _FOR_DOCS: 

561 forward = _EquidistantBase.forward 

562 reverse = _EquidistantBase.reverse 

563 

564 @Property_RO 

565 def geodesic(self): 

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

567 ''' 

568 return self.datum.ellipsoid.geodesicx 

569 

570 

571class EquidistantGeodSolve(_EquidistantBase): 

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

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

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

575 I{for testing purposes only}. 

576 

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

578 ''' 

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

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

581 

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

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

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

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

586 radius (C{meter}). 

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

588 

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

590 ''' 

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

592 

593 if _FOR_DOCS: 

594 forward = _EquidistantBase.forward 

595 reverse = _EquidistantBase.reverse 

596 

597 @Property_RO 

598 def geodesic(self): 

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

600 ''' 

601 return self.datum.ellipsoid.geodsolve 

602 

603 

604class EquidistantKarney(_EquidistantBase): 

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

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

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

608 

609 @see: L{EquidistantExact}. 

610 ''' 

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

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

613 

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

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

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

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

618 radius (C{meter}). 

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

620 

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

622 not installed or not found. 

623 

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

625 ''' 

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

627 

628 if _FOR_DOCS: 

629 forward = _EquidistantBase.forward 

630 reverse = _EquidistantBase.reverse 

631 

632 @Property_RO 

633 def geodesic(self): 

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

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

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

637 package is installed. 

638 ''' 

639 return self.datum.ellipsoid.geodesic 

640 

641 

642_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve, 

643 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI 

644 

645 

646class Gnomonic(_AzimuthalBase): 

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

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

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

650 ''' 

651 if _FOR_DOCS: 

652 __init__ = _AzimuthalBase.__init__ 

653 

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

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

656 

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

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

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

660 

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

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

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

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

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

666 

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

668 ''' 

669 def _k_t(c): 

670 t = c > EPS 

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

672 return k, t 

673 

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

675 

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

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

678 

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

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

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

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

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

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

685 

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

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

688 

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

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

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

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

693 perpendicular to this. 

694 ''' 

695 def _c(c): 

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

697 

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

699 

700 

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

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

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

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

705 

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

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

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

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

710 radius (C{meter}). 

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

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

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

714 

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

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

717 

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

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

720 

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

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

723 

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

725 ''' 

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

727 if G is Gnomonic: 

728 try: 

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

730 except ImportError: 

731 pass 

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

733 

734 

735class _GnomonicBase(_AzimuthalGeodesic): 

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

737 and L{GnomonicKarney}. 

738 ''' 

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

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

741 ''' 

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

743 

744 g = self.geodesic 

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

746 

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

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

749 and northing. 

750 

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

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

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

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

755 the location lies over the horizon. 

756 

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

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

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

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

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

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

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

764 

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

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

767 ''' 

768 self._iteration = 0 

769 

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

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

772 M = r.M21 

773 if M > EPS0: 

774 q = r.m12 / M # .M12 

775 e, n = sincos2d(r.azi1) 

776 e *= q 

777 n *= q 

778 elif raiser: # PYCHOK no cover 

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

780 else: # PYCHOK no cover 

781 e = n = NAN 

782 

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

784 t._iteraton = 0 

785 return t 

786 

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

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

789 

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

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

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

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

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

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

796 

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

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

799 

800 @raise AzimuthalError: No convergence. 

801 

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

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

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

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

806 to this. 

807 ''' 

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

809 

810 d = a = self.equatoradius 

811 s = a * atan(q / a) 

812 if q > a: # PYCHOK no cover 

813 def _d(r, q): 

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

815 

816 q = _1_0 / q 

817 else: # little == True 

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

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

820 

821 a *= _EPS_K 

822 m = self._mask 

823 g = self.geodesic 

824 

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

826 _S2 = Fsum(s).fsum2_ 

827 for i in range(1, _TRIPS): 

828 r = _P(s, outmask=m) 

829 if fabs(d) < a: 

830 break 

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

832 else: # PYCHOK no cover 

833 self._iteration = _TRIPS 

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

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

836 

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

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

839 t._iteration = self._iteration = i 

840 return t 

841 

842 

843class GnomonicExact(_GnomonicBase): 

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

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

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

847 

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

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

850 ''' 

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

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

853 

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

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

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

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

858 radius (C{meter}). 

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

860 

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

862 ''' 

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

864 

865 if _FOR_DOCS: 

866 forward = _GnomonicBase.forward 

867 reverse = _GnomonicBase.reverse 

868 

869 @Property_RO 

870 def geodesic(self): 

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

872 ''' 

873 return self.datum.ellipsoid.geodesicx 

874 

875 

876class GnomonicGeodSolve(_GnomonicBase): 

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

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

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

880 intended I{for testing purposes only}. 

881 

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

883 ''' 

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

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

886 

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

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

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

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

891 radius (C{meter}). 

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

893 

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

895 ''' 

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

897 

898 if _FOR_DOCS: 

899 forward = _GnomonicBase.forward 

900 reverse = _GnomonicBase.reverse 

901 

902 @Property_RO 

903 def geodesic(self): 

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

905 ''' 

906 return self.datum.ellipsoid.geodsolve 

907 

908 

909class GnomonicKarney(_GnomonicBase): 

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

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

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

913 

914 @see: L{GnomonicExact}. 

915 ''' 

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

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

918 

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

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

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

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

923 radius (C{meter}). 

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

925 

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

927 not installed or not found. 

928 

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

930 ''' 

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

932 

933 if _FOR_DOCS: 

934 forward = _GnomonicBase.forward 

935 reverse = _GnomonicBase.reverse 

936 

937 @Property_RO 

938 def geodesic(self): 

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

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

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

942 package is installed. 

943 ''' 

944 return self.datum.ellipsoid.geodesic 

945 

946 

947class LambertEqualArea(_AzimuthalBase): 

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

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

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

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

952 ''' 

953 if _FOR_DOCS: 

954 __init__ = _AzimuthalBase.__init__ 

955 

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

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

958 

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

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

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

962 

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

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

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

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

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

968 

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

970 ''' 

971 def _k_t(c): 

972 c += _1_0 

973 t = c > EPS0 

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

975 return k, t 

976 

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

978 

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

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

981 

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

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

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

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

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

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

988 

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

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

991 

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

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

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

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

996 perpendicular to this. 

997 ''' 

998 def _c(c): 

999 c *= _0_5 

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

1001 

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

1003 

1004 

1005class Orthographic(_AzimuthalBase): 

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

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

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

1009 ''' 

1010 if _FOR_DOCS: 

1011 __init__ = _AzimuthalBase.__init__ 

1012 

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

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

1015 

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

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

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

1019 

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

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

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

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

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

1025 

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

1027 ''' 

1028 def _k_t(c): 

1029 return _1_0, (c >= 0) 

1030 

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

1032 

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

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

1035 

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

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

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

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

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

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

1042 

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

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

1045 

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

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

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

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

1050 perpendicular to this. 

1051 ''' 

1052 def _c(c): 

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

1054 

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

1056 

1057 

1058class Stereographic(_AzimuthalBase): 

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

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

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

1062 ''' 

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

1064 _k02 = _2_0 # double ._k0 

1065 

1066 if _FOR_DOCS: 

1067 __init__ = _AzimuthalBase.__init__ 

1068 

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

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

1071 

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

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

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

1075 

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

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

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

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

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

1081 

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

1083 ''' 

1084 def _k_t(c): 

1085 c += _1_0 

1086 t = isnon0(c) 

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

1088 return k, t 

1089 

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

1091 

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

1093 def k0(self): 

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

1095 ''' 

1096 return self._k0 

1097 

1098 @k0.setter # PYCHOK setter! 

1099 def k0(self, factor): 

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

1101 ''' 

1102 n = Stereographic.k0.fget.__name__ 

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

1104 self._k02 = self._k0 * _2_0 

1105 

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

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

1108 

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

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

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

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

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

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

1115 

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

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

1118 

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

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

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

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

1123 perpendicular to this. 

1124 ''' 

1125 def _c(c): 

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

1127 

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

1129 

1130 

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

1132 

1133# **) MIT License 

1134# 

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

1136# 

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

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

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

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

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

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

1143# 

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

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

1146# 

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

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

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

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

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

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

1153# OTHER DEALINGS IN THE SOFTWARE.