Coverage for pygeodesy/albers.py: 97%

423 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-03-03 11:31 -0500

1 

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

3 

4u'''Albers Equal-Area projections. 

5 

6Classes L{AlbersEqualArea}, L{AlbersEqualArea2}, L{AlbersEqualArea4}, 

7L{AlbersEqualAreaCylindrical}, L{AlbersEqualAreaNorth}, L{AlbersEqualAreaSouth} 

8and L{AlbersError}, transcoded from I{Charles Karney}'s C++ class U{AlbersEqualArea 

9<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AlbersEqualArea.html>}. 

10 

11See also I{Albers Equal-Area Conic Projection} in U{John P. Snyder, "Map Projections 

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

13and the Albers Conical Equal-Area examples on pp 291-294. 

14''' 

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

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

17 

18from pygeodesy.basics import neg, neg_ 

19from pygeodesy.constants import EPS0, EPS02, _EPSqrt as _TOL, \ 

20 _0_0, _0_5, _1_0, _N_1_0, _2_0, \ 

21 _N_2_0, _4_0, _6_0, _90_0, _N_90_0 

22from pygeodesy.datums import _ellipsoidal_datum, _WGS84 

23from pygeodesy.errors import _ValueError, _xkwds 

24from pygeodesy.fmath import hypot, hypot1, sqrt3 

25from pygeodesy.fsums import Fsum, fsum1f_ 

26from pygeodesy.interns import NN, _COMMASPACE_, _datum_, _gamma_, _k0_, \ 

27 _lat_, _lat1_, _lat2_, _lon_, _name_, _not_, \ 

28 _negative_, _scale_, _SPACE_, _x_, _y_ 

29from pygeodesy.karney import _diff182, _norm180, _signBit 

30from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY 

31from pygeodesy.named import _NamedBase, _NamedTuple, _Pass 

32from pygeodesy.props import deprecated_Property_RO, Property_RO, _update_all 

33from pygeodesy.streprs import Fmt, strs, unstr 

34from pygeodesy.units import Bearing, Float_, Lat, Lat_, Lon, Meter, Scalar_ 

35from pygeodesy.utily import atan1, atan1d, degrees360, sincos2, sincos2d, \ 

36 sincos2d_ 

37 

38from math import atan2, atanh, degrees, fabs, radians, sqrt 

39 

40__all__ = _ALL_LAZY.albers 

41__version__ = '23.12.01' 

42 

43_k1_ = 'k1' 

44_NUMIT = 8 # XXX 4? 

45_NUMIT0 = 41 # XXX 21? 

46_TERMS = 21 # XXX 16? 

47_TOL0 = sqrt3(_TOL) 

48 

49 

50def _ct2(s, c): 

51 '''(INTERNAL) Avoid singularities at poles. 

52 ''' 

53 c = max(EPS0, c) 

54 return c, (s / c) 

55 

56 

57def _Fsum1(a, b, c=_0_0): # floats=True 

58 '''(INTERNAL) C{Fsum} 1-primed. 

59 ''' 

60 S = Fsum() 

61 S._facc_(_1_0, a, b, c, _N_1_0, up=False) 

62 return S 

63 

64 

65def _Ks(**name_k): 

66 '''(INTERNAL) Scale C{B{k} >= EPS0}. 

67 ''' 

68 return Scalar_(Error=AlbersError, low=EPS0, **name_k) # > 0 

69 

70 

71def _Lat(*lat, **Error_name_lat): 

72 '''(INTERNAL) Latitude C{-90 <= B{lat} <= 90}. 

73 ''' 

74 kwds = _xkwds(Error_name_lat, Error=AlbersError) 

75 return Lat_(*lat, **kwds) 

76 

77 

78def _qZx(albs): 

79 '''(INTERNAL) Set C{albs._qZ} and C{albs._qx}. 

80 ''' 

81 E = albs._datum.ellipsoid # _AlbersBase 

82 albs._qZ = qZ = _1_0 + E.e21 * _atanheE(_1_0, E) 

83 albs._qx = qZ / (_2_0 * E.e21) 

84 return qZ 

85 

86 

87class AlbersError(_ValueError): 

88 '''An L{AlbersEqualArea}, L{AlbersEqualArea2}, L{AlbersEqualArea4}, 

89 L{AlbersEqualAreaCylindrical}, L{AlbersEqualAreaNorth}, 

90 L{AlbersEqualAreaSouth} or L{Albers7Tuple} issue. 

91 ''' 

92 pass 

93 

94 

95class _AlbersBase(_NamedBase): 

96 '''(INTERNAL) Base class for C{AlbersEqualArea...} projections. 

97 

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

99 html/classGeographicLib_1_1AlbersEqualArea.html>}, method C{Init}. 

100 ''' 

101 _datum = _WGS84 

102 _k = NN # or _k0_ or _k1_ 

