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.07.10'

75

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

77_over_horizon_ = 'over horizon'

78_TRIPS = 21 # numit, 4 sufficient

79

80

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 the

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

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

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

445 '''

446 _mask = 0

447

448 @Property_RO

449 def geodesic(self): # PYCHOK no cover

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

451 '''

452 notOverloaded(self)

453

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

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

456 '''

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

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

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

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

461 name=name or self.name)

462

463

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

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

470 L{EquidistantKarney} projection.

471 '''

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

473

474 g = self.geodesic

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

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

477

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

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

480

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

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

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

484

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

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

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

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

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

490

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

492

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

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

495 '''

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

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

498 x, y = sincos2d(r.azi1)

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

500

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

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

503

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

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

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

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

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

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

510

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

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

513

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

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

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

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

518 to this.

519 '''

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

521

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

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

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

525

526

527 class EquidistantExact(_EquidistantBase):

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

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

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

531

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

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

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

535 clockwise from true North.

536

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

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

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

540 geographic and projected coordinates.

541 '''

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

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

544

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

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

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

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

549 radius (C{meter}).

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

551

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

553 '''

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

555

556 if _FOR_DOCS:

557 forward = _EquidistantBase.forward

558 reverse = _EquidistantBase.reverse

559

560 @Property_RO

561 def geodesic(self):

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

563 '''

564 return self.datum.ellipsoid.geodesicx

565

566

567 class EquidistantGeodSolve(_EquidistantBase):

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

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

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

571 I{for testing purposes only}.

572

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

574 '''

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

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

577

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

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

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

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

582 radius (C{meter}).

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

584

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

586 '''

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

588

589 if _FOR_DOCS:

590 forward = _EquidistantBase.forward

591 reverse = _EquidistantBase.reverse

592

593 @Property_RO

594 def geodesic(self):

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

596 '''

597 return self.datum.ellipsoid.geodsolve

598

599

600 class EquidistantKarney(_EquidistantBase):

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

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

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

604

605 @see: L{EquidistantExact}.

606 '''

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

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

609

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

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

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

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

614 radius (C{meter}).

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

616

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

618 not installed or not found.

619

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

621 '''

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

623

624 if _FOR_DOCS:

625 forward = _EquidistantBase.forward

626 reverse = _EquidistantBase.reverse

627

628 @Property_RO

629 def geodesic(self):

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

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

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

633 package is installed.

634 '''

635 return self.datum.ellipsoid.geodesic

636

637

638_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve,

639 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI

640

641

642 class Gnomonic(_AzimuthalBase):

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

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

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

646 '''

647 if _FOR_DOCS:

648 __init__ = _AzimuthalBase.__init__

649

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

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

652

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

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

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

656

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

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

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

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

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

662

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

664 '''

665 def _k_t(c):

666 t = c > EPS

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

668 return k, t

669

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

671

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

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

674

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

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

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

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

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

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

681

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

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

684

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

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

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

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

689 perpendicular to this.

690 '''

691 def _c(c):

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

693

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

695

696

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

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

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

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

701

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

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

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

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

706 radius (C{meter}).

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

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

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

710

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

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

713

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

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

716

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

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

719

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

721 '''

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

723 if G is Gnomonic:

724 try:

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

726 except ImportError:

727 pass

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

729

730

731 class _GnomonicBase(_AzimuthalGeodesic):

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

733 and L{GnomonicKarney}.

734 '''

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

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

737 '''

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

739

740 g = self.geodesic

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

742

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

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

745 and northing.

746

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

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

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

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

751 the location lies over the horizon.

752

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

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

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

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

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

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

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

760

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

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

763 '''

764 self._iteration = 0

765

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

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

768 M = r.M21

769 if M > EPS0:

770 q = r.m12 / M # .M12

771 e, n = sincos2d(r.azi1)

772 e *= q

773 n *= q

774 elif raiser: # PYCHOK no cover

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

776 else: # PYCHOK no cover

777 e = n = NAN

778

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

780 t._iteraton = 0

781 return t

782

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

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

785

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

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

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

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

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

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

792

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

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

795

796 @raise AzimuthalError: No convergence.

797

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

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

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

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

802 to this.

