Coverage for pygeodesy/albers.py: 97%

375 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-05 15:46 -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, _0_0, \ 

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

21 _3_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 Fmt, Fsum, fsum1_ 

26from pygeodesy.interns import NN, _datum_, _gamma_, _lat_, _lat1_, \ 

27 _lat2_, _lon_, _scale_, _x_, _y_ 

28from pygeodesy.karney import _diff182, _norm180, _signBit 

29from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY 

30from pygeodesy.named import _NamedBase, _NamedTuple, _Pass 

31from pygeodesy.props import deprecated_Property_RO, Property_RO, _update_all 

32# from pygeodesy.streprs import Fmt # from .fsums 

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

34from pygeodesy.utily import atand, atan2d, degrees360, sincos2, \ 

35 sincos2d, sincos2d_ 

36 

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

38 

39__all__ = _ALL_LAZY.albers 

40__version__ = '23.03.19' 

41 

42_NUMIT = 8 # XXX 4? 

43_NUMIT0 = 41 # XXX 21? 

44_TERMS = 21 # XXX 16? 

45_TOL0 = sqrt3(_TOL) 

46 

47 

48class AlbersError(_ValueError): 

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

50 L{AlbersEqualAreaCylindrical}, L{AlbersEqualAreaNorth}, 

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

52 ''' 

53 pass 

54 

55 

56def _1_x21(x): 

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

58 ''' 

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

60 

61 

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

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

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

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

66 

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

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

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

