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 .datums

46 from 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

49 from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB, \

50 _spherical_datum, _WGS84, _xinstanceof

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

52 from pygeodesy.errors import _ValueError, _xdatum, _xkwds

53 from pygeodesy.fmath import euclid, Fsum, hypot

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

55# from pygeodesy.formy import antipode # from latlonBase

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

57 _scale_, _SPACE_, _x_, _y_

58 from pygeodesy.karney import _norm180

59 from pygeodesy.latlonBase import antipode, LatLonBase as _LLB

60 from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS

61 from pygeodesy.named import _NamedBase, _NamedTuple, notOverloaded, _Pass

62 from pygeodesy.namedTuples import LatLon2Tuple, LatLon4Tuple

63 from pygeodesy.props import deprecated_Property_RO, Property_RO, \

64 property_doc_, _update_all

65 from pygeodesy.streprs import Fmt, _fstrLL0, unstr

66 from pygeodesy.units import Bearing, Easting, Lat_, Lon_, \

67 Northing, Scalar, Scalar_

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

69 sincos2d, sincos2d_

70

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

72

73__all__ = _ALL_LAZY.azimuthal

74__version__ = '23.03.21'

75

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

77_over_horizon_ = 'over horizon'

78_TRIPS = 21 # numit, 4 sufficient

79

80

81 def _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

90 class _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

308 class 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

316 class 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

338 class 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

406 def 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

440 class _AzimuthalGeodesic(_AzimuthalBase):

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

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

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

444 '''

445 _mask = 0

446

447 @Property_RO

448 def geodesic(self): # PYCHOK no cover

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

450 '''

451 notOverloaded(self)

452

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

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

455 '''

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

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

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

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

460 name=name or self.name)

461

462

463 class _EquidistantBase(_AzimuthalGeodesic):

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

465 and L{EquidistantKarney}.

466 '''

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

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

469 L{EquidistantKarney} projection.

470 '''

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

472

473 g = self.geodesic

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

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

476

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

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

479

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

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

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

483

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

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

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

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

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

489

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

491

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

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

494 '''

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

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

497 x, y = sincos2d(r.azi1)

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

499

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

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

502

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

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

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

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

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

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

509

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

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

512

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

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

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

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

517 to this.

518 '''

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

520

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

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

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

524

525

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

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

543

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

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

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

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

548 radius (C{meter}).

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

550

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

552 '''

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

554

555 if _FOR_DOCS:

556 forward = _EquidistantBase.forward

557 reverse = _EquidistantBase.reverse

558

559 @Property_RO

560 def geodesic(self):

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

562 '''

563 return self.datum.ellipsoid.geodesicx

564

565

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

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

576

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

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

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

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

581 radius (C{meter}).

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

583

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

585 '''

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

587

588 if _FOR_DOCS:

589 forward = _EquidistantBase.forward

590 reverse = _EquidistantBase.reverse

591

592 @Property_RO

593 def geodesic(self):

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

595 '''

596 return self.datum.ellipsoid.geodsolve

597

598

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

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

608

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

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

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

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

613 radius (C{meter}).

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

615

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

617 not installed or not found.

618

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

620 '''

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

622

623 if _FOR_DOCS:

624 forward = _EquidistantBase.forward

625 reverse = _EquidistantBase.reverse

626

627 @Property_RO

628 def geodesic(self):

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

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

