Coverage for pygeodesy/albers.py: 97%

421 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-21 13:14 -0400

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, isnear0, \ 

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

21 _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, fsum1_ 

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

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

28 _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 atand, atan2d, degrees360, sincos2, \ 

36 sincos2d, sincos2d_ 

37 

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

39 

40__all__ = _ALL_LAZY.albers 

41__version__ = '23.04.21' 

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 _Ks(**name_k): 

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

52 ''' 

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

54 

55 

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

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

58 ''' 

59 kwds = _xkwds(Error_name_lat, Error=AlbersError) 

60 return Lat_(*lat, **kwds) 

61 

62 

63def _qZx(aobj): 

64 '''(INTERNAL) Set C{aobj._qZ} and C{aobj._qx}. 

65 ''' 

66 E = aobj._datum.ellipsoid 

67 aobj._qZ = qZ = _1_0 + E.e21 * _atanhee(_1_0, E) 

68 aobj._qx = qZ / (_2_0 * E.e21) 

69 return qZ 

70 

71 

72class AlbersError(_ValueError): 

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

74 L{AlbersEqualAreaCylindrical}, L{AlbersEqualAreaNorth}, 

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

76 ''' 

77 pass 

78 

79 

80class _AlbersBase(_NamedBase): 

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

82 

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

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

85 ''' 

86 _datum = _WGS84 

87 _k = NN # or _k0_ or _k1_ 

88 _k0 = _Ks(k0=_1_0) 

89 _k0n0 = None # (INTERNAL) k0 * no 

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

91 _k02n0 = None # (INTERNAL) k02 * n0 

92 _lat0 = None # lat origin 

93 _lat1 = None # let 1st parallel 

94 _lat2 = None # lat 2nd parallel 

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

96 _m02 = None # (INTERNAL) cached 

97 _n0 = None # (INTERNAL) cached 

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

99 _polar = False 

100 _qx = None # (INTERNAL) see _qZx 

101 _qZ = None # (INTERNAL) see _qZx 

102 _scxi0_ = None # (INTERNAL) sec(xi) / (qZ...) 

103 _sign = +1 

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

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

106 

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

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

109 ''' 

110 qZ = self._qZ 

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

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

113 qZ = _qZx(self) 

114 elif qZ is None: 

115 qZ = _qZx(_AlbersBase) 

116 if name: 

117 self.name = name 

118 

119 c = min(ca1, ca2) 

120 if _signBit(c): 

121 raise AlbersError(clat1=ca1, clat2=ca2) 

122 polar = c < EPS0 # == 0 

123 E = self.ellipsoid 

124 

125 # determine hemisphere of tangent latitude 

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

127 self._sign = -1 

128 # internally, tangent latitude positive 

129 sa1, sa2 = neg_(sa1, sa2) 

130 if sa1 > sa2: # make phi1 < phi2 

131 sa1, sa2 = sa2, sa1 

132 ca1, ca2 = ca2, ca1 

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

134 raise AlbersError(slat1=sa1, slat2=sa2) 

135 # avoid singularities at poles 

136 ca1, ca2 = max(EPS0, ca1), max(EPS0, ca2) 

137 ta1, ta2 = (sa1 / ca1), (sa2 / ca2) 

138 

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

140 if par1 or polar: 

141 ta0, C = ta2, _1_0 

142 else: 

143 s1_qZ, C = self._s1_qZ_C2(ca1, sa1, ta1, ca2, sa2, ta2) 

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

145 

146 self._lat0 = _Lat(lat0=self._sign * atand(ta0)) 

147 self._m02 = m02 = _1_x21(E.b_a * ta0) 

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

149 if polar: 

150 self._polar = True 

151# self._nrho0 = self._m0 = _0_0 

152 else: 

153 self._m0 = sqrt(m02) # == nrho0 / E.a 

154 self._nrho0 = E.a * self._m0 # == E.a * sqrt(m02) 

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

156 h = hypot1(t) 

157 s = self._sxi0 = t / h 

158 if par1: 

159 self._k0n0 = self._k02n0 = n0 

160 else: 

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

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

163 

164 def _a_b_sxi3(self, *ca_sa_ta_scb_4s): 

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

166 ''' 

167 _1 = _1_0 

168 a = b = s = _0_0 

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

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

171 if sa > 0: 

172 sa += _1 

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

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

175 else: 

176 sa = _1 - sa 

177 a += (_1 - sxi) / sa 

178 b += scb * sa 

179 s += sxi 

180 return a, b, s 

181 

182 def _azik(self, t, ta): 

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

184 ''' 

185 E = self.ellipsoid 

186 t *= self._k0 

187 return _Ks(k=t * hypot1(E.b_a * ta) / E.a) 

188 

189 def _cstxif3(self, ta): 

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

191 ''' 