70 and U{AlbersEqualArea.hpp 

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

72 ''' 

73 # sx = x / hypot1(x) 

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

75 if t > 0: 

76 s = sx + sy 

77 if s: 

78 t = sx * sy / t 

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

80 elif x != y: 

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

82 return d 

83 

84 

85def _Ks(**name_k): 

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

87 ''' 

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

89 

90 

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

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

93 ''' 

94 kwds = _xkwds(Error_name_lat, Error=AlbersError) 

95 return Lat_(*lat, **kwds) 

96 

97 

98def _tol(tol, x): 

99 '''(INTERNAL) Converge tolerance. 

100 ''' 

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

102 

103 

104class _AlbersBase(_NamedBase): 

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

106 

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

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

109 ''' 

110 _datum = _WGS84 

111 _k0 = None # scale 

112 _k0n0 = None # (INTERNAL) k0 * no 

113 _k02 = None # (INTERNAL) k0**2 

114 _k02n0 = None # (INTERNAL) k02 * n0 

115 _lat0 = None # lat origin 

116 _lat1 = None # let 1st parallel 

117 _lat2 = None # lat 2nd parallel 

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

119 _m02 = None # (INTERNAL) cached 

120 _n0 = None # (INTERNAL) cached 

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

122 _polar = False 

123 _qx = None # (INTERNAL) cached 

124 _qZ = None # (INTERNAL) cached 

125 _qZa2 = None # (INTERNAL) qZ * E.a**2 

126 _scxi0 = None # (INTERNAL) sec(xi) 

127 _sign = +1 

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

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

130 

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

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

133 ''' 

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

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

136 if name: 

137 self.name = name 

138 

139 E = self.datum.ellipsoid 

140 b_a = E.b_a # fm = 1 - E.f 

141 e2 = E.e2 

142 e21 = E.e21 # e2m = 1 - E.e2 

143 

144 self._qZ = qZ = _1_0 + e21 * self._atanhee(1) 

145 self._qZa2 = qZ * E.a2 

146 self._qx = qZ / (_2_0 * e21) 

147 

148 c = min(ca1, ca2) 

149 if _signBit(c): 

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

151 polar = c < EPS0 # == 0 

152 # determine hemisphere of tangent latitude 

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

154 self._sign = -1 

155 # internally, tangent latitude positive 

156 sa1, sa2 = neg_(sa1, sa2) 

157 if sa1 > sa2: # make phi1 < phi2 

158 sa1, sa2 = sa2, sa1 

159 ca1, ca2 = ca2, ca1 

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

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

162 # avoid singularities at poles 

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

164 ta1, ta2 = sa1 / ca1, sa2 / ca2 

165 

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

167 if par1 or polar: 

168 C, ta0 = _1_0, ta2 

169 else: 

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

171 

172 ta0 = (ta2 + ta1) * _0_5 

173 Ta02_ = Fsum(ta0).fsum2_ 

174 tol = _tol(_TOL0, ta0) 

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

176 ta02 = ta0**2 

177 sca02 = ta02 + _1_0 

178 sca0 = sqrt(sca02) 

179 sa0 = ta0 / sca0 

180 sa01 = sa0 + _1_0 

181 sa02 = sa0**2 

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

183 sa0m = _1_0 / (sca0 * (ta0 + sca0)) # scb0^2 * sa0 

184 g = (_1_0 + (b_a * ta0)**2) * sa0 

185 dg = e21 * sca02 * (_1_0 + 2 * ta02) + e2 

186 D = sa0m * (_1_0 - e2 * (_1_0 + sa01 * 2 * sa0)) / (e21 * sa01) # dD/dsa0 

187 dD = _2_0 * (_1_0 - e2 * sa02 * (_3_0 + 2 * sa0)) / (e21 * sa01**2) 

188 sa02_ = _1_0 - e2 * sa02 

189 sa0m_ = sa0m / (_1_0 - e2 * sa0) 

190 BA = sa0m_ * (self._atanhx1(e2 * sa0m_**2) * e21 - e2 * sa0m) \ 

191 - sa0m**2 * e2 * (2 + (_1_0 + e2) * sa0) / (e21 * sa02_) # == B + A 

192 dAB = 2 * e2 * (2 - e2 * (_1_0 + sa02)) / (e21 * sa02_**2 * sca02) 

193 u_du = fsum1_(s1_qZ * g, -D, g * BA) \ 

194 / fsum1_(s1_qZ * dg, dD, dg * BA, g * dAB, floats=True) # == u/du 

195 ta0, d = Ta02_(-u_du * (sca0 * sca02)) 

196 if fabs(d) < tol: 

197 break 

198 else: 

199 raise AlbersError(Fmt.no_convergence(d, tol), txt=str(self)) 

200 

201 self._txi0 = txi0 = self._txif(ta0) 

202 self._scxi0 = hypot1(txi0) 

203 self._sxi0 = sxi0 = txi0 / self._scxi0 

204 self._m02 = m02 = _1_x21(b_a * ta0) 

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

206 if polar: 

207 self._polar = True 

208 self._nrho0 = self._m0 = _0_0 

209 else: 

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

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

212 self._k0s(_1_0 if par1 else (k * sqrt(C / (m02 + n0 * qZ * sxi0)))) 

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

214 

215 def _a_b_sxi3(self, *ca_sa_ta_scb_4s): 

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

217 ''' 

218 a = b = s = _0_0 

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

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

221 if sa > 0: 

222 sa += _1_0 

223 a += (cxi / ca)**2 * sa / (sxi + _1_0) 

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

225 else: 

226 sa = _1_0 - sa 

227 a += (_1_0 - sxi) / sa 

228 b += scb * sa 

229 s += sxi 

230 return a, b, s 

231 

232 def _atanhee(self, x): 

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

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

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

236 x if f = 0. 

237 ''' 

238 E = self.datum.ellipsoid 

239 if E.f > 0: # .isOblate 

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

241 elif E.f < 0: # .isProlate 

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

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

244# else: 

245# pass # x = x 

246 return x 

247 

248 def _atanhx1(self, x): 

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

250 ''' 

251 s = fabs(x) 

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

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

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

255 y, k = x, 3 

256 S2_ = Fsum(y / k).fsum2_ 

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

258 y *= x # x**n 

259 k += 2 # 2*n + 1 

260 s, d = S2_(y / k) 

261 if not d: 

262 break 

263 else: 

264 s = sqrt(s) 

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

266 return s 

267 

268 def _azik(self, t, ta): 

269 '''(INTERNAL) Compute the azimuthal scale C{_Ks(k0=k0)}. 

270 ''' 

271 E = self.datum.ellipsoid 

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

273 

274 def _cstxif3(self, ta): 

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

276 ''' 

277 t = self._txif(ta) 

278 c = _1_0 / hypot1(t) 

279 s = c * t 

280 return c, s, t 

281 

282 def _Datanhee(self, x, y): 

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

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

285 ''' 

286 e2 = self.datum.ellipsoid.e2 

287 d = _1_0 - e2 * x * y 

288 if d: 

289 d = _1_0 / d 

290 t = x - y 

291 if t: 

292 d = self._atanhee(t * d) / t 

293 else: 

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

295 return d 

296 

297 def _D2atanhee(self, x, y): 

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

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

300 ''' 

301 e2 = self.datum.ellipsoid.e2 

302 

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

304 e = z = _1_0 

305 k = 1 

306 T = Fsum() # Taylor expansion 

307 _T = T.fsum_ 

308 _C = Fsum().fsum_ 

309 _S2 = Fsum().fsum2_ 

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

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

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

313 e *= e2 

314 k += 2 

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

316 if not d: 

317 break 

318 elif (_1_0 - x): 

319 s = (self._Datanhee(_1_0, y) - self._Datanhee(x, y)) / (_1_0 - x) 

320 else: 

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

322 return s 

323 

324 @Property_RO 

325 def datum(self): 

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

327 ''' 

328 return self._datum 

329 

330 @Property_RO 

331 def equatoradius(self): 

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

333 ''' 

334 return self.datum.ellipsoid.a 

335 

336 @Property_RO 

337 def flattening(self): 

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

339 ''' 

340 return self.datum.ellipsoid.f 

341 

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

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

344 

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

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

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

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

349 

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

351 

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

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

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

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

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

357 ''' 

358 E = self.datum.ellipsoid 

359 s = self._sign 

360 

361 k0 = self._k0 

362 n0 = self._n0 

363 nrho0 = self._nrho0 

364 txi0 = self._txi0 

365 

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

367 ca = max(EPS0, ca) 

368 ta = sa / ca 

369 

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

371 dq = self._qZ * _Dsn(txi, txi0, sxi, self._sxi0) * (txi - txi0) 

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

373 

374 lon, _ = _diff182(lon0, lon) 

375 b = radians(lon) 

376 

377 th = self._k02n0 * b 

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

379 if n0: 

380 x = sth / n0 

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

382 y *= nrho0 / n0 

383 else: 

384 x = self._k02 * b 

385 y = _0_0 

386 t = nrho0 + n0 * drho 

387 x = t * x / k0 

388 y = s * (y - drho * cth) / k0 

389 

390 g = degrees360(s * th) 

391 if t: 

392 k0 = self._azik(t, ta) 

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

394 name=name or self.name) 

395 

396 @Property_RO 

397 def ispolar(self): 

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

399 ''' 

400 return self._polar 

401 

402 isPolar = ispolar # synonym 

403 

404 def _k0s(self, k): 

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

406 ''' 

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

408 self._k02 = k2 = k**2 

409 self._k0n0 = k * self._n0 

410 self._k02n0 = k2 * self._n0 

411 

412 @Property_RO 

413 def lat0(self): 

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

415 

416 This is the latitude of minimum azimuthal scale and 

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

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

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

420 projections. 

421 ''' 

422 return self._lat0 

423 

424 @Property_RO 

425 def lat1(self): 

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

427 ''' 

428 return self._lat1 

429 

430 @Property_RO 

431 def lat2(self): 

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

433 

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

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

436 projections. 

437 ''' 

438 return self._lat2 

439 

440 @deprecated_Property_RO 

441 def majoradius(self): # PYCHOK no cover 

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

443 return self.equatoradius 

444 

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

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

447 

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

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

450 

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

452 

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

454 ''' 

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

456 if self._k0 != k0: 

457 _update_all(self) 

458 self._k0s(k0) 

459 

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

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

462 

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

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

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

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

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

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

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

470 

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

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

473 

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

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

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

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

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

479 space the nearest pole is returned. 

480 ''' 

481 x = Meter(x=x) 

482 y = Meter(y=y) 

483 

484 k0 = self._k0 

485 n0 = self._n0 

486 k0n0 = self._k0n0 

487 nrho0 = self._nrho0 

488 txi0 = self._txi0 

489 

490 y_ = self._sign * y 

491 nx = k0n0 * x 

492 ny = k0n0 * y_ 

493 y1 = nrho0 - ny 

494 

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

496 if den: 

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

498 # dsxia = scxi0 * dsxi 

499 dsxia = -self._scxi0 * (nrho0 * _2_0 + n0 * drho) * drho / self._qZa2 

500 t = _1_0 - dsxia * (txi0 * _2_0 + dsxia) 

501 txi = (txi0 + dsxia) / (sqrt(t) if t > EPS02 else EPS0) 

502 

503 ta = self._tanf(txi) 

504 lat = atand(ta * self._sign) 

505 

506 th = atan2(nx, y1) 

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

508 if lon0: 

509 lon += _norm180(lon0) 

510 lon = _norm180(lon) 

511 

512 n = name or self.name 

513 if LatLon is None: 

514 g = degrees360(self._sign * th) 

515 if den: 

516 k0 = self._azik(nrho0 + n0 * drho, ta) 

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

518 else: # PYCHOK no cover 

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

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

521 return r 

522 

523 @Property_RO 

524 def scale0(self): 

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

526 

527 This is the azimuthal scale on the latitude of origin 

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

529 ''' 

530 return self._k0 

531 

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

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

534 ''' 

535 E = self.datum.ellipsoid 

536 b_a = E.b_a 

537 e2 = E.e2 

538 

539 tb1 = b_a * ta1 

540 tb2 = b_a * ta2 

541 dtb12 = b_a * (tb1 + tb2) 

542 scb12 = _1_0 + tb1**2 

543 scb22 = _1_0 + tb2**2 

544 

545 esa1_2 = (_1_0 - e2 * sa1**2) \ 

546 * (_1_0 - e2 * sa2**2) 

547 esa12 = _1_0 + e2 * sa1 * sa2 

548 

549 dsn = _Dsn(ta2, ta1, sa2, sa1) 

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

551 (ca2, sa2, ta2, scb22)) 

