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

75

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

77_over_horizon_ = 'over horizon'

78_TRIPS = 21 # numit, 4 sufficient

79

80

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

137

138 @Property_RO

139 def flattening(self):

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

141 '''

142 return self.datum.ellipsoid.f

143

144 f = flattening

145

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

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

148 '''

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

150

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

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

153 '''

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

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

156 s0, c0 = self._sc0

157

158 cb *= ca

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

160 if t:

161 r = k * self.radius

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

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

164 else: # 0 or 180

165 e = n = z = _0_0

166

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

168 name=name or self.name)

169 return t

170

171 def _forwards(self, *lls):

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

173 '''

174 _fwd = self.forward

175 for ll in lls:

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

177

178 @Property_RO

179 def lat0(self):

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

181 '''

182 return self._latlon0.lat

183

184 @property

185 def latlon0(self):

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

187 '''

188 return self._latlon0

189

190 @latlon0.setter # PYCHOK setter!

191 def latlon0(self, latlon0):

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

193

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

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

196 '''

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

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

199 if hasattr(latlon0, _datum_):

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

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

202

203 @Property_RO

204 def lon0(self):

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

206 '''

207 return self._latlon0.lon

208

209 @deprecated_Property_RO

210 def majoradius(self): # PYCHOK no cover

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

212 return self.equatoradius

213

214 @Property_RO

215 def radius(self):

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

217 '''

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

219

220 def reset(self, lat0, lon0):

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

222

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

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

225

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

227 '''

228 _update_all(self) # zap caches

229 self._reset(lat0, lon0)

230

231 def _reset(self, lat0, lon0):

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

233 '''

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

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

236 self._sc0 = sincos2d(self.lat0)

237

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

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

240 '''

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

242

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

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

245 '''

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

247

248 c = _c(r / self.radius)

249 if c is None:

250 lat, lon = self.latlon0

251 k, z = _1_0, _0_0

252 else:

253 s0, c0 = self._sc0

254 sc, cc = sincos2(c)

255 k = c / sc

256 s = s0 * cc

257 if r > EPS0:

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

259 lat = degrees(asin1(s))

260 if lea or fabs(c0) > EPS:

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

262 else:

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

264 lon = _norm180(self.lon0 + d)

265

266 if LatLon is None:

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

268 name=name or self.name)

269 else:

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

271 return t

272

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

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

275 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne.

276 '''

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

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

279 return t, d

280

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

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

283 '''

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

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

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

287 _xinstanceof(B, LatLon=r)

288 return r

289

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

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

292

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

294

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

296 (C{str}).

297 '''

298 return _fstrLL0(self, prec, True)

299

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

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

302

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

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

305

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

307 '''

308 t = _fstrLL0(self, prec, False)

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

310

311

312 class AzimuthalError(_ValueError):

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

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

315 L{Azimuthal7Tuple} issue.

316 '''

317 pass

318

319

320 class Azimuthal7Tuple(_NamedTuple):

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

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

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

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

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

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

327 '''

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

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

330

331 def antipodal(self, azimuth=None):

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

333

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

335 (C{compass degrees360}).

336 '''

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

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

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

340

341

342 class Equidistant(_AzimuthalBase):

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

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

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

346

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

348 L{EquidistantGeodSolve} or L{EquidistantKarney} projection

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

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

351 '''

352 if _FOR_DOCS:

353 __init__ = _AzimuthalBase.__init__

354

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

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

357

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

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

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

361

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

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

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

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

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

367

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

369

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

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

372 '''

373 def _k_t(c):

374 k = _N_1_0 if c < 0 else _1_0

375 t = fabs(c) < EPS1

376 if t:

377 a = acos(c)

378 s = sin(a)

379 if s:

380 k = a / s

381 return k, t

382

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

384

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

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

387

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

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

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

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

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

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

394

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

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

397

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

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

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

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

402 perpendicular to this.

403 '''

404 def _c(c):

405 return c if c > EPS else None

406

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

408

409

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

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

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

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

414

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

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

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

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

419 radius (C{meter}).

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

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

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

423

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

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

426

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

428

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

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

431

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

433 '''

434

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

436 if E is Equidistant:

437 try:

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

439 except ImportError:

440 pass

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

442

443

444 class _AzimuthalGeodesic(_AzimuthalBase):

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

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

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

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

