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

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

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 Fsum, _Fsum1f_, fsum1f_

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

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

28 _negative_, _scale_, _SPACE_, _x_, _y_

29 from pygeodesy.karney import _diff182, _norm180, _signBit

30 from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY

31 from pygeodesy.named import _NamedBase, _NamedTuple, _Pass

32 from pygeodesy.props import deprecated_Property_RO, Property_RO, _update_all

33 from pygeodesy.streprs import Fmt, strs, unstr

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

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

36 sincos2d_

37

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

39

40__all__ = _ALL_LAZY.albers

41__version__ = '24.04.14'

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

50 def _ct2(s, c):

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

52 '''

53 c = max(EPS0, c)

54 return c, (s / c)

55

56

57 def _Ks(**name_k):

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

59 '''

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

61

62

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

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

65 '''

66 kwds = _xkwds(Error_name_lat, Error=AlbersError)

67 return Lat_(*lat, **kwds)

68

69

70 def _qZx(albs):

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

72 '''

73 E = albs._datum.ellipsoid # _AlbersBase

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

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

76 return qZ

77

78

79 class AlbersError(_ValueError):

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

81 L{AlbersEqualAreaCylindrical}, L{AlbersEqualAreaNorth},

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

83 '''

84 pass

85

86

87 class _AlbersBase(_NamedBase):

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

89

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

91 C++/doc/classGeographicLib_1_1AlbersEqualArea.html>}, method C{Init}.

92 '''

93 _datum = _WGS84

94 _k = NN # or _k0_ or _k1_

95 _k0 = _Ks(k0=_1_0)

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

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

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

99# _lat0 = None # lat origin

100 _lat1 = None # let 1st parallel

101 _lat2 = None # lat 2nd parallel

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

103# _m02 = None # (INTERNAL) cached

104# _n0 = None # (INTERNAL) cached

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

106 _polar = False

107 _qx = None # (INTERNAL) see _qZx

108 _qZ = None # (INTERNAL) see _qZx

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

110 _sign = +1

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

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

113

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

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

116 '''

117 qZ = self._qZ

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

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

120 qZ = _qZx(self)

121 elif qZ is None:

122 qZ = _qZx(_AlbersBase)

123 if name:

124 self.name = name

125

126 E = self.ellipsoid

127 c = min(ca1, ca2)

128 if _signBit(c):

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

130 polar = c < EPS0 # == 0

131

132 # determine hemisphere of tangent latitude

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

134 self._sign = -1

135 # internally, tangent latitude positive

136 sa1, sa2 = neg_(sa1, sa2)

137 if sa1 > sa2: # make phi1 < phi2

138 sa1, sa2 = sa2, sa1

139 ca1, ca2 = ca2, ca1

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

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

142 ca1, ta1 = _ct2(sa1, ca1)

143 ca2, ta2 = _ct2(sa2, ca2)

144

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

146 if par1 or polar:

147 ta0, C = ta2, _1_0

148 else:

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

150

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

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

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

154 if polar:

155 self._polar = True

156# self._nrho0 = self._m0 = _0_0

157 else: # m0 = nrho0 / E.a

158 self._m0 = sqrt(m02)

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

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

161 h = hypot1(t)

162 s = self._sxi0 = t / h

163 if par1:

164 self._k0n0 = self._k02n0 = n0

165 else:

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

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

168

169 def _a_b_sxi3(self, *ca_sa_ta_scb_4s):

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

171 '''

172 _1 = _1_0

173 a = b = s = _0_0

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

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

176 if sa > 0:

177 sa += _1

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

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

180 else:

181 sa = _1 - sa

182 a += (_1 - sxi) / sa

183 b += scb * sa

184 s += sxi

185 return a, b, s

186

187 def _azik(self, t, ta):

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

189 '''

190 E = self.ellipsoid

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

192

193 def _cstxif3(self, ta):

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

195 '''

196 t = self._txif(ta)

197 c = _1_0 / hypot1(t)

198 s = c * t

199 return c, s, t

200

201 @Property_RO

202 def datum(self):

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

204 '''

205 return self._datum

206

207 @Property_RO

208 def ellipsoid(self):

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

210 '''

211 return self.datum.ellipsoid

212

213 @Property_RO

214 def equatoradius(self):

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

216 '''

217 return self.ellipsoid.a

218

219 a = equatoradius

220

221 @Property_RO

222 def flattening(self):

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

224 '''

225 return self.ellipsoid.f

226

227 f = flattening

228

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

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

231

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

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

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

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

236

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

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

239

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

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

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

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

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

245 '''

246 a = self.ellipsoid.a

247 s = self._sign

248

249 k0 = self._k0

250 n0 = self._n0

251 nrho0 = self._nrho0

252 txi0 = self._txi0

253

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

255

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

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

258 (txi - txi0) * self._qZ

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

260

261 lon, _ = _diff182(lon0, lon)

262 x = radians(lon)

263 th = self._k02n0 * x

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

265 if n0:

266 x = sth / n0

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

268 y *= nrho0 / n0

269 else:

270 x *= self._k02

271 y = _0_0

272 t = nrho0 - n0 * drho

273 x = t * x / k0

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

275

276 g = degrees360(s * th)

277 if t:

278 k0 = self._azik(t, ta)

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

280 name=name or self.name)

281

282 @Property_RO

283 def ispolar(self):

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

285 '''