552 

553 dsxi = (esa12 / esa1_2 + self._Datanhee(sa2, sa1)) * dsn / (self._qx * _2_0) 

554 C = fsum1_(sxi * dtb12 / dsxi, scb22, scb12) / (scb22 * scb12 * _2_0) 

555 

556 sa12 = fsum1_(sa1, sa2, sa1 * sa2) 

557 axi *= (e2 * sa12 + _1_0) / (sa12 + _1_0) 

558 bxi *= e2 * fsum1_(sa1, sa2, esa12) / esa1_2 + E.e21 * self._D2atanhee(sa1, sa2) 

559 s1_qZ = dsn * (axi * self._qZ - bxi) / (dtb12 * _2_0) 

560 return s1_qZ, C 

561 

562 def _tanf(self, txi): # called by .Ellipsoid.auxAuthalic 

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

564 ''' 

565 tol = _tol(_TOL, txi) 

566 

567 e2 = self.datum.ellipsoid.e2 

568 qx = self._qx 

569 

570 ta = txi 

571 _Ta2 = Fsum(ta).fsum2_ 

572 _txif = self._txif 

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

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

575 ta2 = ta**2 

576 sca2 = ta2 + _1_0 

577 eta2 = (_1_0 - e2 * ta2 / sca2)**2 

578 txia = _txif(ta) 

579 s3qx = sqrt3(sca2 * _1_x21(txia)) * qx 

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

581 if fabs(d) < tol: 

582 return ta 

583 raise AlbersError(Fmt.no_convergence(d, tol), txt=str(self)) 

584 

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

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

587 ''' 

