1

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

3

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

5

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

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

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

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

10 and L{gnomonic}.

11

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

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

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

15 classGeographicLib_1_1GeodesicExact.html>}, respectively U{GeodesicLineExact

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

17

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

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

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

21 L{geodsolve} for more details.

22

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

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

25

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

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

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

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

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

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

32

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

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

35

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

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

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

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

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

41'''

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

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

44

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

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

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

48 from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB, \

49 _xinstanceof

50 from pygeodesy.datums import _spherical_datum, _WGS84

51 from pygeodesy.errors import _ValueError, _xdatum, _xkwds

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

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

54# from pygeodesy.formy import antipode # _MODS

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

56 _SPACE_, _x_, _y_

57 from pygeodesy.karney import _norm180

58 from pygeodesy.latlonBase import _MODS, LatLonBase as _LLB

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

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

61 from pygeodesy.namedTuples import LatLon2Tuple, LatLon4Tuple

62 from pygeodesy.props import deprecated_Property_RO, Property_RO, \

63 property_doc_, _update_all

64 from pygeodesy.streprs import Fmt, _fstrLL0, unstr

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

66 Scalar, Scalar_

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

68 sincos2d, sincos2d_

69

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

71

72__all__ = _ALL_LAZY.azimuthal

73__version__ = '24.06.11'

74

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

76_over_horizon_ = 'over horizon'

77_TRIPS = 21 # numit, 4 sufficient

78

79

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

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

82 '''

83 e = Easting( x=x)

84 n = Northing(y=y)

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

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

87

88

89 class _AzimuthalBase(_NamedBase):

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

91

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

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

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

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

96 '''

97 _datum = _WGS84 # L{Datum}

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

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

100

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

102 '''New azimuthal projection.

103

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

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

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

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

108 radius (C{meter}).

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

110

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

112

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

114 '''

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

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

117 if name:

118 self.name = name

119

120 if lat0 or lon0: # often both 0

121 self._reset(lat0, lon0)

122

123 @Property_RO

124 def datum(self):

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

126 '''

127 return self._datum

128

129 @Property_RO

130 def equatoradius(self):

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

132 '''

133 return self.datum.ellipsoid.a

134

135 a = equatoradius

136

137 @Property_RO

138 def flattening(self):

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

140 '''

141 return self.datum.ellipsoid.f

142

143 f = flattening

144

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

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

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

148

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

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

151 '''

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

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

154 s0, c0 = self._sc0

155

156 cb *= ca

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

158 if t:

159 r = k * self.radius

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

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

162 else: # 0 or 180

163 e = n = z = _0_0

164

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

166 name=self._name__(name))

167 return t

168

169 def _forwards(self, *lls):

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

171 '''

172 _fwd = self.forward

173 for ll in lls:

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

175

176 @Property_RO

177 def lat0(self):

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

179 '''

180 return self._latlon0.lat

181

182 @property

183 def latlon0(self):

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

185 '''

186 return self._latlon0

187

188 @latlon0.setter # PYCHOK setter!

189 def latlon0(self, latlon0):

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

191

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

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

194 '''

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

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

197 if hasattr(latlon0, _datum_):

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

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

200

201 @Property_RO

202 def lon0(self):

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

204 '''

205 return self._latlon0.lon

206

207 @deprecated_Property_RO

208 def majoradius(self): # PYCHOK no cover

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

210 return self.equatoradius

211

212 @Property_RO

213 def radius(self):

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

215 '''

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

217

218 def reset(self, lat0, lon0):

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

220

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

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

223

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

225 '''

226 _update_all(self) # zap caches

227 self._reset(lat0, lon0)

228

229 def _reset(self, lat0, lon0):

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

231 '''

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

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

234 self._sc0 = sincos2d(self.lat0)

235

236 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds): # PYCHOK no cover

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

238 self._notOverloaded(x, y, LatLon=LatLon, **name_LatLon_kwds)

239

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

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

242 '''

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

244

245 c = _c(r / self.radius)

246 if c is None:

247 lat, lon = self.latlon0

248 k, z = _1_0, _0_0

249 else:

250 s0, c0 = self._sc0

251 sc, cc = sincos2(c)

252 k = c / sc

253 s = s0 * cc

254 if r > EPS0:

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

256 lat = degrees(asin1(s))

257 if lea or fabs(c0) > EPS:

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

259 else:

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

261 lon = _norm180(self.lon0 + d)

262

263 if LatLon is None:

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

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

266 else:

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

268 return t

269

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

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

272 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne.

273 '''

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

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

276 return t, d

277

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

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

280 '''

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

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

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

284 _xinstanceof(B, LatLon=r)

285 return r

286

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

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

289

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

291

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

293 (C{str}).

294 '''

295 return _fstrLL0(self, prec, True)

296

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

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

299

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

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

302

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

304 '''

305 t = _fstrLL0(self, prec, False)

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

307

308

309 class AzimuthalError(_ValueError):

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

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

312 L{Azimuthal7Tuple} issue.

313 '''

314 pass

315

316

317 class Azimuthal7Tuple(_NamedTuple):

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

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

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

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

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

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

324 '''

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

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

327

328 def antipodal(self, azimuth=None):

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

330

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

332 (C{compass degrees360}).

333 '''

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

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

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

337

338

339 class Equidistant(_AzimuthalBase):

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

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

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

343

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

345 L{EquidistantGeodSolve} or L{EquidistantKarney} projection

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

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

348 '''

349 if _FOR_DOCS:

350 __init__ = _AzimuthalBase.__init__

351

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

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

354

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

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

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

358

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

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

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

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

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

364

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

366

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

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

369 '''

370 def _k_t(c):

371 k = _N_1_0 if c < 0 else _1_0

372 t = fabs(c) < EPS1

373 if t:

374 a = acos(c)

375 s = sin(a)

376 if s:

377 k = a / s

378 return k, t

379

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

381

382 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds):

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

384

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

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

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

388 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} forthe location and

389 optional, additional B{C{LatLon}} keyword arguments,

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

391

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

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

394

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

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

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

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

399 perpendicular to this.

400 '''

401 def _c(c):

402 return c if c > EPS else None

403

404 return self._reverse(x, y, _c, False, LatLon, **name_LatLon_kwds)

405

406

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

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

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

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

411

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

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

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

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

416 radius (C{meter}).

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

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

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

420

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

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

423

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

425

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

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

428

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

430 '''

431

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

433 if E is Equidistant:

434 try:

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

436 except ImportError:

437 pass

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

439

440

441 class _AzimuthalGeodesic(_AzimuthalBase):

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

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

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

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

446 '''

447 _mask = 0

448

449 @Property_RO

450 def geodesic(self): # PYCHOK no cover

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

452 self._notOverloaded()

453

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

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

456 '''

457 s = M

458 if s is None: # reciprocal, azimuthal scale

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

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

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

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

471 L{EquidistantKarney} projection.

472 '''

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

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

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

509 optional, additional B{C{LatLon}} keyword arguments,

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

511

512 @return: The geodetic (C{LatLon}) or if C{B{LatLon} is 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_LatLon_kwds) if LatLon is None else \

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

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

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

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

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

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

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

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

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

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

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

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

682

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

696

697

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

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

727 except ImportError:

728 pass

729 return G(lat0, lon0, datum=datum, **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):

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

738 '''

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

740

741 g = self.geodesic

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

743

744 def forward(self, lat, lon, raiser=True, **name): # 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 raiser: Do or don't throw an error (C{bool}) if

751 the location lies over the horizon.

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

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, _name__(name), M=M)

781 t._iteraton = 0

782 return t

783

784 def reverse(self, x, y, LatLon=None, **name_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 LatLon: Class to use (C{LatLon}) or C{None}.

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

791 optional, additional B{C{LatLon}} keyword arguments,

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

793

794 @return: The geodetic (C{LatLon}) or if C{B{LatLon} is 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).fsum2f_

824 _abs = fabs

825 for i in range(1, _TRIPS):

826 r = _P(s, outmask=m)

827 if _abs(d) < a:

828 break

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

830 else: # PYCHOK no cover

831 self._iteration = _TRIPS

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

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

834

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

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

837 t._iteration = self._iteration = i

838 return t

839

840

841 class GnomonicExact(_GnomonicBase):

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

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

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

845

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

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

848 '''

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

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

851

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

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

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

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

856 radius (C{meter}).

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

858

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

860 '''

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

862

863 if _FOR_DOCS:

864 forward = _GnomonicBase.forward

865 reverse = _GnomonicBase.reverse

866

867 @Property_RO

868 def geodesic(self):

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

870 '''

871 return self.datum.ellipsoid.geodesicx

872

873

874 class GnomonicGeodSolve(_GnomonicBase):

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

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

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

878 intended I{for testing purposes only}.

879

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

881 '''

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

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

884

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

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

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

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

889 radius (C{meter}).

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

891

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

893 '''

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

895

896 if _FOR_DOCS:

897 forward = _GnomonicBase.forward

898 reverse = _GnomonicBase.reverse

899

900 @Property_RO

901 def geodesic(self):

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

903 '''

904 return self.datum.ellipsoid.geodsolve

905

906

907 class GnomonicKarney(_GnomonicBase):

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

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

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

911

912 @see: L{GnomonicExact}.

913 '''

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

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

916

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

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

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

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

921 radius (C{meter}).

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

923

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

925 not installed or not found.

926

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

928 '''

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

930

931 if _FOR_DOCS:

932 forward = _GnomonicBase.forward

933 reverse = _GnomonicBase.reverse

934

935 @Property_RO

936 def geodesic(self):

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

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

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

940 package is installed.

941 '''

942 return self.datum.ellipsoid.geodesic

943

944

945 class LambertEqualArea(_AzimuthalBase):

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

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

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

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

950 '''

951 if _FOR_DOCS:

952 __init__ = _AzimuthalBase.__init__

953

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

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

956

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

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

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

960

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

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

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

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

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

966

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

968 '''

969 def _k_t(c):

970 c += _1_0

971 t = c > EPS0

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

973 return k, t

974

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

976

977 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds):

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