103 _k0 = _Ks(k0=_1_0) 

104# _k0n0 = None # (INTERNAL) k0 * no 

105 _k02 = _1_0 # (INTERNAL) k0**2 

106# _k02n0 = None # (INTERNAL) k02 * n0 

107# _lat0 = None # lat origin 

108 _lat1 = None # let 1st parallel 

109 _lat2 = None # lat 2nd parallel 

110 _m0 = _0_0 # if polar else sqrt(m02) 

111# _m02 = None # (INTERNAL) cached 

112# _n0 = None # (INTERNAL) cached 

113 _nrho0 = _0_0 # if polar else m0 * E.a 

114 _polar = False 

115 _qx = None # (INTERNAL) see _qZx 

116 _qZ = None # (INTERNAL) see _qZx 

117# _scxi0_ = None # (INTERNAL) sec(xi) / (qZ * E.a2) 

118 _sign = +1 

119# _sxi0 = None # (INTERNAL) sin(xi) 

120# _txi0 = None # (INTERNAL) tan(xi) 

121 

122 def __init__(self, sa1, ca1, sa2, ca2, k, datum, name): 

123 '''(INTERNAL) New C{AlbersEqualArea...} instance. 

124 ''' 

125 qZ = self._qZ 

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

127 self._datum = _ellipsoidal_datum(datum, name=name) 

128 qZ = _qZx(self) 

129 elif qZ is None: 

130 qZ = _qZx(_AlbersBase) 

131 if name: 

132 self.name = name 

133 

134 E = self.ellipsoid 

135 c = min(ca1, ca2) 

136 if _signBit(c): 

137 raise AlbersError(clat1=ca1, clat2=ca2, txt=_negative_) 

138 polar = c < EPS0 # == 0 

139 

140 # determine hemisphere of tangent latitude 

141 if sa1 < 0: # and sa2 < 0: 

142 self._sign = -1 

143 # internally, tangent latitude positive 

144 sa1, sa2 = neg_(sa1, sa2) 

145 if sa1 > sa2: # make phi1 < phi2 

146 sa1, sa2 = sa2, sa1 

147 ca1, ca2 = ca2, ca1 

148 if sa1 < 0: # or sa2 < 0: 

149 raise AlbersError(slat1=sa1, slat2=sa2, txt=_negative_) 

150 ca1, ta1 = _ct2(sa1, ca1) 

151 ca2, ta2 = _ct2(sa2, ca2) 

152 

153 par1 = fabs(ta1 - ta2) < EPS02 # ta1 == ta2 

154 if par1 or polar: 

155 ta0, C = ta2, _1_0 

156 else: 

157 ta0, C = self._ta0C2(ca1, sa1, ta1, ca2, sa2, ta2) 

158 

159 self._lat0 = _Lat(lat0=self._sign * atan1d(ta0)) 

160 self._m02 = m02 = _1_x21(E.f1 * ta0) 

161 self._n0 = n0 = ta0 / hypot1(ta0) 

162 if polar: 

163 self._polar = True 

164# self._nrho0 = self._m0 = _0_0 

165 else: # m0 = nrho0 / E.a 

166 self._m0 = sqrt(m02) 

167 self._nrho0 = self._m0 * E.a 

168 t = self._txi0 = self._txif(ta0) 

169 h = hypot1(t) 

170 s = self._sxi0 = t / h 

171 if par1: 

172 self._k0n0 = self._k02n0 = n0 

173 else: 

174 self._k0s(k * sqrt(C / (m02 + n0 * qZ * s))) 

175 self._scxi0_ = h / (qZ * E.a2) 

176 

177 def _a_b_sxi3(self, *ca_sa_ta_scb_4s): 

178 '''(INTERNAL) Sum of C{sm1} terms and C{sin(xi)}s for ._ta0C2. 

179 ''' 

180 _1 = _1_0 

181 a = b = s = _0_0 

182 for ca, sa, ta, scb in ca_sa_ta_scb_4s: 

183 cxi, sxi, _ = self._cstxif3(ta) 

184 if sa > 0: 

185 sa += _1 

186 a += (cxi / ca)**2 * sa / (sxi + _1) 

187 b += scb * ca**2 / sa 

188 else: 

189 sa = _1 - sa 

190 a += (_1 - sxi) / sa 

191 b += scb * sa 

192 s += sxi 

193 return a, b, s 

194 

195 def _azik(self, t, ta): 

196 '''(INTERNAL) Compute the azimuthal scale C{_Ks(k=k)}. 

197 ''' 

198 E = self.ellipsoid 

199 return _Ks(k=hypot1(E.b_a * ta) * self._k0 * t / E.a) 

200 

201 def _cstxif3(self, ta): 

202 '''(INTERNAL) Get 3-tuple C{(cos, sin, tan)} of M{xi(ta)}. 

203 ''' 