286 return self._polar

287

288 isPolar = ispolar # synonym

289

290 def _k0s(self, k0):

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

292 '''

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

294 self._k02 = k02 = k0**2

295 self._k0n0 = k0 * self._n0

296 self._k02n0 = k02 * self._n0

297

298 @Property_RO

299 def lat0(self):

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

301

302 This is the latitude of minimum azimuthal scale and

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

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

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

306 projections.

307 '''

308 return self._lat0

309

310 @Property_RO

311 def lat1(self):

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

313 '''

314 return self._lat1

315

316 @Property_RO

317 def lat2(self):

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

319

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

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

322 projections.

323 '''

324 return self._lat2

325

326 @deprecated_Property_RO

327 def majoradius(self): # PYCHOK no cover

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

329 return self.equatoradius

330

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

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

333

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

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

336

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

338

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

340 '''

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

342 if self._k0 != k0:

343 _update_all(self)

344 self._k0s(k0)

345

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

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

348

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

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

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

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

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

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

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

356

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

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

359

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

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

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

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

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

365 space the nearest pole is returned.

366 '''

367 k0 = self._k0

368 n0 = self._n0

369 k0n0 = self._k0n0

370 s = self._sign

371 txi = self._txi0

372

373 x = Meter(x=x)

374 nx = k0n0 * x

375 y = Meter(y=y)

376 y_ = s * y

377 ny = k0n0 * y_

378 t = nrho0 = self._nrho0

379 y1 = nrho0 - ny

380

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

382 if den:

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

384 # dsxia = scxi0 * dsxi

385 t += drho * n0 # k0 below

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

387 t_ = txi - d_

388 d_ = (txi + t_) * d_ + _1_0

389 txi = t_ / (sqrt(d_) if d_ > EPS02 else EPS0)

390

391 ta = self._tanf(txi)

392 lat = atan1d(s * ta)

393

394 th = atan2(nx, y1)

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

396 if lon0:

397 lon += _norm180(lon0)

398 lon = _norm180(lon)

399

400 n = name or self.name

401 if LatLon is None:

402 g = degrees360(s * th)

403 if den:

404 k0 = self._azik(t, ta)

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

406 else: # PYCHOK no cover

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

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

409 return r

410

411 @Property_RO

412 def scale0(self):

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

414

415 This is the azimuthal scale on the latitude of origin

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

417 '''

418 return self._k0

419

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

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

422 '''

423 e2 = E.e2

424 e21 = E.e21

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

426 tol = _tol(_TOL0, ta0)

427 _Ta02 = Fsum(ta0).fsum2f_

428 _fabs = fabs

429 _fsum1 = fsum1f_

430 _sqrt = sqrt

431 _1, _2 = _1_0, _2_0

432 _4, _6 = _4_0, _6_0

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

434 ta02 = ta0**2

435 sca02 = ta02 + _1

436 sca0 = _sqrt(sca02)

437 sa0 = ta0 / sca0

438 sa01 = sa0 + _1

439 sa02 = sa0**2

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

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

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

443 sa021 = _1 - e2 * sa02

444

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

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

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

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

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

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

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

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

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

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

455 if _fabs(d) < tol:

456 return ta0

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

458

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

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

461 '''

462 E = self.ellipsoid

463 f1, e2 = E.f1, E.e2

464 _1 = _1_0

465

466 tb1 = f1 * ta1

467 tb2 = f1 * ta2

468 dtb12 = f1 * (tb1 + tb2)

469 scb12 = _1 + tb1**2

470 scb22 = _1 + tb2**2

471

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

473 sa12 = sa1 * sa2

474

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

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

477 esa12 = _1 + e2 * sa12

478

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

480 (ca2, sa2, ta2, scb22))

481

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

483 C = _Fsum1f_(sxi * dtb12 / dsxi, scb22, scb12).fover(scb22 * scb12 * _2_0)

484

485 S = _Fsum1f_(sa1, sa2, sa12)

486 axi *= (S * e2 + _1).fover(S + _1, raiser=False)

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

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

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

490 return ta0, C

491

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

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

494 '''

495 tol = _tol(_TOL, txi)

496 e2 = self.ellipsoid.e2

497 qx = self._qx

498

499 ta = txi

500 _Ta2 = Fsum(ta).fsum2f_

501 _fabs = fabs

502 _sqrt3 = sqrt3

503 _txif = self._txif

504 _1 = _1_0

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

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

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

508 ta2 = ta**2

509 sca2 = _1 + ta2

510 txia = _txif(ta)

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

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

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

514 if _fabs(d) < tol:

515 return ta

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

517

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

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

520

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

522

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

524 (C{str}).

525 '''

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

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

528

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

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

531

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

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

534

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

536 '''

537 k = self._k

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

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

540 (self.lat1,))

541 t = strs(t, prec=prec)

542 if k:

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

544 if self.datum != _WGS84:

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

546 if self.name:

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

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