979

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

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

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

983 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} for the location

984 and optional, additional B{C{LatLon}} keyword

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

986

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

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

989

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

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

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

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

994 perpendicular to this.

995 '''

996 def _c(c):

997 c *= _0_5

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

999

1000 return self._reverse(x, y, _c, True, LatLon, **name_LatLon_kwds)

1001

1002

1003 class Orthographic(_AzimuthalBase):

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

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

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

1007 '''

1008 if _FOR_DOCS:

1009 __init__ = _AzimuthalBase.__init__

1010

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

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

1013

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

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

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

1017

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

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

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

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

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

1023

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

1025 '''

1026 def _k_t(c):

1027 return _1_0, (c >= 0)

1028

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

1030

1031 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds):

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

1033

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

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

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

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

1038 optional, additional B{C{LatLon}} keyword arguments,

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

1040

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

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

1043

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

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

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

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

1048 perpendicular to this.

1049 '''

1050 def _c(c):

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

1052

1053 return self._reverse(x, y, _c, False, LatLon, **name_LatLon_kwds)

1054

1055

1056 class Stereographic(_AzimuthalBase):

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

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

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

1060 '''

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

1062 _k02 = _2_0 # double ._k0

1063

1064 if _FOR_DOCS:

1065 __init__ = _AzimuthalBase.__init__

1066

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

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

1069

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

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

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

1073

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

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

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

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

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

1079

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

1081 '''

1082 def _k_t(c):

1083 c += _1_0

1084 t = isnon0(c)

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

1086 return k, t

1087

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

1089

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

1091 def k0(self):

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

1093 '''

1094 return self._k0

1095

1096 @k0.setter # PYCHOK setter!

1097 def k0(self, factor):

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

1099 '''

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

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

1102 self._k02 = self._k0 * _2_0

1103

1104 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds):

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

1106

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

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

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

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

1111 optional, additional B{C{LatLon}} keyword arguments,

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

1113

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

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

1116

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

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

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

1120 is C{1 / reciprocal} in the direction 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, _c, False, LatLon, **name_LatLon_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.