204 t = self._txif(ta) 

205 c = _1_0 / hypot1(t) 

206 s = c * t 

207 return c, s, t 

208 

209 @Property_RO 

210 def datum(self): 

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

212 ''' 

213 return self._datum 

214 

215 @Property_RO 

216 def ellipsoid(self): 

217 '''Get the datum's ellipsoid (L{Ellipsoid}). 

218 ''' 

219 return self.datum.ellipsoid 

220 

221 @Property_RO 

222 def equatoradius(self): 

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

224 ''' 

225 return self.ellipsoid.a 

226 

227 a = equatoradius 

228 

229 @Property_RO 

230 def flattening(self): 

231 '''Get the C{ellipsoid}'s flattening (C{scalar}). 

232 ''' 

233 return self.ellipsoid.f 

234 

235 f = flattening 

236 

237 def forward(self, lat, lon, lon0=0, name=NN): 

238 '''Convert a geodetic location to east- and northing. 

239 

240 @arg lat: Latitude of the location (C{degrees}). 

241 @arg lon: Longitude of the location (C{degrees}). 

242 @kwarg lon0: Optional central meridian longitude (C{degrees}). 

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

244 

245 @return: An L{Albers7Tuple}C{(x, y, lat, lon, gamma, scale, datum)}, 

246 with C{lon} offset by B{C{lon0}} and reduced C{[-180,180]}. 

247 

248 @note: The origin latitude is returned by C{property lat0}. No 

249 false easting or northing is added. The value of B{C{lat}} 

250 should be in the range C{[-90..90] degrees}. The returned 

251 values C{x} and C{y} will be large but finite for points 

252 projecting to infinity, i.e. one or both of the poles. 

253 ''' 

254 a = self.ellipsoid.a 

255 s = self._sign 

256 

257 k0 = self._k0 

258 n0 = self._n0 

259 nrho0 = self._nrho0 

260 txi0 = self._txi0 

261 

262 _, ta = _ct2(*sincos2d(s * _Lat(lat=lat))) 

263 

264 _, sxi, txi = self._cstxif3(ta) 

265 dq = _Dsn(txi, txi0, sxi, self._sxi0) * \ 

266 (txi - txi0) * self._qZ 

267 drho = a * dq / (sqrt(self._m02 - n0 * dq) + self._m0) 

268 

269 lon, _ = _diff182(lon0, lon) 

270 x = radians(lon) 

271 th = self._k02n0 * x 

272 sth, cth = sincos2(th) # XXX sin, cos 

273 if n0: 

274 x = sth / n0 

275 y = (_1_0 - cth) if cth < 0 else (sth**2 / (cth + _1_0)) 

276 y *= nrho0 / n0 

277 else: 

278 x *= self._k02 

279 y = _0_0 

280 t = nrho0 - n0 * drho 

281 x = t * x / k0 

282 y = s * (y + drho * cth) / k0 

283 

284 g = degrees360(s * th) 

285 if t: 

286 k0 = self._azik(t, ta) 

287 return Albers7Tuple(x, y, lat, lon, g, k0, self.datum, 

288 name=name or self.name) 

289 

290 @Property_RO 

291 def ispolar(self): 

292 '''Is this projection polar (C{bool})? 

293 ''' 

294 return self._polar 

295 

296 isPolar = ispolar # synonym 

297 

298 def _k0s(self, k0): 

299 '''(INTERNAL) Set C{._k0}, C{._k02}, etc. 

300 ''' 

301 self._k0 = k0 = _Ks(k0=k0) 

302 self._k02 = k02 = k0**2 

303 self._k0n0 = k0 * self._n0 

304 self._k02n0 = k02 * self._n0 

305 

306 @Property_RO 

307 def lat0(self): 

308 '''Get the latitude of the projection origin (C{degrees}). 

309 

310 This is the latitude of minimum azimuthal scale and 

311 equals the B{C{lat}} in the 1-parallel L{AlbersEqualArea} 

312 and lies between B{C{lat1}} and B{C{lat2}} for the 

313 2-parallel L{AlbersEqualArea2} and L{AlbersEqualArea4} 

314 projections. 

315 ''' 

316 return self._lat0 

317 

318 @Property_RO 

319 def lat1(self): 

320 '''Get the latitude of the first parallel (C{degrees}). 

321 ''' 

322 return self._lat1 

323 

324 @Property_RO 

325 def lat2(self): 

326 '''Get the latitude of the second parallel (C{degrees}). 

327 

328 @note: The second and first parallel latitudes are the 

329 same instance for 1-parallel C{AlbersEqualArea*} 

330 projections. 

331 ''' 

332 return self._lat2 

333 

334 @deprecated_Property_RO 

335 def majoradius(self): # PYCHOK no cover 

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

337 return self.equatoradius 

338 

