Coverage for pygeodesy/albers.py: 97%
421 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-21 13:14 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-21 13:14 -0400
2# -*- coding: utf-8 -*-
4u'''Albers Equal-Area projections.
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>}.
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
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_
38from math import atan, atan2, atanh, degrees, fabs, radians, sqrt
40__all__ = _ALL_LAZY.albers
41__version__ = '23.04.21'
43_k1_ = 'k1'
44_NUMIT = 8 # XXX 4?
45_NUMIT0 = 41 # XXX 21?
46_TERMS = 21 # XXX 16?
47_TOL0 = sqrt3(_TOL)
50def _Ks(**name_k):
51 '''(INTERNAL) Scale C{B{k} >= EPS0}.
52 '''
53 return Scalar_(Error=AlbersError, low=EPS0, **name_k) # > 0
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)
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
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
80class _AlbersBase(_NamedBase):
81 '''(INTERNAL) Base class for C{AlbersEqualArea...} projections.
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)
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
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
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)
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)
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)
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
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)
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
197 @Property_RO
198 def datum(self):
199 '''Get the datum (L{Datum}).
200 '''
201 return self._datum
203 @Property_RO
204 def ellipsoid(self):
205 '''Get the datum's ellipsoid (L{Ellipsoid}).
206 '''
207 return self.datum.ellipsoid
209 @Property_RO
210 def equatoradius(self):
211 '''Get the geodesic's equatorial radius, semi-axis (C{meter}).
212 '''
213 return self.ellipsoid.a
215 @Property_RO
216 def flattening(self):
217 '''Get the geodesic's flattening (C{float}).
218 '''
219 return self.ellipsoid.f
221 def forward(self, lat, lon, lon0=0, name=NN):
222 '''Convert a geodetic location to east- and northing.
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}).
229 @return: An L{Albers7Tuple}C{(x, y, lat, lon, gamma, scale, datum)}.
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
240 k0 = self._k0
241 n0 = self._n0
242 nrho0 = self._nrho0
243 txi0 = self._txi0
245 sa, ca = sincos2d(s * _Lat(lat=lat))
246 ta = sa / max(EPS0, ca)
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)
253 lon, _ = _diff182(lon0, lon)
254 x = radians(lon)
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
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)
275 @Property_RO
276 def ispolar(self):
277 '''Is this projection polar (C{bool})?
278 '''
279 return self._polar
281 isPolar = ispolar # synonym
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
291 @Property_RO
292 def lat0(self):
293 '''Get the latitude of the projection origin (C{degrees}).
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
303 @Property_RO
304 def lat1(self):
305 '''Get the latitude of the first parallel (C{degrees}).
306 '''
307 return self._lat1
309 @Property_RO
310 def lat2(self):
311 '''Get the latitude of the second parallel (C{degrees}).
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
319 @deprecated_Property_RO
320 def majoradius(self): # PYCHOK no cover
321 '''DEPRECATED, use property C{equatoradius}.'''
322 return self.equatoradius
324 def rescale0(self, lat, k=1): # PYCHOK no cover
325 '''Set the azimuthal scale for this projection.
327 @arg lat: Northern latitude (C{degrees}).
328 @arg k: Azimuthal scale at latitude B{C{lat}} (C{scalar}).
330 @raise AlbersError: Invalid B{C{lat}} or B{C{k}}.
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)
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.
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}.
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)}.
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
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
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)
383 ta = self._tanf(txi)
384 lat = atand(s * ta)
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)
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
403 @Property_RO
404 def scale0(self):
405 '''Get the central scale for the projection (C{float}).
407 This is the azimuthal scale on the latitude of origin
408 of the projection, see C{property lat0}.
409 '''
410 return self._k0
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
420 tb1 = f1 * ta1
421 tb2 = f1 * ta2
422 dtb12 = f1 * (tb1 + tb2)
423 scb12 = _1 + tb1**2
424 scb22 = _1 + tb2**2
426 dsn_2 = _Dsn(ta2, ta1, sa2, sa1) * _0_5
427 sa12 = sa1 * sa2
429 esa1_2 = (_1 - e2 * sa1**2) \
430 * (_1 - e2 * sa2**2)
431 esa12 = _1 + e2 * sa12
433 axi, bxi, sxi = self._a_b_sxi3((ca1, sa1, ta1, scb12),
434 (ca2, sa2, ta2, scb22))
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)
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
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
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))
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
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))
509 def toRepr(self, prec=6, **unused): # PYCHOK expected
510 '''Return a string representation of this projection.
512 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
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)
520 def toStr(self, prec=6, sep=_SPACE_, **unused): # PYCHOK expected
521 '''Return a string representation of this projection.
523 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
524 @kwarg sep: Separator to join (C{str}).
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)
541 def _txif(self, ta): # in .Ellipsoid.auxAuthalic
542 '''(INTERNAL) Function M{tan-xi from tan-phi}.
543 '''
544 E = self.ellipsoid
546 ca2 = _1_x21(ta)
547 sa = sqrt(ca2) * fabs(ta) # enforce odd parity
548 sa1 = _1_0 + sa
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
561class AlbersEqualArea(_AlbersBase):
562 '''An Albers equal-area (authalic) projection with a single standard parallel.
564 @see: L{AlbersEqualArea2} and L{AlbersEqualArea4}.
565 '''
566 _k = _k0_
568 def __init__(self, lat, k0=1, datum=_WGS84, name=NN):
569 '''New L{AlbersEqualArea} projection.
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}).
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)
584class AlbersEqualArea2(_AlbersBase):
585 '''An Albers equal-area (authalic) projection with two standard parallels.
587 @see: L{AlbersEqualArea} and L{AlbersEqualArea4}.
588 '''
589 _k = _k1_
591 def __init__(self, lat1, lat2, k1=1, datum=_WGS84, name=NN):
592 '''New L{AlbersEqualArea2} projection.
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}).
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)
609class AlbersEqualArea4(_AlbersBase):
610 '''An Albers equal-area (authalic) projection specified by the C{sin}
611 and C{cos} of both standard parallels.
613 @see: L{AlbersEqualArea} and L{AlbersEqualArea2}.
614 '''
615 _k = _k1_
617 def __init__(self, slat1, clat1, slat2, clat2, k1=1, datum=_WGS84, name=NN):
618 '''New L{AlbersEqualArea4} projection.
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}).
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)
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)
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)
649 def __init__(self, lat=_0_0, datum=_WGS84, name=NN):
650 '''New L{AlbersEqualAreaCylindrical} projection.
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)
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)
667 def __init__(self, lat=_90_0, datum=_WGS84, name=NN):
668 '''New L{AlbersEqualAreaNorth} projection.
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)
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)
685 def __init__(self, lat=_N_90_0, datum=_WGS84, name=NN):
686 '''New L{AlbersEqualAreaSouth} projection.
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)
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)
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
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
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
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
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))}.
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
811def _tol(tol, x):
812 '''(INTERNAL) Converge tolerance.
813 '''
814 return tol * max(_1_0, fabs(x))
817def _1_x21(x):
818 '''(INTERNAL) Return M{1 / (x**2 + 1)}.
819 '''
820 return _1_0 / (x**2 + _1_0)
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))
831__all__ += _ALL_DOCS(_AlbersBase)
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.