192 t = self._txif(ta) 

193 c = _1_0 / hypot1(t) 

194 s = c * t 

195 return c, s, t 

196 

197 @Property_RO 

198 def datum(self): 

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

200 ''' 

201 return self._datum 

202 

203 @Property_RO 

204 def ellipsoid(self): 

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

206 ''' 

207 return self.datum.ellipsoid 

208 

209 @Property_RO 

210 def equatoradius(self): 

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

212 ''' 

213 return self.ellipsoid.a 

214 

215 @Property_RO 

216 def flattening(self): 

217 '''Get the geodesic's flattening (C{float}). 

218 ''' 

219 return self.ellipsoid.f 

220 

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

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

223 

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

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

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

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

228 

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

230 

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

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

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

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

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

236 ''' 

237 a = self.ellipsoid.a 

238 s = self._sign 

239 

240 k0 = self._k0 

241 n0 = self._n0 

242 nrho0 = self._nrho0 

243 txi0 = self._txi0 

244 

245 sa, ca = sincos2d(s * _Lat(lat=lat)) 

246 ta = sa / max(EPS0, ca) 

247 

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

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

250 (txi - txi0) * self._qZ 

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

252 

253 lon, _ = _diff182(lon0, lon) 

254 x = radians(lon) 

255 

256 th = self._k02n0 * x 

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

258 if n0: 

259 x = sth / n0 

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

261 y *= nrho0 / n0 

262 else: 

263 x *= self._k02 

264 y = _0_0 

265 t = nrho0 - n0 * drho 

266 x = t * x / k0 

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

268 

269 g = degrees360(s * th) 

270 if t: 

271 k0 = self._azik(t, ta) 

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

273 name=name or self.name) 

274 

275 @Property_RO 

276 def ispolar(self): 

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

278 ''' 

279 return self._polar 

280 

281 isPolar = ispolar # synonym 

282 

283 def _k0s(self, k0): 

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

285 ''' 

286 self._k0 = k = _Ks(k0=k0) 

287 self._k02 = k2 = k**2 

288 self._k0n0 = k * self._n0 

289 self._k02n0 = k2 * self._n0 

290 

291 @Property_RO 

292 def lat0(self): 

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

294 

295 This is the latitude of minimum azimuthal scale and 

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

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

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

299 projections. 

300 ''' 

301 return self._lat0 

302 

303 @Property_RO 

304 def lat1(self): 

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

306 ''' 

307 return self._lat1 

308 

309 @Property_RO 

310 def lat2(self): 

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

312 

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

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

315 projections. 

316 ''' 

317 return self._lat2 

318 

319 @deprecated_Property_RO 

320 def majoradius(self): # PYCHOK no cover 

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

322 return self.equatoradius 

323 

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

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

326 

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

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

329 

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

331 

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

333 ''' 

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

335 if self._k0 != k0: 

336 _update_all(self) 

337 self._k0s(k0) 

338 

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

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

341 

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

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

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

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

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

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

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

349 

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

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

352 

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

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

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

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

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

358 space the nearest pole is returned. 

359 ''' 

360 k0 = self._k0 

361 n0 = self._n0 

362 k0n0 = self._k0n0 

363 s = self._sign 

364 txi = self._txi0 

365 

366 x = Meter(x=x) 

367 nx = k0n0 * x 

368 y = Meter(y=y) 

369 y_ = s * y 

370 ny = k0n0 * y_ 

371 t = nrho0 = self._nrho0 

372 y1 = nrho0 - ny 

373 

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

375 if den: 

376 drho = fsum1_(x * nx, y_ * nrho0 * _N_2_0, y_ * ny) * k0 / den 

377 # dsxia = scxi0 * dsxi 

378 t += drho * n0 

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

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

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

382 

383 ta = self._tanf(txi) 

384 lat = atand(s * ta) 

385 

386 th = atan2(nx, y1) 

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

388 if lon0: 

389 lon += _norm180(lon0) 

390 lon = _norm180(lon) 

391 

392 n = name or self.name 

393 if LatLon is None: 

394 g = degrees360(s * th) 

395 if den: 

396 k0 = self._azik(t, ta) 

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

398 else: # PYCHOK no cover 

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

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

401 return r 

402 

403 @Property_RO 

404 def scale0(self): 

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

406 

407 This is the azimuthal scale on the latitude of origin 

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

409 ''' 

410 return self._k0 

411 

412 def _s1_qZ_C2(self, ca1, sa1, ta1, ca2, sa2, ta2): 

413 '''(INTERNAL) Compute C{sm1 / (s / qZ)} and C{C} for .__init__. 

414 ''' 

415 _1 = _1_0 

416 _fsum1 = fsum1_ 

417 E = self.ellipsoid 

418 f1, e2 = E.b_a, E.e2 

