1

2# -*- coding: utf-8 -*-

3

4 u'''Equidistant, Equal-Area, and other Azimuthal projections.

5

6 Classes L{Equidistant}, L{EquidistantExact}, L{EquidistantGeodSolve},

7 L{EquidistantKarney}, L{Gnomonic}, L{GnomonicExact}, L{GnomonicKarney},

8 L{LambertEqualArea}, L{Orthographic} and L{Stereographic}, classes

9 L{AzimuthalError}, L{Azimuthal7Tuple} and functions L{equidistant}

10 and L{gnomonic}.

11

12 L{EquidistantExact} and L{GnomonicExact} are based on exact geodesic classes

13 L{GeodesicExact} and L{GeodesicLineExact}, Python versions of I{Charles Karney}'s

14 C++ original U{GeodesicExact<https://GeographicLib.SourceForge.io/C++/doc/

15 classGeographicLib_1_1GeodesicExact.html>}, respectively U{GeodesicLineExact

16<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicLineExact.html>}.

17

18 Using L{EquidistantGeodSolve} requires I{Karney}'s utility U{GeodSolve

19<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>} to be

20 executable and set in env variable C{PYGEODESY_GEODSOLVE}, see module

21 L{geodsolve} for more details.

22

23 L{EquidistantKarney} and L{GnomonicKarney} require I{Karney}'s Python package

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

25

26 Other 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,

28 spherical and ellipsoidal. The radius used for the latter is the ellipsoid's

29 I{mean radius of curvature} at the latitude of the projection center point. For

30 further justification, see the first paragraph under U{Snyder's FORMULAS FOR THE

31 ELLIPSOID, page 197<https://Pubs.USGS.gov/pp/1395/report.pdf>}.

32

33 Page 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

36 See also U{here<https://WikiPedia.org/wiki/Azimuthal_equidistant_projection>},

37 especially the U{Comparison of the Azimuthal equidistant projection and some

38 azimuthal projections centred on 90° N at the same scale, ordered by projection

39 altitude 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

43 from __future__ import division as _; del _ # PYCHOK semicolon

44

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

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

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

48 from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB, \

49 _xinstanceof

50 from pygeodesy.datums import _spherical_datum, _WGS84

51 from pygeodesy.errors import _ValueError, _xdatum, _xkwds

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

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

54# from pygeodesy.formy import antipode # _MODS

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

56 _SPACE_, _x_, _y_

57 from pygeodesy.karney import _norm180

58 from pygeodesy.latlonBase import _MODS, LatLonBase as _LLB

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

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

61 from pygeodesy.namedTuples import LatLon2Tuple, LatLon4Tuple

62 from pygeodesy.props import deprecated_Property_RO, Property_RO, \

63 property_doc_, _update_all

64 from pygeodesy.streprs import Fmt, _fstrLL0, unstr

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

66 Scalar, Scalar_

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

68 sincos2, sincos2d, sincos2d_

69

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

71

72__all__ = _ALL_LAZY.azimuthal

73__version__ = '24.11.24'

74

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

76_over_horizon_ = 'over horizon'

77_TRIPS = 21 # numit, 4 sufficient

78

79

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

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

82 '''

83 e = Easting( x=x)

84 n = Northing(y=y)

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

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

87

88

89 class _AzimuthalBase(_NamedBase):

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

91

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

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

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

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

96 '''

97 _datum = _WGS84 # L{Datum}

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

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

100

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

102 '''New azimuthal projection.

103

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

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

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

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

108 radius (C{meter}).

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

110

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

112

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

114 '''

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

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

117 if name:

118 self.name = name

119

120 if lat0 or lon0: # often both 0

121 self._reset(lat0, lon0)

122

123 @Property_RO

124 def datum(self):

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

126 '''

127 return self._datum

128

129 @Property_RO

130 def equatoradius(self):

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

132 '''

133 return self.datum.ellipsoid.a

134

135 a = equatoradius

136

137 @Property_RO

138 def flattening(self):

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

140 '''

141 return self.datum.ellipsoid.f

142

143 f = flattening

144

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

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

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

148

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

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

151 '''

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

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

154 s0, c0 = self._sc0

155

156 cb *= ca

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

158 if t:

159 r = k * self.radius

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

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

162 else: # 0 or 180

163 e = n = z = _0_0

164

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

166 name=self._name__(name))

167 return t

168

169 def _forwards(self, *lls):

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

171 '''

172 _fwd = self.forward

173 for ll in lls:

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

175

176 @Property_RO

177 def lat0(self):

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

179 '''

180 return self._latlon0.lat

181

182 @property

183 def latlon0(self):

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

185 '''