339 def rescale0(self, lat, k=1): # PYCHOK no cover 

340 '''Set the azimuthal scale for this projection. 

341 

342 @arg lat: Northern latitude (C{degrees}). 

343 @arg k: Azimuthal scale at latitude B{C{lat}} (C{scalar}). 

344 

345 @raise AlbersError: Invalid B{C{lat}} or B{C{k}}. 

346 

347 @note: This allows a I{latitude of conformality} to be specified. 

348 ''' 

349 k0 = _Ks(k=k) / self.forward(lat, _0_0).scale 

350 if self._k0 != k0: 

351 _update_all(self) 

352 self._k0s(k0) 

353 

354 def reverse(self, x, y, lon0=0, name=NN, LatLon=None, **LatLon_kwds): 

355 '''Convert an east- and northing location to geodetic lat- and longitude. 

356 

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

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

359 @kwarg lon0: Optional central meridian longitude (C{degrees}). 

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

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

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

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

364 

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

366 L{Albers7Tuple}C{(x, y, lat, lon, gamma, scale, datum)}. 

367 

368 @note: The origin latitude is returned by C{property lat0}. No 

369 false easting or northing is added. The returned value of 

370 C{lon} is in the range C{[-180..180] degrees} and C{lat} 

371 is in the range C{[-90..90] degrees}. If the given 

372 B{C{x}} or B{C{y}} point is outside the valid projected 

373 space the nearest pole is returned. 

374 ''' 

375 k0 = self._k0 

376 n0 = self._n0 

377 k0n0 = self._k0n0 

378 s = self._sign 

379 txi = self._txi0 

380 

381 x = Meter(x=x) 

382 nx = k0n0 * x 

383 y = Meter(y=y) 

384 y_ = s * y 

385 ny = k0n0 * y_ 

386 t = nrho0 = self._nrho0 

387 y1 = nrho0 - ny 

388 

389 den = hypot(nx, y1) + nrho0 # 0 implies origin with polar aspect 

390 if den: 

391 drho = _Fsum1(x * nx, y_ * nrho0 * _N_2_0, y_ * ny).fover(den / k0) 

392 # dsxia = scxi0 * dsxi 

393 t += drho * n0 

394 d_ = (nrho0 + t) * drho * self._scxi0_ # / (qZ * E.a2) 

395 d_2 = (txi * _2_0 - d_) * d_ + _1_0 

396 txi = (txi - d_) / (sqrt(d_2) if d_2 > EPS02 else EPS0) 

397 

398 ta = self._tanf(txi) 

399 lat = atan1d(s * ta) 

400 

401 th = atan2(nx, y1) 

402 lon = degrees((th / self._k02n0) if n0 else (x / (y1 * k0))) 

403 if lon0: 

404 lon += _norm180(lon0) 

405 lon = _norm180(lon) 

406 

407 n = name or self.name 

408 if LatLon is None: 

409 g = degrees360(s * th) 

410 if den: 

411 k0 = self._azik(t, ta) 

412 r = Albers7Tuple(x, y, lat, lon, g, k0, self.datum, name=n) 

413 else: # PYCHOK no cover 

414 kwds = _xkwds(LatLon_kwds, datum=self.datum, name=n) 

415 r = LatLon(lat, lon, **kwds) 

416 return r 

417 

418 @Property_RO 

419 def scale0(self): 

420 '''Get the central scale for the projection (C{float}). 

421 

422 This is the azimuthal scale on the latitude of origin 

423 of the projection, see C{property lat0}. 

424 ''' 

425 return self._k0 

426 

427 def _ta0(self, s1_qZ, ta0, E): 

428 '''(INTERNAL) Refine C{ta0} for C{._ta0C2}. 

429 ''' 

430 e2 = E.e2 

431 e21 = E.e21 

432 e22 = E.e22 # == e2 / e21 

433 tol = _tol(_TOL0, ta0) 

434 _Ta02 = Fsum(ta0).fsum2_ 

435 _fabs = fabs 

436 _fsum1 = fsum1f_ 

437 _sqrt = sqrt 

438 _1, _2 = _1_0, _2_0 

439 _4, _6 = _4_0, _6_0 

440 for self._iteration in range(1, _NUMIT0): # 4 trips 

441 ta02 = ta0**2 

442 sca02 = ta02 + _1 

443 sca0 = _sqrt(sca02) 

444 sa0 = ta0 / sca0 

445 sa01 = sa0 + _1 

446 sa02 = sa0**2 

447 # sa0m = 1 - sa0 = 1 / (sec(a0) * (tan(a0) + sec(a0))) 

448 sa0m = _1 / (sca0 * (ta0 + sca0)) # scb0^2 * sa0 

449 sa0m1 = sa0m / (_1 - e2 * sa0) 

450 sa021 = _1 - e2 * sa02 

451 

452 g = (_1 + ta02 * e21) * sa0 