419 

420 tb1 = f1 * ta1 

421 tb2 = f1 * ta2 

422 dtb12 = f1 * (tb1 + tb2) 

423 scb12 = _1 + tb1**2 

424 scb22 = _1 + tb2**2 

425 

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

427 sa12 = sa1 * sa2 

428 

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

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

431 esa12 = _1 + e2 * sa12 

432 

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

434 (ca2, sa2, ta2, scb22)) 

435 

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

437 C = _fsum1(sxi * dtb12 / dsxi, scb22, scb12) / (scb22 * scb12 * _2_0) 

438 

439 sa12 = _fsum1(sa1, sa2, sa12) 

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

441 bxi *= _fsum1(sa1, sa2, esa12) * e2 / esa1_2 + E.e21 * _D2atanhee(sa1, sa2, E) 

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

443 return s1_qZ, C 

444 

445 def _ta0(self, s1_qZ, ta0): 

446 '''(INTERNAL) Refine C{ta0}. 

447 ''' 

448 e2 = self.ellipsoid.e2 

449 e21 = self.ellipsoid.e21 

450 tol = _tol(_TOL0, ta0) 

451 _Ta02 = Fsum(ta0).fsum2_ 

452 _fabs = fabs 

453 _fsum1 = fsum1_ 

454 _sqrt = sqrt 

455 _1, _2 = _1_0, _2_0 

456 _4, _6 = _4_0, _6_0 

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

458 ta02 = ta0**2 

459 sca02 = ta02 + _1 

460 sca0 = _sqrt(sca02) 

461 sa0 = ta0 / sca0 

462 sa01 = sa0 + _1 

463 sa02 = sa0**2 

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

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

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

467 sa021 = _1 - e2 * sa02 

468 

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

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

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

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

473 BA = (_atanhs1(e2 * sa0m1**2) * e21 - e2 * sa0m) * sa0m1 \ 

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

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

476 u = _fsum1(s1_qZ * g, -D, g * BA, floats=True) 

477 du = _fsum1(s1_qZ * dg, dD, dg * BA, g * d, floats=True) 

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

479 if _fabs(d) < tol: 

480 return ta0 

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

482 

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

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

485 ''' 

486 tol = _tol(_TOL, txi) 

487 e2 = self.ellipsoid.e2 

488 qx = self._qx 

489 

490 ta = txi 

491 _Ta2 = Fsum(ta).fsum2_ 

492 _fabs = fabs 

493 _sqrt3 = sqrt3 

494 _txif = self._txif 

495 _1 = _1_0 

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

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

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

499 ta2 = ta**2 

500 sca2 = _1 + ta2 

501 txia = _txif(ta) 

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

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

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

505 if _fabs(d) < tol: 

506 return ta 

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

508 

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

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

511 

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

513 

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

515 (C{str}). 

516 ''' 

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

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

519 

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

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

522 

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

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

525 

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

527 ''' 

528 k = self._k 

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

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

531 (self.lat1,)) 

532 t = strs(t, prec=prec) 

533 if k: 

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

535 if self.datum != _WGS84: 

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

537 if self.name: 

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

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

540 

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

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

543 ''' 

544 E = self.ellipsoid 

545 

546 ca2 = _1_x21(ta) 

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

548 sa1 = _1_0 + sa 

549 

550 es1 = sa * E.e2 

551 es1m1 = sa1 * (_1_0 - es1) 

552 es1p1 = sa1 / (_1_0 + es1) 

553 es2m1 = _1_0 - sa * es1 

554 es2m1a = es2m1 * E.e21 # e2m 

555 s = sqrt((ca2 / (es1p1 * es2m1a) + _atanhee(ca2 / es1m1, E)) 

556 * (es1m1 / es2m1a + _atanhee(es1p1, E))) 

557 t = (sa / es2m1 + _atanhee(sa, E)) / s 

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

559 

560 

561class AlbersEqualArea(_AlbersBase): 

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

563 

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

565 ''' 

566 _k = _k0_ 

567 

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

569 '''New L{AlbersEqualArea} projection. 

570 

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

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

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

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

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

576 

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

578 ''' 

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

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

581 _AlbersBase.__init__(self, *args) 

582 

583 

584class AlbersEqualArea2(_AlbersBase): 

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

586 

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

588 ''' 

589 _k = _k1_ 

590 

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

592 '''New L{AlbersEqualArea2} projection. 

593 

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

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

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

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

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

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

600 

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

602 or no convergence. 

603 ''' 

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

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

606 _AlbersBase.__init__(self, *args) 

607 

608 

609class AlbersEqualArea4(_AlbersBase): 

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

611 and C{cos} of both standard parallels. 

612 

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

614 ''' 

615 _k = _k1_ 

616 

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

618 '''New L{AlbersEqualArea4} projection. 