186 return self._latlon0

187

188 @latlon0.setter # PYCHOK setter!

189 def latlon0(self, latlon0):

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

191

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

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

194 '''

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

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

197 if hasattr(latlon0, _datum_):

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

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

200

201 @Property_RO

202 def lon0(self):

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

204 '''

205 return self._latlon0.lon

206

207 @deprecated_Property_RO

208 def majoradius(self): # PYCHOK no cover

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

210 return self.equatoradius

211

212 @Property_RO

213 def radius(self):

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

215 '''

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

217

218 def reset(self, lat0, lon0):

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

220

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

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

223

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

225 '''

226 _update_all(self) # zap caches

227 self._reset(lat0, lon0)

228

229 def _reset(self, lat0, lon0):

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

231 '''

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

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

234 self._sc0 = sincos2d(self.lat0)

235

236 def reverse(self, x, y, **name_LatLon_and_kwds):

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

238 self._notOverloaded(x, y, **name_LatLon_and_kwds) # PYCHOK no cover

239

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

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

242 '''

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

244

245 c = _c(r / self.radius)

246 if c is None:

247 lat, lon = self.latlon0

248 k, z = _1_0, _0_0

249 else:

250 s0, c0 = self._sc0

251 sc, cc = sincos2(c)

252 k = c / sc

253 s = s0 * cc

254 if r > EPS0:

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

256 lat = degrees(asin1(s))

257 if lea or fabs(c0) > EPS:

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

259 else:

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

261 lon = _norm180(self.lon0 + d)

262

263 if LatLon is None:

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

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

266 else:

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

268 return t

269

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

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

272 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne.

273 '''

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

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

276 return t, d

277

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

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

280 '''

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

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

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

284 _xinstanceof(B, LatLon=r)

285 return r

286

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

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

289

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

291

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

293 (C{str}).

294 '''

295 return _fstrLL0(self, prec, True)

296

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

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

299

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

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

302

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

304 '''

305 t = _fstrLL0(self, prec, False)

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

307

308

309 class AzimuthalError(_ValueError):

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

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

312 L{Azimuthal7Tuple} issue.

313 '''

314 pass

315

316

317 class Azimuthal7Tuple(_NamedTuple):

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

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

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

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

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

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

324 '''

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

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

327

328 def antipodal(self, azimuth=None):

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

330

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

332 (C{compass degrees360}).

333 '''

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

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

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

337

338

339 class Equidistant(_AzimuthalBase):

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

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

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

343

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

345 L{EquidistantGeodSolve} or L{EquidistantKarney} projection

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

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

348 '''

349 if _FOR_DOCS:

350 __init__ = _AzimuthalBase.__init__

351

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

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

354

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

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

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

358

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

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

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

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

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

364

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

366

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

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

369 '''

370 def _k_t(c):

371 k = _N_1_0 if c < 0 else _1_0

372 t = fabs(c) < EPS1

373 if t:

374 a = acos(c)

375 s = sin(a)

376 if s:

377 k = a / s

378 return k, t

379

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

381

382 def reverse(self, x, y, **name_LatLon_and_kwds):

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

384

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

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

387 @kwarg name_LatLon_and_kwds: Optional C{B{name}=NN} and class C{B{LatLon}=None}

388 to use and optionally, additional B{C{LatLon}} keyword arguments,

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

390

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

404

405

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

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

435 except ImportError:

436 pass

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

438

439

440 class _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 '''I{Must be overloaded}.'''

451 self._notOverloaded()

452

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

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

455 '''

456 s = M

457 if s is None: # reciprocal, azimuthal scale

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

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

462

463

464 class _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):

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

470 L{EquidistantKarney} projection.

471 '''

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

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

507 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} and optionally, additional

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

509

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

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

512

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

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

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

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

517 to this.

518 '''

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

520

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

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

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

524

525

526 class EquidistantExact(_EquidistantBase):

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

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

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

530

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

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

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

534 clockwise from true North.

535

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

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

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

539 geographic and projected coordinates.

540 '''

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

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

543

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

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

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

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

548 radius (C{meter}).

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

550

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

552 '''

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

554

555 if _FOR_DOCS:

556 forward = _EquidistantBase.forward

557 reverse = _EquidistantBase.reverse

558

559 @Property_RO

560 def geodesic(self):

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

562 '''

563 return self.datum.ellipsoid.geodesicx

564

565

566 class EquidistantGeodSolve(_EquidistantBase):

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

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

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

570 I{for testing purposes only}.

571

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

573 '''

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

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

576

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

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

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

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

581 radius (C{meter}).

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

583

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