588 E = self.datum.ellipsoid 

589 

590 ca2 = _1_x21(ta) 

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

592 

593 es1 = E.e2 * sa 

594 es2m1 = _1_0 - sa * es1 

595 sp1 = _1_0 + sa 

596 es1p1 = sp1 / (_1_0 + es1) 

597 es1m1 = sp1 * (_1_0 - es1) 

598 es2m1a = es2m1 * E.e21 # e2m 

599 s = sqrt((ca2 / (es1p1 * es2m1a) + self._atanhee(ca2 / es1m1)) 

600 * (es1m1 / es2m1a + self._atanhee(es1p1))) 

601 t = (sa / es2m1 + self._atanhee(sa)) / s 

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

603 

604 

605class AlbersEqualArea(_AlbersBase): 

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

607 

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

609 ''' 

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

611 '''New L{AlbersEqualArea} projection. 

612 

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

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

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

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

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

618 

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

620 ''' 

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

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

623 _AlbersBase.__init__(self, *args) 

624 

625 

626class AlbersEqualArea2(_AlbersBase): 

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

628 

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

630 ''' 

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

632 '''New L{AlbersEqualArea2} projection. 

633 

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

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

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

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

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

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

640 

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

642 or no convergence. 

643 ''' 

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

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

646 _AlbersBase.__init__(self, *args) 