549

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

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

552 '''

553 E = self.ellipsoid

554 _1 = _1_0

555

556 ca2 = _1_x21(ta)

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

558 sa1 = _1 + sa

559

560 es1 = sa * E.e2

561 es1m1 = sa1 * (_1 - es1)

562 es1p1 = sa1 / (_1 + es1)

563 es2m1 = _1 - sa * es1

564 es2m1a = es2m1 * E.e21 # e2m

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

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

567 t = _Fsum1f_(sa / es2m1, _atanheE(sa, E)).fover(s)

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

569

570

571 class AlbersEqualArea(_AlbersBase):

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

573

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

575 '''

576 _k = _k0_

577

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

579 '''New L{AlbersEqualArea} projection.

580

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

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

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

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

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

586

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

588 '''

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

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

591 _AlbersBase.__init__(self, *args)

592

593

594 class AlbersEqualArea2(_AlbersBase):

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

596

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

598 '''

599 _k = _k1_

600

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

602 '''New L{AlbersEqualArea2} projection.

603

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

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

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

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

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

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

610

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

612 or no convergence.

613 '''

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

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

616 _AlbersBase.__init__(self, *args)

617

618

619 class AlbersEqualArea4(_AlbersBase):

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

621 and C{cos} of both standard parallels.

622

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

624 '''

625 _k = _k1_

626

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

628 '''New L{AlbersEqualArea4} projection.

629

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

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

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

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

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

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

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

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

638

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

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

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

642 '''

643 def _Lat_s_c3(name, s, c):

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

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

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

647

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

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

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

651

652

653 class AlbersEqualAreaCylindrical(_AlbersBase):

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

655 to the cylindrical-equal-area projection.

656 '''

657 _lat1 = _lat2 = _Lat(lat1=_0_0)

658

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

660 '''New L{AlbersEqualAreaCylindrical} projection.

661

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

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 _xlat(lat, _0_0, AlbersEqualAreaCylindrical)

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

669

670

671 class AlbersEqualAreaNorth(_AlbersBase):

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

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

674 '''

675 _lat1 = _lat2 = _Lat(lat1=_90_0)

676

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

678 '''New L{AlbersEqualAreaNorth} projection.

679

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

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

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

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

684 '''

685 _xlat(lat, _90_0, AlbersEqualAreaNorth)

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

687

688

689 class AlbersEqualAreaSouth(_AlbersBase):

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

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

692 '''

693 _lat1 = _lat2 = _Lat(lat1=_N_90_0)

694

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

696 '''New L{AlbersEqualAreaSouth} projection.

697

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

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

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

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

702 '''

703 _xlat(lat, _N_90_0, AlbersEqualAreaSouth)

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

705

706

707 class Albers7Tuple(_NamedTuple):

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

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

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

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

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

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

714 the reciprocal C{1 / scale}.

715 '''

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

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

718

719

720 def _atanh1(x):

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

722 '''

723 s = fabs(x)

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

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

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

727 y, k = x, 3

728 _S2 = Fsum(y / k).fsum2f_

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

730 y *= x # x**n

731 k += 2 # 2*n + 1

732 s, d = _S2(y / k)

733 if not d:

734 break

735 elif s:

736 s = sqrt(s)

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

738 return s

739

740

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

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

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

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

745 x if f = 0.

746 '''

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

748 if e and x:

749 if E.f > 0: # .isOblate

750 x = atanh(x * e) / e

751 elif E.f < 0: # .isProlate

752 x = atan1(x * e) / e

753 return x

754

755

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

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

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

759 '''

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

761 if e:

762 d = x - y

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

764 return e

765

766

767 def _D2atanheE(x, y, E):

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

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

770 '''

771 s, e2 = _0_0, E.e2

772 if e2:

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

774 e = z = _1_0

775 k = 1

776 T = Fsum() # Taylor expansion

777 _T = T.Fsumf_

778 _C = Fsum().Fsum_

779 _S2 = Fsum().fsum2_

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

781 T *= y; P = _T(z); z *= x # PYCHOK ;

782 T *= y; Q = _T(z); z *= x # PYCHOK ;

783 e *= e2

784 k += 2

785 s, d = _S2(_C(P, Q) * e / k)

786 if not d:

787 break

788 else: # PYCHOK no cover

789 s = _1_0 - x

790 if s:

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

792 return s

793

794

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

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

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

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

799

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

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

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

803 and U{AlbersEqualArea.hpp

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

805 '''

806 # sx = x / hypot1(x)

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

808 if t > 0:

809 s = sx + sy

810 if s:

811 t = sx * sy / t

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

813 elif x != y:

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

815 return d

816

817

818 def _tol(tol, x):

819 '''(INTERNAL) Converge tolerance.

820 '''

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

822

823

824 def _1_x21(x):

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

826 '''

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

828

829

830 def _xlat(lat, f, where):

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

832 '''

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

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

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

836

837

838__all__ += _ALL_DOCS(_AlbersBase)

839

840# **) MIT License

841#

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

843#

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

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

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

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

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

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

850#

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

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

853#

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

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

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

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

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

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

860# OTHER DEALINGS IN THE SOFTWARE.