585 '''

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

587

588 if _FOR_DOCS:

589 forward = _EquidistantBase.forward

590 reverse = _EquidistantBase.reverse

591

592 @Property_RO

593 def geodesic(self):

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

595 '''

596 return self.datum.ellipsoid.geodsolve

597

598

599 class EquidistantKarney(_EquidistantBase):

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

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

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

603

604 @see: L{EquidistantExact}.

605 '''

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

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

608

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

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

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

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

613 radius (C{meter}).

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

615

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

617 not installed or not found.

618

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

620 '''

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

622

623 if _FOR_DOCS:

624 forward = _EquidistantBase.forward

625 reverse = _EquidistantBase.reverse

626

627 @Property_RO

628 def geodesic(self):

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

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

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

632 package is installed.

633 '''

634 return self.datum.ellipsoid.geodesic

635

636

637_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve,

638 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI

639

640

641 class Gnomonic(_AzimuthalBase):

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

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

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

645 '''

646 if _FOR_DOCS:

647 __init__ = _AzimuthalBase.__init__

648

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

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

651

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

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

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

655

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

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

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

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

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

661

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

663 '''

664 def _k_t(c):

665 t = c > EPS

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

667 return k, t

668

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

670

671 def reverse(self, x, y, **name_LatLon_and_kwds):

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

673

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

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

676 @kwarg name_LatLon_and_kwds: Optional C{B{name}=NN} and class C{B{LatLon}=None}

677 for the location and optionally, additional B{C{LatLon}} keyword

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

679

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

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

682

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

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

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

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

687 perpendicular to this.

688 '''

689 def _c(c):

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

691

692 return self._reverse(x, y, _c, False, **name_LatLon_and_kwds)

693

694

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

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

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

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

699

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

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

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

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

704 radius (C{meter}).

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

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

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

708

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

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

711

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

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

714

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

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

717

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

719 '''

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

721 if G is Gnomonic:

722 try:

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

724 except ImportError:

725 pass

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

727

728

729 class _GnomonicBase(_AzimuthalGeodesic):

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

731 and L{GnomonicKarney}.

732 '''

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

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

735 '''

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

737

738 g = self.geodesic

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

740

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

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

743 and northing.

744

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

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

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

748 the location lies over the horizon.

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

750

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

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

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

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

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

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

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

758

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

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

761 '''

762 self._iteration = 0

763

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

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

766 M = r.M21

767 if M > EPS0:

768 q = r.m12 / M # .M12

769 e, n = sincos2d(r.azi1)

770 e *= q

771 n *= q

772 elif raiser: # PYCHOK no cover

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

774 else: # PYCHOK no cover

775 e = n = NAN

776

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

778 t._iteraton = 0

779 return t

780

781 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds): # PYCHOK signature

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

783

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

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

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

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

788 additional B{C{LatLon}} keyword arguments, ignored if C{B{LatLon} is None}.

789

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

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

792

793 @raise AzimuthalError: No convergence.

794

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

796 C{[-180..180] degrees}. The C{azimuth} is clockwise from true North. The

797 scale is C{1 / reciprocal**2} in C{radial} direction and C{1 / reciprocal}

798 in the direction perpendicular to this.

799 '''

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

801

802 d = a = self.equatoradius

803 s = a * atan1(q, a)

804 if q > a: # PYCHOK no cover

805 def _d(r, q):

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

807

808 q = _1_0 / q

809 else: # little == True

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

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

812

813 a *= _EPS_K

814 m = self._mask

815 g = self.geodesic

816

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

818 _S2 = Fsum(s).fsum2f_

819 _abs = fabs

820 for i in range(1, _TRIPS):

821 r = _P(s, outmask=m)

822 if _abs(d) < a:

823 break

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

825 else: # PYCHOK no cover

826 self._iteration = _TRIPS

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

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

829

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

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

832 t._iteration = self._iteration = i

833 return t

834

835

836 class GnomonicExact(_GnomonicBase):

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

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

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

840

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

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

843 '''

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

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

846

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

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

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

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

851 radius (C{meter}).

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

853

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

855 '''

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

857

858 if _FOR_DOCS:

859 forward = _GnomonicBase.forward

860 reverse = _GnomonicBase.reverse

861

862 @Property_RO

863 def geodesic(self):

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

865 '''

866 return self.datum.ellipsoid.geodesicx

867

868

869 class GnomonicGeodSolve(_GnomonicBase):

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

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

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

873 intended I{for testing purposes only}.

874

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

876 '''

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

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

879

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

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

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

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

884 radius (C{meter}).

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

886

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

888 '''

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

890

891 if _FOR_DOCS:

892 forward = _GnomonicBase.forward

893 reverse = _GnomonicBase.reverse

