1

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

3

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

5

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

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

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

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

10 and L{gnomonic}.

11

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

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

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

15 classGeographicLib_1_1GeodesicExact.html>}, respectively U{GeodesicLineExact

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

17

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

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

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

21 L{geodsolve} for more details.

22

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

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

25

26 Other azimuthal classes implement only (***) U{Snyder's FORMULAS FOR THE SPHERE

27<https://Pubs.USGS.gov/pp/1395/report.pdf>} and use those for any datum,

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

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

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

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

32

33 Page numbers in C{Snyder} references apply to U{John P. Snyder, "Map Projections

34-- A Working Manual", 1987<https://Pubs.USGS.gov/pp/1395/report.pdf>}.

35

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

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

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

39 altitude in Earth radii<https://WikiPedia.org/wiki/Azimuthal_equidistant_projection

40#/media/File:Comparison_azimuthal_projections.svg>}.

41'''

42# make sure int/int division yields float quotient, see .basics

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

44

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

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

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 _xinstanceof

51 from pygeodesy.datums import _spherical_datum, _WGS84

52 from pygeodesy.errors import _ValueError, _xdatum, _xkwds

53# from pygeodesy.fmath import euclid, Fsum, hypot as _hypot # from .karney

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

55# from pygeodesy.formy import antipode # _MODS

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

57 _scale_, _SPACE_, _x_, _y_

58 from pygeodesy.karney import _norm180, euclid, Fsum, _hypot

59 from pygeodesy.latlonBase import _MODS, LatLonBase as _LLB

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

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_, Northing, \

67 Scalar, Scalar_

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

69 sincos2d, sincos2d_

70

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

72

73__all__ = _ALL_LAZY.azimuthal

74__version__ = '23.12.12'

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

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

149

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

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

152 '''

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

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

155 s0, c0 = self._sc0

156

157 cb *= ca

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

159 if t:

160 r = k * self.radius

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

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

163 else: # 0 or 180

164 e = n = z = _0_0

165

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

167 name=name or self.name)

168 return t

169

170 def _forwards(self, *lls):

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

172 '''

173 _fwd = self.forward

174 for ll in lls:

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

176

177 @Property_RO

178 def lat0(self):

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

180 '''

181 return self._latlon0.lat

182

183 @property

184 def latlon0(self):

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

186 '''

187 return self._latlon0

188

189 @latlon0.setter # PYCHOK setter!

190 def latlon0(self, latlon0):

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

192

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

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

195 '''

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

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

198 if hasattr(latlon0, _datum_):

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

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

201

202 @Property_RO

203 def lon0(self):

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

205 '''

206 return self._latlon0.lon

207

208 @deprecated_Property_RO

209 def majoradius(self): # PYCHOK no cover

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

211 return self.equatoradius

212

213 @Property_RO

214 def radius(self):

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

216 '''

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

218

219 def reset(self, lat0, lon0):

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

221

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

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

224

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

226 '''

227 _update_all(self) # zap caches

228 self._reset(lat0, lon0)

229

230 def _reset(self, lat0, lon0):

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

232 '''

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

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

235 self._sc0 = sincos2d(self.lat0)

236

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

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

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

240

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

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

243 '''

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

245

246 c = _c(r / self.radius)

247 if c is None:

248 lat, lon = self.latlon0

249 k, z = _1_0, _0_0

250 else:

251 s0, c0 = self._sc0

252 sc, cc = sincos2(c)

253 k = c / sc

254 s = s0 * cc

255 if r > EPS0:

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

257 lat = degrees(asin1(s))

258 if lea or fabs(c0) > EPS:

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

260 else:

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

262 lon = _norm180(self.lon0 + d)

263

264 if LatLon is None:

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

266 name=name or self.name)

267 else:

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

269 return t

270

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

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

273 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne.

274 '''

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

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

277 return t, d

278

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

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

281 '''

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

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

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

285 _xinstanceof(B, LatLon=r)

286 return r

287

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

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

290

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

292

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

294 (C{str}).

295 '''

296 return _fstrLL0(self, prec, True)

297

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

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

300

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

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

303

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

305 '''

306 t = _fstrLL0(self, prec, False)

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

308

309

310 class AzimuthalError(_ValueError):

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

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

313 L{Azimuthal7Tuple} issue.

314 '''

315 pass

316

317

318 class Azimuthal7Tuple(_NamedTuple):

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

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

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

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

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

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

325 '''

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

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

328

329 def antipodal(self, azimuth=None):

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

331

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

333 (C{compass degrees360}).

334 '''

335 a = _MODS.formy.antipode(self.lat, self.lon) # PYCHOK _namedTuple

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

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

338

339

340 class Equidistant(_AzimuthalBase):

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

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

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

344

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

346 L{EquidistantGeodSolve} or L{EquidistantKarney} projection

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

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

349 '''

350 if _FOR_DOCS:

351 __init__ = _AzimuthalBase.__init__

352

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

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

355

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

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

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

359

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

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

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

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

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

365

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

367

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

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

370 '''

371 def _k_t(c):

372 k = _N_1_0 if c < 0 else _1_0

373 t = fabs(c) < EPS1

374 if t:

375 a = acos(c)

376 s = sin(a)

377 if s:

378 k = a / s

379 return k, t

380

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

382

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

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

385

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

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

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

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

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

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

392

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

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

395

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

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

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

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

400 perpendicular to this.

401 '''

402 def _c(c):

403 return c if c > EPS else None

404

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

406

407

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

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

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

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

412

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

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

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

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

417 radius (C{meter}).

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

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

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

421

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

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

424

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

426

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

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

429

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

431 '''

432

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

434 if E is Equidistant:

435 try:

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

437 except ImportError:

438 pass

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

440

441

442 class _AzimuthalGeodesic(_AzimuthalBase):

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

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

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

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

447 '''

448 _mask = 0

449

450 @Property_RO

451 def geodesic(self): # PYCHOK no cover

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

453 notOverloaded(self)

454

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

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

457 '''

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

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

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

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

462 name=name or self.name)

463

464

465 class _EquidistantBase(_AzimuthalGeodesic):

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

467 and L{EquidistantKarney}.

468 '''

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

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

471 L{EquidistantKarney} projection.

472 '''

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

474

475 g = self.geodesic

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

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

478

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

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

481

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

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

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

485

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

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

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

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

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

491

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

493

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

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

496 '''

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

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

499 x, y = sincos2d(r.azi1)

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

501

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

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

504

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

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

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

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

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

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

511

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

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

514

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

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

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

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

519 to this.

520 '''

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

522

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

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

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

526

527

528 class EquidistantExact(_EquidistantBase):

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

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

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

532

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

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

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

536 clockwise from true North.

537

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

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

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

541 geographic and projected coordinates.

542 '''

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

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

545

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

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

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

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

550 radius (C{meter}).

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

552

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

554 '''

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

556

557 if _FOR_DOCS:

558 forward = _EquidistantBase.forward

559 reverse = _EquidistantBase.reverse

560

561 @Property_RO

562 def geodesic(self):

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

564 '''

565 return self.datum.ellipsoid.geodesicx

566

567

568 class EquidistantGeodSolve(_EquidistantBase):

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

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

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

572 I{for testing purposes only}.

573

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

575 '''

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

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

578

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

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

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

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

583 radius (C{meter}).

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

585

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

587 '''

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

589

590 if _FOR_DOCS:

591 forward = _EquidistantBase.forward

592 reverse = _EquidistantBase.reverse

593

594 @Property_RO

595 def geodesic(self):

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

597 '''

598 return self.datum.ellipsoid.geodsolve

599

600

601 class EquidistantKarney(_EquidistantBase):

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

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

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

605

606 @see: L{EquidistantExact}.

607 '''

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

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

610

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

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

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

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

615 radius (C{meter}).

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

617

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

619 not installed or not found.

620

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

622 '''

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

624

625 if _FOR_DOCS:

626 forward = _EquidistantBase.forward

627 reverse = _EquidistantBase.reverse

628

629 @Property_RO

630 def geodesic(self):

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

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

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

634 package is installed.

635 '''

636 return self.datum.ellipsoid.geodesic

637

638

639_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve,

640 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI

641

642

643 class Gnomonic(_AzimuthalBase):

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

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

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

647 '''

648 if _FOR_DOCS:

649 __init__ = _AzimuthalBase.__init__

650

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

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

653

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

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

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

657

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

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

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

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

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

663

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

665 '''

666 def _k_t(c):

667 t = c > EPS

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

669 return k, t

670

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

672

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

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

675

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

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

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

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

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

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

682

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

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

685

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

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

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

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

690 perpendicular to this.

691 '''

692 def _c(c):

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

694

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

696

697

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

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

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

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

702

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

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

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

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

707 radius (C{meter}).

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

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

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

711

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

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

714

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

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

717

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

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

720

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

722 '''

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

724 if G is Gnomonic:

725 try:

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

727 except ImportError:

728 pass

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

730

731

732 class _GnomonicBase(_AzimuthalGeodesic):

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

734 and L{GnomonicKarney}.

735 '''

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

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

738 '''

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

740

741 g = self.geodesic

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

743

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

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

746 and northing.

747

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

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

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

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

752 the location lies over the horizon.

753

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

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

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

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

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

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

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

761

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

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

764 '''

765 self._iteration = 0

766

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

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

769 M = r.M21

770 if M > EPS0:

771 q = r.m12 / M # .M12

772 e, n = sincos2d(r.azi1)

773 e *= q

774 n *= q

775 elif raiser: # PYCHOK no cover

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

777 else: # PYCHOK no cover

778 e = n = NAN

779

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

781 t._iteraton = 0

782 return t

783

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

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

786

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

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

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

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

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

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

793

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

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

796

797 @raise AzimuthalError: No convergence.

798

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

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

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

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

803 to this.

804 '''

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

806

807 d = a = self.equatoradius

808 s = a * atan1(q, a)

809 if q > a: # PYCHOK no cover

810 def _d(r, q):

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

812

813 q = _1_0 / q

814 else: # little == True

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

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

817

818 a *= _EPS_K

819 m = self._mask

820 g = self.geodesic

821

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

823 _S2 = Fsum(s).fsum2_

824 for i in range(1, _TRIPS):

825 r = _P(s, outmask=m)

826 if fabs(d) < a:

827 break

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

829 else: # PYCHOK no cover

830 self._iteration = _TRIPS

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

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

833

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

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

836 t._iteration = self._iteration = i

837 return t

838

839

840 class GnomonicExact(_GnomonicBase):

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

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

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

844

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

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

847 '''

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

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

850

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

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

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

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

855 radius (C{meter}).

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

857

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

859 '''

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

861

862 if _FOR_DOCS:

863 forward = _GnomonicBase.forward

864 reverse = _GnomonicBase.reverse

865

866 @Property_RO

867 def geodesic(self):

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

869 '''

870 return self.datum.ellipsoid.geodesicx

871

872

873 class GnomonicGeodSolve(_GnomonicBase):

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

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

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

877 intended I{for testing purposes only}.

878

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

880 '''

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

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

883

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

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

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

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

888 radius (C{meter}).

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

890

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

892 '''

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

894

895 if _FOR_DOCS:

896 forward = _GnomonicBase.forward

897 reverse = _GnomonicBase.reverse

898

899 @Property_RO

900 def geodesic(self):

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

902 '''

903 return self.datum.ellipsoid.geodsolve

904

905

906 class GnomonicKarney(_GnomonicBase):

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

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

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

910

911 @see: L{GnomonicExact}.

912 '''

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

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

915

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

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

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

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

920 radius (C{meter}).

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

922

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

924 not installed or not found.

925

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

927 '''

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

929

930 if _FOR_DOCS:

931 forward = _GnomonicBase.forward

932 reverse = _GnomonicBase.reverse

933

934 @Property_RO

935 def geodesic(self):

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

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

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

939 package is installed.

940 '''

941 return self.datum.ellipsoid.geodesic

942

943

944 class LambertEqualArea(_AzimuthalBase):

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

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

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

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

949 '''

950 if _FOR_DOCS:

951 __init__ = _AzimuthalBase.__init__

952

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

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

955

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

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

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

959

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

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

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

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

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

965

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

967 '''

968 def _k_t(c):

969 c += _1_0

970 t = c > EPS0

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

972 return k, t

973

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

975

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

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

978

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

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

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

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

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

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

985

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

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

988

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

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

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

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

993 perpendicular to this.

994 '''

995 def _c(c):

996 c *= _0_5

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

998

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

1000

1001

1002 class Orthographic(_AzimuthalBase):

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

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

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

1006 '''

1007 if _FOR_DOCS:

1008 __init__ = _AzimuthalBase.__init__

1009

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

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

1012

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

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

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

1016

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

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

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

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

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

1022

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

1024 '''

1025 def _k_t(c):

1026 return _1_0, (c >= 0)

1027

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

1029

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

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

1032

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

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

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

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

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

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

1039

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

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

1042

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

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

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

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

1047 perpendicular to this.

1048 '''

1049 def _c(c):

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

1051

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

1053

1054

1055 class Stereographic(_AzimuthalBase):

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

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

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

1059 '''

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

1061 _k02 = _2_0 # double ._k0

1062

1063 if _FOR_DOCS:

1064 __init__ = _AzimuthalBase.__init__

1065

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

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

1068

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

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

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

1072

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

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

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

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

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

1078

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

1080 '''

1081 def _k_t(c):

1082 c += _1_0

1083 t = isnon0(c)

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

1085 return k, t

1086

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

1088

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

1090 def k0(self):

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

1092 '''

1093 return self._k0

1094

1095 @k0.setter # PYCHOK setter!

1096 def k0(self, factor):

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

1098 '''

1099 n = Stereographic.k0.fget.__name__

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

1101 self._k02 = self._k0 * _2_0

1102

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

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

1105

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

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

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

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

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

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

1112

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

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

1115

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

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

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

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

1120 perpendicular to this.

1121 '''

1122 def _c(c):

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

1124

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

1126

1127

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

1129

1130# **) MIT License

1131#

1132# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

1133#

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

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

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

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

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

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

1140#

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

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

1143#

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

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

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

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

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

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

1150# OTHER DEALINGS IN THE SOFTWARE.