647 

648 

649class AlbersEqualArea4(_AlbersBase): 

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

651 and C{cos} of both standard parallels. 

652 

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

654 ''' 

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

656 '''New L{AlbersEqualArea4} projection. 

657 

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

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

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

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

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

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

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

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

666 

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

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

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

670 ''' 

671 def _Lat_s_c3(name, s, c): 

672 L = Lat_(atan2d(s, c), name=name, Error=AlbersError) 

673 r = _1_0 / Float_(hypot(s, c), name=name, Error=AlbersError) 

674 return L, s * r, c * r 

675 

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

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

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

679 

680 

681class AlbersEqualAreaCylindrical(_AlbersBase): 

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

683 to the cylindrical-equal-area projection. 

684 ''' 

685 _lat1 = _lat2 = _Lat(lat1=_0_0) 

686 

687 def __init__(self, datum=_WGS84, name=NN): 

688 '''New L{AlbersEqualAreaCylindrical} projection. 

689 

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

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

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

693 ''' 

694 _AlbersBase.__init__(self, _0_0, _1_0, _0_0, _1_0, _1_0, datum, name) 

695 

696 

697class AlbersEqualAreaNorth(_AlbersBase): 

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

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

700 ''' 

701 _lat1 = _lat2 = _Lat(lat1=_90_0) 

702 

703 def __init__(self, datum=_WGS84, name=NN): 

704 '''New L{AlbersEqualAreaNorth} projection. 

705 

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 _AlbersBase.__init__(self, _1_0, _0_0, _1_0, _0_0, _1_0, datum, name) 

711 

712 

713class AlbersEqualAreaSouth(_AlbersBase): 

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

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

716 ''' 

717 _lat1 = _lat2 = _Lat(lat1=_N_90_0) 

718 

719 def __init__(self, datum=_WGS84, name=NN): 

720 '''New L{AlbersEqualAreaSouth} projection. 

721 

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

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

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

725 ''' 

726 _AlbersBase.__init__(self, _N_1_0, _0_0, _N_1_0, _0_0, _1_0, datum, name) 

727 

728 

729class Albers7Tuple(_NamedTuple): 

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

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

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

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

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

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

736 the reciprocal C{1 / scale}. 

737 ''' 

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

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

740 

741 

742__all__ += _ALL_DOCS(_AlbersBase) 

743 

744# **) MIT License 

745# 

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

747# 

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

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

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

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

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

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

754# 

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

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

757# 

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

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

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

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

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

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

764# OTHER DEALINGS IN THE SOFTWARE.