894

895 @Property_RO

896 def geodesic(self):

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

898 '''

899 return self.datum.ellipsoid.geodsolve

900

901

902 class GnomonicKarney(_GnomonicBase):

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

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

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

906

907 @see: L{GnomonicExact}.

908 '''

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

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

911

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

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

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

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

916 radius (C{meter}).

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

918

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

920 not installed or not found.

921

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

923 '''

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

925

926 if _FOR_DOCS:

927 forward = _GnomonicBase.forward

928 reverse = _GnomonicBase.reverse

929

930 @Property_RO

931 def geodesic(self):

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

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

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

935 package is installed.

936 '''

937 return self.datum.ellipsoid.geodesic

938

939

940 class LambertEqualArea(_AzimuthalBase):

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

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

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

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

945 '''

946 if _FOR_DOCS:

947 __init__ = _AzimuthalBase.__init__

948

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

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

951

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

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

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

955

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

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

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

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

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

961

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

963 '''

964 def _k_t(c):

965 c += _1_0

966 t = c > EPS0

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

968 return k, t

969

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

971

972 def reverse(self, x, y, **name_LatLon_and_kwds):

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

974

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

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

977 @kwarg name_LatLon_and_kwds: Optional C{B{name}=NN} and class C{B{LatLon}=None}

978 to use and optionally, additional B{C{LatLon}} keyword arguments,

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

980

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

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

983

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

985 range C{[-180..180] degrees}. The C{scale} of the projection is C{1}

986 in I{radial} direction, C{azimuth} clockwise from true North and is C{1

987 / reciprocal} in the direction perpendicular to this.

988 '''

989 def _c(c):

990 c *= _0_5

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

992

993 return self._reverse(x, y, _c, True, **name_LatLon_and_kwds)

994

995

996 class Orthographic(_AzimuthalBase):

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

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

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

1000 '''

1001 if _FOR_DOCS:

1002 __init__ = _AzimuthalBase.__init__

1003

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

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

1006

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

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

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

1010

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

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

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

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

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

1016

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

1018 '''

1019 def _k_t(c):

1020 return _1_0, (c >= 0)

1021

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

1023

1024 def reverse(self, x, y, **name_LatLon_and_kwds):

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

1026

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

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

1029 @kwarg name_LatLon_and_kwds: Optional C{B{name}=NN} and class C{B{LatLon}=None}

1030 to use and optionally, additional B{C{LatLon}} keyword arguments,

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

1032

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

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

1035

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

1037 range C{[-180..180] degrees}. The C{scale} of the projection is C{1}

1038 in I{radial} direction, C{azimuth} clockwise from true North and is C{1

1039 / reciprocal} in the direction perpendicular to this.

1040 '''

1041 def _c(c):

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

1043

1044 return self._reverse(x, y, _c, False, **name_LatLon_and_kwds)

1045

1046

1047 class Stereographic(_AzimuthalBase):

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

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

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

1051 '''

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

1053 _k02 = _2_0 # double ._k0

1054

1055 if _FOR_DOCS:

1056 __init__ = _AzimuthalBase.__init__

1057

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

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

1060

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

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

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

1064

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

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

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

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

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

1070

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

1072 '''

1073 def _k_t(c):

1074 c += _1_0

1075 t = isnon0(c)

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

1077 return k, t

1078

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

1080

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

1082 def k0(self):

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

1084 '''

1085 return self._k0

1086

1087 @k0.setter # PYCHOK setter!

1088 def k0(self, factor):

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

1090 '''

1091 n = Stereographic.k0.fget.__name__ # 'k0', name__=Stereographic.k0.fget

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

1093 self._k02 = self._k0 * _2_0

1094

1095 def reverse(self, x, y, **name_LatLon_and_kwds):

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

1097

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

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

1100 @kwarg name_LatLon_and_kwds: Optional C{B{name}=NN} and class C{B{LatLon}=None}

1101 to use and optionally, additional B{C{LatLon}} keyword arguments,

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

1103

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

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

1106

1107 @note: The C{lat} will be in range C{[-90..90] degrees}, C{lon} in range

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

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

1110 / reciprocal} in the direction perpendicular to this.

1111 '''

1112 def _c(c):

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

1114

1115 return self._reverse(x, y, _c, False, **name_LatLon_and_kwds)

1116

1117

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

1119

1120# **) MIT License

1121#

1122# Copyright (C) 2016-2025 -- mrJean1 at Gmail -- All Rights Reserved.

1123#

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

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

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

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

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

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

1130#

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

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

1133#

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

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

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

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

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

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

1140# OTHER DEALINGS IN THE SOFTWARE.