619 

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

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

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

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

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

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

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

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

628 

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

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

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

632 ''' 

633 def _Lat_s_c3(name, s, c): 

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

635 L = _Lat(atan2d(s, c), name=name) 

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

637 

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

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

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

641 

642 

643class AlbersEqualAreaCylindrical(_AlbersBase): 

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

645 to the cylindrical-equal-area projection. 

646 ''' 

647 _lat1 = _lat2 = _Lat(lat1=_0_0) 

648 

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

650 '''New L{AlbersEqualAreaCylindrical} projection. 

651 

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

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

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

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

656 ''' 

657 _xlat(lat, _0_0, AlbersEqualAreaCylindrical) 

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

659 

660 

661class AlbersEqualAreaNorth(_AlbersBase): 

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

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

664 ''' 

665 _lat1 = _lat2 = _Lat(lat1=_90_0) 

666 

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

668 '''New L{AlbersEqualAreaNorth} projection. 

669 

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

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

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

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

674 ''' 

675 _xlat(lat, _90_0, AlbersEqualAreaNorth) 

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

677 

678 

679class AlbersEqualAreaSouth(_AlbersBase): 

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

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

682 ''' 

683 _lat1 = _lat2 = _Lat(lat1=_N_90_0) 

684 

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

686 '''New L{AlbersEqualAreaSouth} projection. 

687 

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

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

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

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

692 ''' 

693 _xlat(lat, _N_90_0, AlbersEqualAreaSouth) 

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

695 

696 

697class Albers7Tuple(_NamedTuple): 

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

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

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

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

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

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

704 the reciprocal C{1 / scale}. 

705 ''' 

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

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

708 

709 

710def _atanhee(x, E): # see Ellipsoid._es_atanh 

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

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

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

714 x if f = 0. 

715 ''' 

716 if E.f > 0: # .isOblate 

717 x = atanh( x * E.e) / E.e 

718 elif E.f < 0: # .isProlate 

719 x = (atan2(-x * E.e, _N_1_0) if x < 0 else 

720 atan2( x * E.e, _1_0)) / E.e 

721 return x 

722 

723 

724def _atanhs1(x): 

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

726 ''' 

727 s = fabs(x) 

728 if s < _0_5: # for typical ... 

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

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

731 y, k = x, 3 

732 _S2 = Fsum(y / k).fsum2_ 

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

734 y *= x # x**n 

735 k += 2 # 2*n + 1 

736 s, d = _S2(y / k) 

737 if not d: 

738 break 

739 else: 

740 s = sqrt(s) 

741 s = (atanh(s) if x > 0 else atan(s)) / s - _1_0 

742 return s 

743 

744 

745def _Datanhee(x, y, E): # see .rhumbx._DeatanhE 

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

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

748 ''' 

749 d = _1_0 - E.e2 * x * y 

750 if d: 

751 d = _1_0 / d 

752 t = x - y 

753 if t: 

754 d = _atanhee(t * d, E) / t 

755 else: 

756 raise AlbersError(x=x, y=y, txt=_Datanhee.__name__) 

757 return d 

758 

759 

760def _D2atanhee(x, y, E): 

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

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

763 ''' 

764 e2 = E.e2 

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

766 e = z = _1_0 

767 k = 1 

768 T = Fsum() # Taylor expansion 

769 _T = T.fsum_ 

770 _C = Fsum().fsum_ 

771 _S2 = Fsum().fsum2_ 

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

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

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

775 e *= e2 

776 k += 2 

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

778 if not d: 

779 break 

780 else: # PYCHOK no cover 

781 d = _1_0 - x 

782 if isnear0(d): 

783 raise AlbersError(x=x, y=y, txt=_D2atanhee.__name__) 

784 s = (_Datanhee(_1_0, y, E) - _Datanhee(x, y, E)) / d 

785 return s 

786 

787 

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

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

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

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

792 

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

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

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

796 and U{AlbersEqualArea.hpp 

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

798 ''' 

799 # sx = x / hypot1(x) 

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

801 if t > 0: 

802 s = sx + sy 

803 if s: 

804 t = sx * sy / t 

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

806 elif x != y: 

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

808 return d 

809 

810 

811def _tol(tol, x): 

812 '''(INTERNAL) Converge tolerance. 

813 ''' 

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

815 

816 

817def _1_x21(x): 

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

819 ''' 

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

821 

822 

823def _xlat(lat, f, where): 

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

825 ''' 

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

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

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

829 

830 

831__all__ += _ALL_DOCS(_AlbersBase) 

832 

833# **) MIT License 

834# 

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

836# 

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

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

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

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

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

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

843# 

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

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

846# 

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

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

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

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

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

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

853# OTHER DEALINGS IN THE SOFTWARE.