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, hypot as _hypot, Fsum

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

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, _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__ = '24.04.07'

75

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

77_over_horizon_ = 'over horizon'

78_TRIPS = 21 # numit, 4 sufficient

79

80

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

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

83 '''

84 e = Easting( x=x)

85 n = Northing(y=y)

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

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

88

89

90 class _AzimuthalBase(_NamedBase):

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

92

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

94 C++/doc/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 self._notOverloaded(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 self._notOverloaded(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 named

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 self._notOverloaded()

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).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, M=r.M12, name=name) if LatLon is None else \

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

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

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

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

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

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 name 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, name=NN, **LatLon_and_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 name: Optional name for the location (C{str}).

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

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

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

986

987 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{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, name, _c, True, **LatLon_and_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=NN):

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 name 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, name=NN, **LatLon_and_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 name: Optional name for the location (C{str}).

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

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

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

1040

1041 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{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, name, _c, False, **LatLon_and_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=NN):

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 name 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__

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, name=NN, **LatLon_and_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 name: Optional name for the location (C{str}).

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

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

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

1113

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

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

1116

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

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

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

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

1121 perpendicular to this.

1122 '''

1123 def _c(c):

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

1125

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

1127

1128

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

1130

1131# **) MIT License

1132#

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

1134#

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

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

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

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

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

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

1141#

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

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

1144#

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

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

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

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

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

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

1151# OTHER DEALINGS IN THE SOFTWARE.