453 dg = (_1 + ta02 * _2) * sca02 * e21 + e2 

454 D = (_1 - (_1 + sa0 * _2 * sa01) * e2) * sa0m / (e21 * sa01) # dD/dsa0 

455 dD = (_2 - (_6 + sa0 * _4) * sa02 * e2) / (e21 * sa01**2) 

456 BA = (_atanh1(e2 * sa0m1**2) * e21 - e2 * sa0m) * sa0m1 \ 

457 - (_2 + (_1 + e2) * sa0) * sa0m**2 * e22 / sa021 # B + A 

458 d = (_4 - (_1 + sa02) * e2 * _2) * e22 / (sa021**2 * sca02) # dAB 

459 u = _fsum1(s1_qZ * g, -D, g * BA) 

460 du = _fsum1(s1_qZ * dg, dD, dg * BA, g * d) 

461 ta0, d = _Ta02(-u / du * (sca0 * sca02)) 

462 if _fabs(d) < tol: 

463 return ta0 

464 raise AlbersError(Fmt.no_convergence(d, tol), txt=repr(self)) 

465 

466 def _ta0C2(self, ca1, sa1, ta1, ca2, sa2, ta2): 

467 '''(INTERNAL) Compute C{ta0} and C{C} for C{.__init__}. 

468 ''' 

469 E = self.ellipsoid 

470 f1, e2 = E.f1, E.e2 

471 _1 = _1_0 

472 

473 tb1 = f1 * ta1 

474 tb2 = f1 * ta2 

475 dtb12 = f1 * (tb1 + tb2) 

476 scb12 = _1 + tb1**2 

477 scb22 = _1 + tb2**2 

478 

479 dsn_2 = _Dsn(ta2, ta1, sa2, sa1) * _0_5 

480 sa12 = sa1 * sa2 

481 

482 esa1_2 = (_1 - e2 * sa1**2) \ 

483 * (_1 - e2 * sa2**2) 

484 esa12 = _1 + e2 * sa12 

485 

486 axi, bxi, sxi = self._a_b_sxi3((ca1, sa1, ta1, scb12), 

487 (ca2, sa2, ta2, scb22)) 

488 

489 dsxi = ((esa12 / esa1_2) + _DatanheE(sa2, sa1, E)) * dsn_2 / self._qx 

490 C = _Fsum1(sxi * dtb12 / dsxi, scb22, scb12).fover(scb22 * scb12 * _2_0) 

491 

492 sa12 = fsum1f_(sa1, sa2, sa12) 

493 axi *= (sa12 * e2 + _1) / (sa12 + _1) 

494 bxi *= _Fsum1(sa1, sa2, esa12).fover(esa1_2) * e2 + _D2atanheE(sa1, sa2, E) * E.e21 

495 s1_qZ = (axi * self._qZ - bxi) * dsn_2 / dtb12 

496 ta0 = self._ta0(s1_qZ, (ta1 + ta2) * _0_5, E) 

497 return ta0, C 

498 

499 def _tanf(self, txi): # in .Ellipsoid.auxAuthalic 

500 '''(INTERNAL) Function M{tan-phi from tan-xi}. 

501 ''' 

502 tol = _tol(_TOL, txi) 

503 e2 = self.ellipsoid.e2 

504 qx = self._qx 

505 

506 ta = txi 

507 _Ta2 = Fsum(ta).fsum2_ 

508 _fabs = fabs 

509 _sqrt3 = sqrt3 

510 _txif = self._txif 

511 _1 = _1_0 

512 for self._iteration in range(1, _NUMIT): # max 2, mean 1.99 

513 # dtxi / dta = (scxi / sca)^3 * 2 * (1 - e^2) 

514 # / (qZ * (1 - e^2 * sa^2)^2) 

515 ta2 = ta**2 

516 sca2 = _1 + ta2 

517 txia = _txif(ta) 

518 s3qx = _sqrt3(sca2 / (txia**2 + _1)) * qx # * _1_x21(txia) 

519 eta2 = (_1 - e2 * ta2 / sca2)**2 

520 ta, d = _Ta2((txi - txia) * s3qx * eta2) 

521 if _fabs(d) < tol: 

522 return ta 

523 raise AlbersError(Fmt.no_convergence(d, tol), txt=repr(self)) 

524 

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

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

527 

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

529 

530 @return: This projection as C{"<classname>(lat1, lat2, ...)"} 

531 (C{str}). 

532 ''' 

533 t = self.toStr(prec=prec, sep=_COMMASPACE_) 

534 return Fmt.PAREN(self.classname, t) 

535 

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

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

538 

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

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

541 

542 @return: This projection as C{"lat1 lat2"} (C{str}). 

543 ''' 

544 k = self._k 

545 t = (self.lat1, self.lat2, self._k0) if k is _k1_ else ( 

546 (self.lat1, self._k0) if k is _k0_ else 

547 (self.lat1,)) 