449 '''

450 _mask = 0

451

452 @Property_RO

453 def geodesic(self): # PYCHOK no cover

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

455 '''

456 notOverloaded(self)

457

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

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

460 '''

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

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

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

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

465 name=name or self.name)

466

467

468 class _EquidistantBase(_AzimuthalGeodesic):

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

470 and L{EquidistantKarney}.

471 '''

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

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

474 L{EquidistantKarney} projection.

475 '''

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

477

478 g = self.geodesic

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

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

481

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

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

484

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

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

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

488

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

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

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

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

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

494

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

496

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

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

499 '''

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

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

502 x, y = sincos2d(r.azi1)

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

504

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

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

507

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

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

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

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

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

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

514

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

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

517

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

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

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

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

522 to this.

523 '''

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

525

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

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

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

529

530

531 class EquidistantExact(_EquidistantBase):

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

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

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

535

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

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

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

539 clockwise from true North.

540

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

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

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

544 geographic and projected coordinates.

545 '''

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

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

548

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

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

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

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

553 radius (C{meter}).

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

555

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

557 '''

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

559

560 if _FOR_DOCS:

561 forward = _EquidistantBase.forward

562 reverse = _EquidistantBase.reverse

563

564 @Property_RO

565 def geodesic(self):

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

567 '''

568 return self.datum.ellipsoid.geodesicx

569

570

571 class EquidistantGeodSolve(_EquidistantBase):

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

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

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

575 I{for testing purposes only}.

576

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

578 '''

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

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

581

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

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

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

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

586 radius (C{meter}).

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

588

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

590 '''

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

592

593 if _FOR_DOCS:

594 forward = _EquidistantBase.forward

595 reverse = _EquidistantBase.reverse

596

597 @Property_RO

598 def geodesic(self):

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

600 '''

601 return self.datum.ellipsoid.geodsolve

602

603

604 class EquidistantKarney(_EquidistantBase):

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

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

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

608

609 @see: L{EquidistantExact}.

610 '''

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

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

613

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

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

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

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

618 radius (C{meter}).

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

620

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

622 not installed or not found.

623

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

625 '''

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

627

628 if _FOR_DOCS:

629 forward = _EquidistantBase.forward

630 reverse = _EquidistantBase.reverse

631

632 @Property_RO

633 def geodesic(self):

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

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

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

637 package is installed.

638 '''

639 return self.datum.ellipsoid.geodesic

640

641

642_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve,

643 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI

644

645

646 class Gnomonic(_AzimuthalBase):

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

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

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

650 '''

651 if _FOR_DOCS:

652 __init__ = _AzimuthalBase.__init__

653

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

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

656

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

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

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

660

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

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

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

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

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

666

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

668 '''

669 def _k_t(c):

670 t = c > EPS

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

672 return k, t

673

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

675

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

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

678

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

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

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

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

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

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

685

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

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

688

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

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

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

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

693 perpendicular to this.

694 '''

695 def _c(c):

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

697

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

699

700

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

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

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

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

705

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

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

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

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

710 radius (C{meter}).

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

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

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

714

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

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

717

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

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

720

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

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

723

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

725 '''

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

727 if G is Gnomonic:

728 try:

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

730 except ImportError:

731 pass

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

733

734

735 class _GnomonicBase(_AzimuthalGeodesic):

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

737 and L{GnomonicKarney}.

738 '''

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

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

741 '''

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

743

744 g = self.geodesic

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

746

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

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

749 and northing.

750

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

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

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

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

755 the location lies over the horizon.

756

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

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

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

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

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

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

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

764

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

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

767 '''

768 self._iteration = 0

769

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

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

772 M = r.M21

773 if M > EPS0:

774 q = r.m12 / M # .M12

775 e, n = sincos2d(r.azi1)

776 e *= q

777 n *= q

778 elif raiser: # PYCHOK no cover

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

780 else: # PYCHOK no cover

781 e = n = NAN

782

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

784 t._iteraton = 0

785 return t

786

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

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

789

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

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

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

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

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

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

796

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

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

799

800 @raise AzimuthalError: No convergence.

801

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

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

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

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

806 to this.

807 '''

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

809

810 d = a = self.equatoradius

811 s = a * atan(q / a)

812 if q > a: # PYCHOK no cover

813 def _d(r, q):

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

815

816 q = _1_0 / q

