Coverage for pygeodesy/azimuthal.py: 98%

312 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-23 12:10 -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.07.10' 

75 

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

77_over_horizon_ = 'over horizon' 

78_TRIPS = 21 # numit, 4 sufficient 

79 

80 

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

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

83 ''' 

84 e = Easting( x=x) 

85 n = Northing(y=y) 

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

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

88 

89 

90class _AzimuthalBase(_NamedBase): 

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

92 

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

94 html/classGeographicLib_1_1AzimuthalEquidistant.html>} and U{Gnomonic 

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

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

97 ''' 

98 _datum = _WGS84 # L{Datum} 

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

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

101 

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

103 '''New azimuthal projection. 

104 

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

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

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

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

109 radius (C{meter}). 

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

111 

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

113 

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

115 ''' 

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

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

118 if name: 

119 self.name = name 

120 

121 if lat0 or lon0: # often both 0 

122 self._reset(lat0, lon0) 

123 

124 @Property_RO 

125 def datum(self): 

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

127 ''' 

128 return self._datum 

129 

130 @Property_RO 

131 def equatoradius(self): 

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

133 ''' 

134 return self.datum.ellipsoid.a 

135 

136 @Property_RO 

137 def flattening(self): 

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

139 ''' 

140 return self.datum.ellipsoid.f 

141 

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

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

144 ''' 

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

146 

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

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

149 ''' 

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

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

152 s0, c0 = self._sc0 

153 

154 cb *= ca 

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

156 if t: 

157 r = k * self.radius 

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

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

160 else: # 0 or 180 

161 e = n = z = _0_0 

162 

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

164 name=name or self.name) 

165 return t 

166 

167 def _forwards(self, *lls): 

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

169 ''' 

170 _fwd = self.forward 

171 for ll in lls: 

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

173 

174 @Property_RO 

175 def lat0(self): 

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

177 ''' 

178 return self._latlon0.lat 

179 

180 @property 

181 def latlon0(self): 

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

183 ''' 

184 return self._latlon0 

185 

186 @latlon0.setter # PYCHOK setter! 

187 def latlon0(self, latlon0): 

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

189 

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

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

192 ''' 

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

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

195 if hasattr(latlon0, _datum_): 

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

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

198 

199 @Property_RO 

200 def lon0(self): 

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

202 ''' 

203 return self._latlon0.lon 

204 

205 @deprecated_Property_RO 

206 def majoradius(self): # PYCHOK no cover 

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

208 return self.equatoradius 

209 

210 @Property_RO 

211 def radius(self): 

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

213 ''' 

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

215 

216 def reset(self, lat0, lon0): 

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

218 

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

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

221 

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

223 ''' 

224 _update_all(self) # zap caches 

225 self._reset(lat0, lon0) 

226 

227 def _reset(self, lat0, lon0): 

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

229 ''' 

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

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

232 self._sc0 = sincos2d(self.lat0) 

233 

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

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

236 ''' 

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

238 

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

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

241 ''' 

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

243 

244 c = _c(r / self.radius) 

245 if c is None: 

246 lat, lon = self.latlon0 

247 k, z = _1_0, _0_0 

248 else: 

249 s0, c0 = self._sc0 

250 sc, cc = sincos2(c) 

251 k = c / sc 

252 s = s0 * cc 

253 if r > EPS0: 

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

255 lat = degrees(asin1(s)) 

256 if lea or fabs(c0) > EPS: 

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

258 else: 

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

260 lon = _norm180(self.lon0 + d) 

261 

262 if LatLon is None: 

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

264 name=name or self.name) 

265 else: 

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

267 return t 

268 

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

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

271 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne. 

272 ''' 

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

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

275 return t, d 

276 

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

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

279 ''' 

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

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

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

283 _xinstanceof(B, LatLon=r) 

284 return r 

285 

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

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

288 

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

290 

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

292 (C{str}). 

293 ''' 

294 return _fstrLL0(self, prec, True) 

295 

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

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

298 

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

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

301 

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

303 ''' 

304 t = _fstrLL0(self, prec, False) 

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

306 

307 

308class AzimuthalError(_ValueError): 

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

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

311 L{Azimuthal7Tuple} issue. 

312 ''' 

313 pass 

314 

315 

316class Azimuthal7Tuple(_NamedTuple): 

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

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

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

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

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

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

323 ''' 

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

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

326 

327 def antipodal(self, azimuth=None): 

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

329 

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

331 (C{compass degrees360}). 

332 ''' 

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

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

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

336 

337 

338class Equidistant(_AzimuthalBase): 

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

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

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

342 

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

344 L{EquidistantGeodSolve} or L{EquidistantKarney} projection 

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

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

347 ''' 

348 if _FOR_DOCS: 

349 __init__ = _AzimuthalBase.__init__ 

350 

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

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

353 

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

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

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

357 

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

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

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

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

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

363 

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

365 

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

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

368 ''' 

369 def _k_t(c): 

370 k = _N_1_0 if c < 0 else _1_0 

371 t = fabs(c) < EPS1 

372 if t: 

373 a = acos(c) 

374 s = sin(a) 

375 if s: 