548 t = strs(t, prec=prec) 

549 if k: 

550 t = t[:-1] + (Fmt.EQUAL(k, t[-1]),) 

551 if self.datum != _WGS84: 

552 t += (Fmt.EQUAL(_datum_, self.datum),) 

553 if self.name: 

554 t += (Fmt.EQUAL(_name_, repr(self.name)),) 

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

556 

557 def _txif(self, ta): # in .Ellipsoid.auxAuthalic 

558 '''(INTERNAL) Function M{tan-xi from tan-phi}. 

559 ''' 

560 E = self.ellipsoid 

561 _1 = _1_0 

562 

563 ca2 = _1_x21(ta) 

564 sa = sqrt(ca2) * fabs(ta) # enforce odd parity 

565 sa1 = _1 + sa 

566 

567 es1 = sa * E.e2 

568 es1m1 = sa1 * (_1 - es1) 

569 es1p1 = sa1 / (_1 + es1) 

570 es2m1 = _1 - sa * es1 

571 es2m1a = es2m1 * E.e21 # e2m 

572 s = sqrt((ca2 / (es1p1 * es2m1a) + _atanheE(ca2 / es1m1, E)) 

573 * (es1m1 / es2m1a + _atanheE(es1p1, E))) 

574 t = _Fsum1(sa / es2m1, _atanheE(sa, E)).fover(s) 

575 return neg(t) if ta < 0 else t 

576 

577 

578class AlbersEqualArea(_AlbersBase): 

579 '''An Albers equal-area (authalic) projection with a single standard parallel. 

580 

581 @see: L{AlbersEqualArea2} and L{AlbersEqualArea4}. 

582 ''' 

583 _k = _k0_ 

584 

585 def __init__(self, lat, k0=1, datum=_WGS84, name=NN): 

586 '''New L{AlbersEqualArea} projection. 

587 

588 @arg lat: Standard parallel (C{degrees}). 

589 @kwarg k0: Azimuthal scale on the standard parallel (C{scalar}). 

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

591 L{Ellipsoid2} or L{a_f2Tuple}). 

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

593 

594 @raise AlbersError: Invalid B{C{lat}}, B{C{k0}} or no convergence. 

595 ''' 

596 self._lat1 = self._lat2 = lat = _Lat(lat1=lat) 

597 args = tuple(sincos2d(lat)) * 2 + (_Ks(k0=k0), datum, name) 

598 _AlbersBase.__init__(self, *args) 

599 

600 

601class AlbersEqualArea2(_AlbersBase): 

602 '''An Albers equal-area (authalic) projection with two standard parallels. 

603 

604 @see: L{AlbersEqualArea} and L{AlbersEqualArea4}. 

605 ''' 

606 _k = _k1_ 

607 

608 def __init__(self, lat1, lat2, k1=1, datum=_WGS84, name=NN): 

609 '''New L{AlbersEqualArea2} projection. 

610 

611 @arg lat1: First standard parallel (C{degrees}). 

612 @arg lat2: Second standard parallel (C{degrees}). 

613 @kwarg k1: Azimuthal scale on the standard parallels (C{scalar}). 

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

615 L{Ellipsoid2} or L{a_f2Tuple}). 

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

617 

618 @raise AlbersError: Invalid B{C{lat1}}m B{C{lat2}}, B{C{k1}} 

619 or no convergence. 

620 ''' 

621 self._lat1, self._lat2 = lats = _Lat(lat1=lat1), _Lat(lat2=lat2) 

622 args = tuple(sincos2d_(*lats)) + (_Ks(k1=k1), datum, name) 

623 _AlbersBase.__init__(self, *args) 

624 

625 

626class AlbersEqualArea4(_AlbersBase): 

627 '''An Albers equal-area (authalic) projection specified by the C{sin} 

628 and C{cos} of both standard parallels. 

629 

630 @see: L{AlbersEqualArea} and L{AlbersEqualArea2}. 

631 ''' 

632 _k = _k1_ 

633 

634 def __init__(self, slat1, clat1, slat2, clat2, k1=1, datum=_WGS84, name=NN): 

635 '''New L{AlbersEqualArea4} projection. 

636 

637 @arg slat1: Sine of first standard parallel (C{scalar}). 

638 @arg clat1: Cosine of first standard parallel (non-negative C{scalar}). 

639 @arg slat2: Sine of second standard parallel (C{scalar}). 

640 @arg clat2: Cosine of second standard parallel (non-negative C{scalar}). 

641 @kwarg k1: Azimuthal scale on the standard parallels (C{scalar}). 

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

643 L{Ellipsoid2} or L{a_f2Tuple}). 

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

645 

646 @raise AlbersError: Negative B{C{clat1}} or B{C{clat2}}, B{C{slat1}} 

647 and B{C{slat2}} have opposite signs (hemispheres), 

648 invalid B{C{k1}} or no convergence. 

649 ''' 

