1

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

3

4 u'''Albers Equal-Area projections.

5

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

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

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

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

10

11 See 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

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

14'''

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

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

17

18 from pygeodesy.basics import neg, neg_

19 from 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

22 from pygeodesy.datums import _ellipsoidal_datum, _WGS84

23 from pygeodesy.errors import _ValueError, _xkwds

24 from pygeodesy.fmath import hypot, hypot1, sqrt3

25 from pygeodesy.fsums import Fmt, Fsum, fsum1_

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

27 _lat2_, _lon_, _scale_, _x_, _y_

28 from pygeodesy.karney import _diff182, _norm180, _signBit

29 from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY

30 from pygeodesy.named import _NamedBase, _NamedTuple, _Pass

31 from pygeodesy.props import deprecated_Property_RO, Property_RO, _update_all

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

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

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

35 sincos2d, sincos2d_

36

37 from 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

48 class 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

56 def _1_x21(x):

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

58 '''

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

60

61

62 def _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

85 def _Ks(**name_k):

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

87 '''

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

89

90

91 def _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

98 def _tol(tol, x):

99 '''(INTERNAL) Converge tolerance.

100 '''

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

102

103

104 class _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

605 class 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

626 class 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

649 class 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

681 class 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

697 class 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

713 class 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

729 class 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.