631 provided package U{geographiclib

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

633 '''

634 return self.datum.ellipsoid.geodesic

635

636

637_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve,

638 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI

639

640

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

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

651

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

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

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

655

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

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

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

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

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

661

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

663 '''

664 def _k_t(c):

665 t = c > EPS

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

667 return k, t

668

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

670

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

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

673

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

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

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

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

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

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

680

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

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

683

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

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

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

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

688 perpendicular to this.

689 '''

690 def _c(c):

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

692

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

694

695

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

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

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

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

700

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

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

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

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

705 radius (C{meter}).

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

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

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

709

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

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

712

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

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

715

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

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

718

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

720 '''

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

722 if G is Gnomonic:

723 try:

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

725 except ImportError:

726 pass

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

728

729

730 class _GnomonicBase(_AzimuthalGeodesic):

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

732 and L{GnomonicKarney}.

733 '''

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

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

736 '''

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

738

739 g = self.geodesic

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

741

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

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

744 and northing.

745

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

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

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

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

750 the location lies over the horizon.

751

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

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

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

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

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

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

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

759

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

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

762 '''

763 self._iteration = 0

764

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

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

767 M = r.M21

768 if M > EPS0:

769 q = r.m12 / M # .M12

770 e, n = sincos2d(r.azi1)

771 e *= q

772 n *= q

773 elif raiser: # PYCHOK no cover

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

775 else: # PYCHOK no cover

776 e = n = NAN

777

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

779 t._iteraton = 0

780 return t

781

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

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

784

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

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

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

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

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

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

791

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

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

794

795 @raise AzimuthalError: No convergence.

796

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

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

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

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

801 to this.

802 '''

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

804

805 d = a = self.equatoradius

806 s = a * atan(q / a)

807 if q > a: # PYCHOK no cover

808 def _d(r, q):

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

810

811 q = _1_0 / q

812 else: # little == True

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

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

815

816 a *= _EPS_K

817 m = self._mask

818 g = self.geodesic

819

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

821 _S2 = Fsum(s).fsum2_

822 for i in range(1, _TRIPS):

823 r = P(s, m)

824 if fabs(d) < a:

825 break

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

827 else: # PYCHOK no cover

828 self._iteration = _TRIPS

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

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

831

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

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

834 t._iteration = self._iteration = i

835 return t

836

837

838 class GnomonicExact(_GnomonicBase):

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

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

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

842

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

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

845 '''

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

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

848

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

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

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

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

853 radius (C{meter}).

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

855

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

857 '''

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

859

860 if _FOR_DOCS:

861 forward = _GnomonicBase.forward

862 reverse = _GnomonicBase.reverse

863

864 @Property_RO

865 def geodesic(self):

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

867 '''

868 return self.datum.ellipsoid.geodesicx

869

870

871 class GnomonicGeodSolve(_GnomonicBase):

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

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

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

875 intended I{for testing purposes only}.

876

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

878 '''

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

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

881

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

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

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

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

886 radius (C{meter}).

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

888

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

890 '''

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

892

893 if _FOR_DOCS:

894 forward = _GnomonicBase.forward

895 reverse = _GnomonicBase.reverse

896

897 @Property_RO

898 def geodesic(self):

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

900 '''

901 return self.datum.ellipsoid.geodsolve

902

903

904 class GnomonicKarney(_GnomonicBase):

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

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

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

908

909 @see: L{GnomonicExact}.

910 '''

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

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

913

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

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

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

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

918 radius (C{meter}).

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

920

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

922 not installed or not found.

923

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

925 '''

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

927

928 if _FOR_DOCS:

929 forward = _GnomonicBase.forward

930 reverse = _GnomonicBase.reverse

931

932 @Property_RO

933 def geodesic(self):

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

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

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

937 '''

938 return self.datum.ellipsoid.geodesic

939

940

941 class LambertEqualArea(_AzimuthalBase):

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

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

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

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

946 '''

947 if _FOR_DOCS:

948 __init__ = _AzimuthalBase.__init__

949

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

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

952

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

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

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

956

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

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

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

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

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

962

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

964 '''

965 def _k_t(c):

966 c += _1_0

967 t = c > EPS0

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

969 return k, t

970

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

972

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

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

975

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

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

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

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

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

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

982

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

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

985

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

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

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

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

990 perpendicular to this.

991 '''

992 def _c(c):

993 c *= _0_5

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

995

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

997

998

999 class Orthographic(_AzimuthalBase):

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

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

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

1003 '''

1004 if _FOR_DOCS:

1005 __init__ = _AzimuthalBase.__init__

1006

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

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

1009

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

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

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

1013

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

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

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

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

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

1019

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

1021 '''

1022 def _k_t(c):

1023 return _1_0, (c >= 0)

1024

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

1026

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

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

1029

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

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

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

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

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

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

1036

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

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

1039

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

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

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

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

1044 perpendicular to this.

1045 '''

1046 def _c(c):

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

1048

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

1050

1051

1052 class Stereographic(_AzimuthalBase):

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

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

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

1056 '''

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

1058 _k02 = _2_0 # double ._k0

1059

1060 if _FOR_DOCS:

1061 __init__ = _AzimuthalBase.__init__

1062

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

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

1065

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

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

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

1069

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

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

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

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

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

1075

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

1077 '''

1078 def _k_t(c):

1079 c += _1_0

1080 t = isnon0(c)

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

1082 return k, t

1083

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

1085

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

1087 def k0(self):

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

1089 '''

1090 return self._k0

1091

1092 @k0.setter # PYCHOK setter!

1093 def k0(self, factor):

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

1095 '''

1096 n = Stereographic.k0.fget.__name__

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

1098 self._k02 = self._k0 * _2_0

1099

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

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

1102

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

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

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

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

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

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

1109

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

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

1112

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

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

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

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

1117 perpendicular to this.

1118 '''

1119 def _c(c):

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

1121

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

1123

1124

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

1126

1127# **) MIT License

1128#

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

1130#

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

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

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

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

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

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

1137#

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

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

1140#

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

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

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

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

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

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

1147# OTHER DEALINGS IN THE SOFTWARE.