650 def _Lat_s_c3(name, s, c): 

651 r = Float_(hypot(s, c), name=name, Error=AlbersError) 

652 L = _Lat(atan1d(s, c), name=name) 

653 return L, (s / r), (c / r) 

654 

655 self._lat1, sa1, ca1 = _Lat_s_c3(_lat1_, slat1, clat1) 

656 self._lat2, sa2, ca2 = _Lat_s_c3(_lat2_, slat2, clat2) 

657 _AlbersBase.__init__(self, sa1, ca1, sa2, ca2, _Ks(k1=k1), datum, name) 

658 

659 

660class AlbersEqualAreaCylindrical(_AlbersBase): 

661 '''An L{AlbersEqualArea} projection at C{lat=0} and C{k0=1} degenerating 

662 to the cylindrical-equal-area projection. 

663 ''' 

664 _lat1 = _lat2 = _Lat(lat1=_0_0) 

665 

666 def __init__(self, lat=_0_0, datum=_WGS84, name=NN): 

667 '''New L{AlbersEqualAreaCylindrical} projection. 

668 

669 @kwarg lat: Standard parallel (C{0 degrees} I{fixed}). 

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

671 L{Ellipsoid2} or L{a_f2Tuple}). 

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

673 ''' 

674 _xlat(lat, _0_0, AlbersEqualAreaCylindrical) 

675 _AlbersBase.__init__(self, _0_0, _1_0, _0_0, _1_0, 1, datum, name) 

676 

677 

678class AlbersEqualAreaNorth(_AlbersBase): 

679 '''An azimuthal L{AlbersEqualArea} projection at C{lat=90} and C{k0=1} 

680 degenerating to the L{azimuthal} L{LambertEqualArea} projection. 

681 ''' 

682 _lat1 = _lat2 = _Lat(lat1=_90_0) 

683 

684 def __init__(self, lat=_90_0, datum=_WGS84, name=NN): 

685 '''New L{AlbersEqualAreaNorth} projection. 

686 

687 @kwarg lat: Standard parallel (C{90 degrees} I{fixed}). 

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

689 L{Ellipsoid2} or L{a_f2Tuple}). 

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

691 ''' 

692 _xlat(lat, _90_0, AlbersEqualAreaNorth) 

693 _AlbersBase.__init__(self, _1_0, _0_0, _1_0, _0_0, 1, datum, name) 

694 

695 

696class AlbersEqualAreaSouth(_AlbersBase): 

697 '''An azimuthal L{AlbersEqualArea} projection at C{lat=-90} and C{k0=1} 

698 degenerating to the L{azimuthal} L{LambertEqualArea} projection. 

699 ''' 

700 _lat1 = _lat2 = _Lat(lat1=_N_90_0) 

701 

702 def __init__(self, lat=_N_90_0, datum=_WGS84, name=NN): 

703 '''New L{AlbersEqualAreaSouth} projection. 

704 

705 @kwarg lat: Standard parallel (C{-90 degrees} I{fixed}). 

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

707 L{Ellipsoid2} or L{a_f2Tuple}). 

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

709 ''' 

710 _xlat(lat, _N_90_0, AlbersEqualAreaSouth) 

711 _AlbersBase.__init__(self, _N_1_0, _0_0, _N_1_0, _0_0, 1, datum, name) 

712 

713 

714class Albers7Tuple(_NamedTuple): 

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

716 C{meter}, C{degrees90}, C{degrees180}, C{degrees360}, C{scalar} and 

717 C{Datum} where C{(x, y)} is the projected, C{(lat, lon)} the geodetic 

718 location, C{gamma} the meridian convergence at point, the bearing of 

719 the y-axis measured clockwise from true North and C{scale} is the 

720 azimuthal scale of the projection at point. The radial scale is 

721 the reciprocal C{1 / scale}. 

722 ''' 

723 _Names_ = (_x_, _y_, _lat_, _lon_, _gamma_, _scale_, _datum_) 

724 _Units_ = ( Meter, Meter, Lat, Lon, Bearing, _Pass, _Pass) 

725 

726 

727def _atanh1(x): 

728 '''(INTERNAL) Function M{atanh(sqrt(x)) / sqrt(x) - 1}. 

729 ''' 

730 s = fabs(x) 

731 if 0 < s < _0_5: # for typical ... 

732 # x < E.e^2 == 2 * E.f use ... 

733 # x / 3 + x^2 / 5 + x^3 / 7 + ... 

734 y, k = x, 3 

735 _S2 = Fsum(y / k).fsum2_ 

736 for _ in range(_TERMS): # 9 terms 

737 y *= x # x**n 

738 k += 2 # 2*n + 1 

739 s, d = _S2(y / k) 

740 if not d: 

741 break 

742 elif s: 

743 s = sqrt(s) 

744 s = (atanh(s) if x > 0 else atan1(s)) / s - _1_0 

745 return s 

746 

747 

748def _atanheE(x, E): # see Ellipsoid._es_atanh, .AuxLat._atanhee 

749 '''(INTERNAL) Function M{atanhee(x)}, defined as ... 