376 k = a / s 

377 return k, t 

378 

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

380 

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

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

383 

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

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

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

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

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

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

390 

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

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

393 

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

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

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

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

398 perpendicular to this. 

399 ''' 

400 def _c(c): 

401 return c if c > EPS else None 

402 

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

404 

405 

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

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

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

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

410 

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

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

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

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

415 radius (C{meter}). 

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

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

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

419 

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

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

422 

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

424 

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

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

427 

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

429 ''' 

430 

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

432 if E is Equidistant: 

433 try: 

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

435 except ImportError: 

436 pass 

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

438 

439 

440class _AzimuthalGeodesic(_AzimuthalBase): 

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

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

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

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

445 ''' 

446 _mask = 0 

447 

448 @Property_RO 

449 def geodesic(self): # PYCHOK no cover 

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

451 ''' 

452 notOverloaded(self) 

453 

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

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

456 ''' 

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

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

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

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

461 name=name or self.name) 

462 

463 

464class _EquidistantBase(_AzimuthalGeodesic): 

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

466 and L{EquidistantKarney}. 

467 ''' 

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

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

470 L{EquidistantKarney} projection. 

471 ''' 

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

473 

474 g = self.geodesic 

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

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

477 

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

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

480 

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

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

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

484 

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

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

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

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

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

490 

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

492 

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

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

495 ''' 

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

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

498 x, y = sincos2d(r.azi1) 

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

500 

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

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

503 

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

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

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

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

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

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

510 

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

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

513 

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

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

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

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

518 to this. 

519 ''' 

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

521 

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

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

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

525 

526 

527class EquidistantExact(_EquidistantBase): 

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

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

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

531 

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

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

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

535 clockwise from true North. 

536 

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

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

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

540 geographic and projected coordinates. 

541 ''' 

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

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

544 

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

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

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

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

549 radius (C{meter}). 

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

551 

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

553 ''' 

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

555 

556 if _FOR_DOCS: 

557 forward = _EquidistantBase.forward 

558 reverse = _EquidistantBase.reverse 

559 

560 @Property_RO 

561 def geodesic(self): 

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

563 ''' 

564 return self.datum.ellipsoid.geodesicx 

565 

566 

567class EquidistantGeodSolve(_EquidistantBase): 

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

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

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

571 I{for testing purposes only}. 

572 

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

574 ''' 

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

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

577 

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

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

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

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

582 radius (C{meter}). 

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

584 

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

586 ''' 

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

588 

589 if _FOR_DOCS: 

590 forward = _EquidistantBase.forward 

591 reverse = _EquidistantBase.reverse 

592 

593 @Property_RO 

594 def geodesic(self): 

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

596 ''' 

597 return self.datum.ellipsoid.geodsolve 

598 

599 

600class EquidistantKarney(_EquidistantBase): 

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

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

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

604 

605 @see: L{EquidistantExact}. 

606 ''' 

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

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

609 

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

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

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

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

614 radius (C{meter}). 

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

616 

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

618 not installed or not found. 

619 

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

621 ''' 

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

623 

624 if _FOR_DOCS: 

625 forward = _EquidistantBase.forward 

626 reverse = _EquidistantBase.reverse 

627 

628 @Property_RO 

629 def geodesic(self): 

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

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

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

633 package is installed. 

634 ''' 

635 return self.datum.ellipsoid.geodesic 

636 

637 

638_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve, 

639 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI 

640 

641 

642class Gnomonic(_AzimuthalBase): 

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

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

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

646 ''' 

647 if _FOR_DOCS: 

648 __init__ = _AzimuthalBase.__init__ 

649 

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

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

652 

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

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

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

656 

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

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

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

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

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

662 

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

664 ''' 

665 def _k_t(c): 

666 t = c > EPS 

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

668 return k, t 

669 

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

671 

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

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

674 

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

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

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

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

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

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

681 

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

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

684 

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

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

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

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

689 perpendicular to this. 

690 ''' 

691 def _c(c): 

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

693 

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

695 

696 

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

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

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

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

701 

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

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

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

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

706 radius (C{meter}). 

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

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

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

710 

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

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

713 

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

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

716 

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

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

719 

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

721 ''' 

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

723 if G is Gnomonic: 

724 try: 

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

726 except ImportError: 

727 pass 

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

729 

730 

731class _GnomonicBase(_AzimuthalGeodesic): 

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

733 and L{GnomonicKarney}. 

734 ''' 

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

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

737 ''' 

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

739 

740 g = self.geodesic 

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

742 

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

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

745 and northing. 

746 

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

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

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

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

751 the location lies over the horizon. 

752 

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

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

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

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

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

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

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

760 

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

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

763 ''' 

764 self._iteration = 0 

765 

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

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

768 M = r.M21 

769 if M > EPS0: 

770 q = r.m12 / M # .M12 

771 e, n = sincos2d(r.azi1) 

772 e *= q 

773 n *= q 

774 elif raiser: # PYCHOK no cover 

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

776 else: # PYCHOK no cover 

777 e = n = NAN 

778 

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

780 t._iteraton = 0 

781 return t 

782 

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

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

785 

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

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

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

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

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

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

792 

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

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

795 

796 @raise AzimuthalError: No convergence. 

797 

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

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

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

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

802 to this. 

803 ''' 

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