803 '''

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

805

806 d = a = self.equatoradius

807 s = a * atan(q / a)

808 if q > a: # PYCHOK no cover

809 def _d(r, q):

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

811

812 q = _1_0 / q

813 else: # little == True

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

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

816

817 a *= _EPS_K

818 m = self._mask

819 g = self.geodesic

820

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

822 _S2 = Fsum(s).fsum2_

823 for i in range(1, _TRIPS):

824 r = _P(s, outmask=m)

825 if fabs(d) < a:

826 break

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

828 else: # PYCHOK no cover

829 self._iteration = _TRIPS

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

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

832

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

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

835 t._iteration = self._iteration = i

836 return t

837

838

839 class GnomonicExact(_GnomonicBase):

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

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

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

843

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

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

846 '''

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

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

849

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

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

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

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

854 radius (C{meter}).

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

856

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

858 '''

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

860

861 if _FOR_DOCS:

862 forward = _GnomonicBase.forward

863 reverse = _GnomonicBase.reverse

864

865 @Property_RO

866 def geodesic(self):

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

868 '''

869 return self.datum.ellipsoid.geodesicx

870

871

872 class GnomonicGeodSolve(_GnomonicBase):

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

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

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

876 intended I{for testing purposes only}.

877

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

879 '''

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

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

882

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

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

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

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

887 radius (C{meter}).

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

889

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

891 '''

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

893

894 if _FOR_DOCS:

895 forward = _GnomonicBase.forward

896 reverse = _GnomonicBase.reverse

897

898 @Property_RO

899 def geodesic(self):

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

901 '''

902 return self.datum.ellipsoid.geodsolve

903

904

905 class GnomonicKarney(_GnomonicBase):

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

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

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

909

910 @see: L{GnomonicExact}.

911 '''

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

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

914

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

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

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

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

919 radius (C{meter}).

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

921

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

923 not installed or not found.

924

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

926 '''

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

928

929 if _FOR_DOCS:

930 forward = _GnomonicBase.forward

931 reverse = _GnomonicBase.reverse

932

933 @Property_RO

934 def geodesic(self):

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

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

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

938 package is installed.

939 '''

940 return self.datum.ellipsoid.geodesic

941

942

943 class LambertEqualArea(_AzimuthalBase):

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

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

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

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

948 '''

949 if _FOR_DOCS:

950 __init__ = _AzimuthalBase.__init__

951

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

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

954

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

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

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

958

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

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

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

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

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

964

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

966 '''

967 def _k_t(c):

968 c += _1_0

969 t = c > EPS0

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

971 return k, t

972

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

974

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

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

977

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

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

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

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

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

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

984

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

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

987

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

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

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

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

992 perpendicular to this.

993 '''

994 def _c(c):

995 c *= _0_5

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

997

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

999

1000

1001 class Orthographic(_AzimuthalBase):

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

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

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

1005 '''

1006 if _FOR_DOCS:

1007 __init__ = _AzimuthalBase.__init__

1008

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

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

1011

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

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

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

1015

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

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

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

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

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

1021

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

1023 '''

1024 def _k_t(c):

1025 return _1_0, (c >= 0)

1026

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

1028

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

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

1031

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

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

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

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

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

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

1038

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

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

1041

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

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

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

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

1046 perpendicular to this.

1047 '''

1048 def _c(c):

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

1050

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

1052

1053

1054 class Stereographic(_AzimuthalBase):

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

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

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

1058 '''

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

1060 _k02 = _2_0 # double ._k0

1061

1062 if _FOR_DOCS:

1063 __init__ = _AzimuthalBase.__init__

1064

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

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

1067

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

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

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

1071

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

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

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

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

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

1077

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

1079 '''

1080 def _k_t(c):

1081 c += _1_0

1082 t = isnon0(c)

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

1084 return k, t

1085

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

1087

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

1089 def k0(self):

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

1091 '''

1092 return self._k0

1093

1094 @k0.setter # PYCHOK setter!

1095 def k0(self, factor):

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

1097 '''

1098 n = Stereographic.k0.fget.__name__

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

1100 self._k02 = self._k0 * _2_0

1101

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

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

1104

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

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

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

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

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

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

1111

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

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

1114

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

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

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

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

1119 perpendicular to this.

1120 '''

1121 def _c(c):

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

1123

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

1125

1126

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

1128

1129# **) MIT License

1130#

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

1132#

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

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

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

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

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

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

1139#

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

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

1142#

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

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

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

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

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

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

1149# OTHER DEALINGS IN THE SOFTWARE.