817 else: # little == True

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

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

820

821 a *= _EPS_K

822 m = self._mask

823 g = self.geodesic

824

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

826 _S2 = Fsum(s).fsum2_

827 for i in range(1, _TRIPS):

828 r = _P(s, outmask=m)

829 if fabs(d) < a:

830 break

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

832 else: # PYCHOK no cover

833 self._iteration = _TRIPS

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

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

836

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

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

839 t._iteration = self._iteration = i

840 return t

841

842

843 class GnomonicExact(_GnomonicBase):

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

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

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

847

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

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

850 '''

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

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

853

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

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

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

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

858 radius (C{meter}).

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

860

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

862 '''

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

864

865 if _FOR_DOCS:

866 forward = _GnomonicBase.forward

867 reverse = _GnomonicBase.reverse

868

869 @Property_RO

870 def geodesic(self):

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

872 '''

873 return self.datum.ellipsoid.geodesicx

874

875

876 class GnomonicGeodSolve(_GnomonicBase):

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

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

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

880 intended I{for testing purposes only}.

881

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

883 '''

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

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

886

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

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

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

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

891 radius (C{meter}).

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

893

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

895 '''

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

897

898 if _FOR_DOCS:

899 forward = _GnomonicBase.forward

900 reverse = _GnomonicBase.reverse

901

902 @Property_RO

903 def geodesic(self):

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

905 '''

906 return self.datum.ellipsoid.geodsolve

907

908

909 class GnomonicKarney(_GnomonicBase):

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

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

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

913

914 @see: L{GnomonicExact}.

915 '''

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

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

918

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

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

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

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

923 radius (C{meter}).

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

925

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

927 not installed or not found.

928

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

930 '''

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

932

933 if _FOR_DOCS:

934 forward = _GnomonicBase.forward

935 reverse = _GnomonicBase.reverse

936

937 @Property_RO

938 def geodesic(self):

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

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

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

942 package is installed.

943 '''

944 return self.datum.ellipsoid.geodesic

945

946

947 class LambertEqualArea(_AzimuthalBase):

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

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

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

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

952 '''

953 if _FOR_DOCS:

954 __init__ = _AzimuthalBase.__init__

955

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

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

958

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

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

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

962

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

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

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

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

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

968

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

970 '''

971 def _k_t(c):

972 c += _1_0

973 t = c > EPS0

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

975 return k, t

976

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

978

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

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

981

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

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

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

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

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

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

988

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

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

991

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

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

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

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

996 perpendicular to this.

997 '''

998 def _c(c):

999 c *= _0_5

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

1001

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

1003

1004

1005 class Orthographic(_AzimuthalBase):

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

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

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

1009 '''

1010 if _FOR_DOCS:

1011 __init__ = _AzimuthalBase.__init__

1012

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

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

1015

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

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

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

1019

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

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

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

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

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

1025

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

1027 '''

1028 def _k_t(c):

1029 return _1_0, (c >= 0)

1030

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

1032

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

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

1035

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

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

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

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

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

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

1042

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

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

1045

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

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

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

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

1050 perpendicular to this.

1051 '''

1052 def _c(c):

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

1054

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

1056

1057

1058 class Stereographic(_AzimuthalBase):

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

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

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

1062 '''

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

1064 _k02 = _2_0 # double ._k0

1065

1066 if _FOR_DOCS:

1067 __init__ = _AzimuthalBase.__init__

1068

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

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

1071

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

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

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

1075

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

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

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

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

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

1081

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

1083 '''

1084 def _k_t(c):

1085 c += _1_0

1086 t = isnon0(c)

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

1088 return k, t

1089

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

1091

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

1093 def k0(self):

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

1095 '''

1096 return self._k0

1097

1098 @k0.setter # PYCHOK setter!

1099 def k0(self, factor):

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

1101 '''

1102 n = Stereographic.k0.fget.__name__

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

1104 self._k02 = self._k0 * _2_0

1105

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

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

1108

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

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

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

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

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

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

1115

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

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

1118

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

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

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

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

1123 perpendicular to this.

1124 '''

1125 def _c(c):

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

1127

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

1129

1130

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

1132

1133# **) MIT License

1134#

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

1136#

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

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

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

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

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

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

1143#

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

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

1146#

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

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

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

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

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

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

1153# OTHER DEALINGS IN THE SOFTWARE.