805 

806 d = a = self.equatoradius 

807 s = a * atan(q / a) 

808 if q > a: # PYCHOK no cover 

809 def _d(r, q): 

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

811 

812 q = _1_0 / q 

813 else: # little == True 

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

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

816 

817 a *= _EPS_K 

818 m = self._mask 

819 g = self.geodesic 

820 

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

822 _S2 = Fsum(s).fsum2_ 

823 for i in range(1, _TRIPS): 

824 r = _P(s, outmask=m) 

825 if fabs(d) < a: 

826 break 

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

828 else: # PYCHOK no cover 

829 self._iteration = _TRIPS 

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

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

832 

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

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

835 t._iteration = self._iteration = i 

836 return t 

837 

838 

839class GnomonicExact(_GnomonicBase): 

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

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

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

843 

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

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

846 ''' 

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

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

849 

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

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

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

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

854 radius (C{meter}). 

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

856 

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

858 ''' 

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

860 

861 if _FOR_DOCS: 

862 forward = _GnomonicBase.forward 

863 reverse = _GnomonicBase.reverse 

864 

865 @Property_RO 

866 def geodesic(self): 

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

868 ''' 

869 return self.datum.ellipsoid.geodesicx 

870 

871 

872class GnomonicGeodSolve(_GnomonicBase): 

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

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

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

876 intended I{for testing purposes only}. 

877 

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

879 ''' 

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

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

882 

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

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

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

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

887 radius (C{meter}). 

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

889 

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

891 ''' 

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

893 

894 if _FOR_DOCS: 

895 forward = _GnomonicBase.forward 

896 reverse = _GnomonicBase.reverse 

897 

898 @Property_RO 

899 def geodesic(self): 

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

901 ''' 

902 return self.datum.ellipsoid.geodsolve 

903 

904 

905class GnomonicKarney(_GnomonicBase): 

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

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

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

909 

910 @see: L{GnomonicExact}. 

911 ''' 

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

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

914 

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

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

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

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

919 radius (C{meter}). 

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

921 

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

923 not installed or not found. 

924 

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

926 ''' 

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

928 

929 if _FOR_DOCS: 

930 forward = _GnomonicBase.forward 

931 reverse = _GnomonicBase.reverse 

932 

933 @Property_RO 

934 def geodesic(self): 

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

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

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

938 package is installed. 

939 ''' 

940 return self.datum.ellipsoid.geodesic 

941 

942 

943class LambertEqualArea(_AzimuthalBase): 

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

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

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

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

948 ''' 

949 if _FOR_DOCS: 

950 __init__ = _AzimuthalBase.__init__ 

951 

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

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

954 

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

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

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

958 

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

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

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

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

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

964 

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

966 ''' 

967 def _k_t(c): 

968 c += _1_0 

969 t = c > EPS0 

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

971 return k, t 

972 

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

974 

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

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

977 

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

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

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

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

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

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

984 

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

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

987 

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

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

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

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

992 perpendicular to this. 

993 ''' 

994 def _c(c): 

995 c *= _0_5 

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

997 

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

999 

1000 

1001class Orthographic(_AzimuthalBase): 

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

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

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

1005 ''' 

1006 if _FOR_DOCS: 

1007 __init__ = _AzimuthalBase.__init__ 

1008 

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

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

1011 

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

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

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

1015 

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

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

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

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

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

1021 

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

1023 ''' 

1024 def _k_t(c): 

1025 return _1_0, (c >= 0) 

1026 

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

1028 

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

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

1031 

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

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

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

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

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

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

1038 

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

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

1041 

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

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

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

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

1046 perpendicular to this. 

1047 ''' 

1048 def _c(c): 

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

1050 

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

1052 

1053 

1054class Stereographic(_AzimuthalBase): 

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

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

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

1058 ''' 

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

1060 _k02 = _2_0 # double ._k0 

1061 

1062 if _FOR_DOCS: 

1063 __init__ = _AzimuthalBase.__init__ 

1064 

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

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

1067 

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

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

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

1071 

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

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

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

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

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

1077 

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

1079 ''' 

1080 def _k_t(c): 

1081 c += _1_0 

1082 t = isnon0(c) 

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

1084 return k, t 

1085 

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

1087 

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

1089 def k0(self): 

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

1091 ''' 

1092 return self._k0 

1093 

1094 @k0.setter # PYCHOK setter! 

1095 def k0(self, factor): 

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

1097 ''' 

1098 n = Stereographic.k0.fget.__name__ 

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

1100 self._k02 = self._k0 * _2_0 

1101 

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

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

1104 

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

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

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

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

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

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

1111 

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

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

1114 

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

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

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

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

1119 perpendicular to this. 

1120 ''' 

1121 def _c(c): 

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

1123 

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

1125 

1126 

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

1128 

1129# **) MIT License 

1130# 

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

1132# 

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

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

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

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

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

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

1139# 

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

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

1142# 

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

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

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

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

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

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

1149# OTHER DEALINGS IN THE SOFTWARE.