750 atanh( E.e * x) / E.e if f > 0 # oblate 

751 atan (sqrt(-E.e2) * x) / sqrt(-E.e2) if f < 0 # prolate 

752 x if f = 0. 

753 ''' 

754 e = E.e # == sqrt(E.e2abs) 

755 if e and x: 

756 if E.f > 0: # .isOblate 

757 x = atanh(x * e) / e 

758 elif E.f < 0: # .isProlate 

759 x = atan1(x * e) / e 

760 return x 

761 

762 

763def _DatanheE(x, y, E): # see .rhumb.ekx._DeatanhE 

764 '''(INTERNAL) Function M{Datanhee(x, y)}, defined as 

765 M{atanhee((x - y) / (1 - E.e^2 * x * y)) / (x - y)}. 

766 ''' 

767 e = _1_0 - E.e2 * x * y 

768 if e: 

769 d = x - y 

770 e = (_atanheE(d / e, E) / d) if d else (_1_0 / e) 

771 return e 

772 

773 

774def _D2atanheE(x, y, E): 

775 '''(INTERNAL) Function M{D2atanhee(x, y)}, defined as 

776 M{(Datanhee(1, y) - Datanhee(1, x)) / (y - x)}. 

777 ''' 

778 s, e2 = _0_0, E.e2 

779 if e2: 

780 if ((fabs(x) + fabs(y)) * e2) < _0_5: 

781 e = z = _1_0 

782 k = 1 

783 T = Fsum() # Taylor expansion 

784 _T = T.fsum_ 

785 _C = Fsum().fsum_ 

786 _S2 = Fsum().fsum2_ 

787 for _ in range(_TERMS): # 15 terms 

788 T *= y; p = _T(z); z *= x # PYCHOK ; 

789 T *= y; t = _T(z); z *= x # PYCHOK ; 

790 e *= e2 

791 k += 2 

792 s, d = _S2(_C(p, t) * e / k) 

793 if not d: 

794 break 

795 else: # PYCHOK no cover 

796 s = _1_0 - x 

797 if s: 

798 s = (_DatanheE(_1_0, y, E) - _DatanheE(x, y, E)) / s 

799 return s 

800 

801 

802def _Dsn(x, y, sx, sy): 

803 '''(INTERNAL) Divided differences, defined as M{Df(x, y) = (f(x) - f(y)) / (x - y)} 

804 with M{sn(x) = x / sqrt(1 + x^2)}: M{Dsn(x, y) = (x + y) / ((sn(x) + sn(y)) * 

805 (1 + x^2) * (1 + y^2))}. 

806 

807 @see: U{W. M. Kahan and R. J. Fateman, "Sympbolic Computation of Divided 

808 Differences"<https://People.EECS.Berkeley.EDU/~fateman/papers/divdiff.pdf>}, 

809 U{ACM SIGSAM Bulletin 33(2), 7-28 (1999)<https://DOI.org/10.1145/334714.334716>} 

810 and U{AlbersEqualArea.hpp 

811 <https://GeographicLib.SourceForge.io/C++/doc/AlbersEqualArea_8hpp_source.html>}. 

812 ''' 

813 # sx = x / hypot1(x) 

814 d, t = _1_0, (x * y) 

815 if t > 0: 

816 s = sx + sy 

817 if s: 

818 t = sx * sy / t 

819 d = t**2 * (x + y) / s 

820 elif x != y: 

821 d = (sx - sy) / (x - y) 

822 return d 

823 

824 

825def _tol(tol, x): 

826 '''(INTERNAL) Converge tolerance. 

827 ''' 

828 return tol * max(_1_0, fabs(x)) 

829 

830 

831def _1_x21(x): 

832 '''(INTERNAL) Return M{1 / (x**2 + 1)}. 

833 ''' 

834 return _1_0 / (x**2 + _1_0) 

835 

836 

837def _xlat(lat, f, where): 

838 '''(INTERNAL) check fixed C{lat}. 

839 ''' 

840 if lat is not f and _Lat(lat=lat) != f: 

841 t = unstr(where.__name__, lat=lat) 

842 raise AlbersError(t, txt=_not_(f)) 

843 

844 

845__all__ += _ALL_DOCS(_AlbersBase) 

846 

847# **) MIT License 

848# 

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

850# 

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

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

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

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

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

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

857# 

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

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

860# 

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

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

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

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

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

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

867# OTHER DEALINGS IN THE SOFTWARE.