Coverage for pygeodesy/albers.py: 97%
375 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-02 08:40 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-02 08:40 -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, _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_
37from math import atan, atan2, atanh, degrees, fabs, radians, sqrt
39__all__ = _ALL_LAZY.albers
40__version__ = '23.03.19'
42_NUMIT = 8 # XXX 4?
43_NUMIT0 = 41 # XXX 21?
44_TERMS = 21 # XXX 16?
45_TOL0 = sqrt3(_TOL)
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
56def _1_x21(x):
57 '''(INTERNAL) Return M{1 / (x**2 + 1)}.
58 '''
59 return _1_0 / (x**2 + _1_0)
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))}.
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
85def _Ks(**name_k):
86 '''(INTERNAL) Scale C{B{k} >= EPS0}.
87 '''
88 return Scalar_(Error=AlbersError, low=EPS0, **name_k) # > 0
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)
98def _tol(tol, x):
99 '''(INTERNAL) Converge tolerance.
100 '''
101 return tol * max(_1_0, fabs(x))
104class _AlbersBase(_NamedBase):
105 '''(INTERNAL) Base class for C{AlbersEqualArea...} projections.
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)
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
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
144 self._qZ = qZ = _1_0 + e21 * self._atanhee(1)
145 self._qZa2 = qZ * E.a2
146 self._qx = qZ / (_2_0 * e21)
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
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)
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))
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))
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
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
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
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)
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
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
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
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
324 @Property_RO
325 def datum(self):
326 '''Get the datum (L{Datum}).
327 '''
328 return self._datum
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
336 @Property_RO
337 def flattening(self):
338 '''Get the geodesic's flattening (C{float}).
339 '''
340 return self.datum.ellipsoid.f
342 def forward(self, lat, lon, lon0=0, name=NN):
343 '''Convert a geodetic location to east- and northing.
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}).
350 @return: An L{Albers7Tuple}C{(x, y, lat, lon, gamma, scale, datum)}.
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
361 k0 = self._k0
362 n0 = self._n0
363 nrho0 = self._nrho0
364 txi0 = self._txi0
366 sa, ca = sincos2d(_Lat(lat=lat) * s)
367 ca = max(EPS0, ca)
368 ta = sa / ca
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)
374 lon, _ = _diff182(lon0, lon)
375 b = radians(lon)
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
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)
396 @Property_RO
397 def ispolar(self):
398 '''Is this projection polar (C{bool})?
399 '''
400 return self._polar
402 isPolar = ispolar # synonym
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
412 @Property_RO
413 def lat0(self):
414 '''Get the latitude of the projection origin (C{degrees}).
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
424 @Property_RO
425 def lat1(self):
426 '''Get the latitude of the first parallel (C{degrees}).
427 '''
428 return self._lat1
430 @Property_RO
431 def lat2(self):
432 '''Get the latitude of the second parallel (C{degrees}).
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
440 @deprecated_Property_RO
441 def majoradius(self): # PYCHOK no cover
442 '''DEPRECATED, use property C{equatoradius}.'''
443 return self.equatoradius
445 def rescale0(self, lat, k=1): # PYCHOK no cover
446 '''Set the azimuthal scale for this projection.
448 @arg lat: Northern latitude (C{degrees}).
449 @arg k: Azimuthal scale at latitude B{C{lat}} (C{scalar}).
451 @raise AlbersError: Invalid B{C{lat}} or B{C{k}}.
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)
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.
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}.
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)}.
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)
484 k0 = self._k0
485 n0 = self._n0
486 k0n0 = self._k0n0
487 nrho0 = self._nrho0
488 txi0 = self._txi0
490 y_ = self._sign * y
491 nx = k0n0 * x
492 ny = k0n0 * y_
493 y1 = nrho0 - ny
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)
503 ta = self._tanf(txi)
504 lat = atand(ta * self._sign)
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)
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
523 @Property_RO
524 def scale0(self):
525 '''Get the central scale for the projection (C{float}).
527 This is the azimuthal scale on the latitude of origin
528 of the projection, see C{property lat0}.
529 '''
530 return self._k0
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
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
545 esa1_2 = (_1_0 - e2 * sa1**2) \
546 * (_1_0 - e2 * sa2**2)
547 esa12 = _1_0 + e2 * sa1 * sa2
549 dsn = _Dsn(ta2, ta1, sa2, sa1)
550 axi, bxi, sxi = self._a_b_sxi3((ca1, sa1, ta1, scb12),
551 (ca2, sa2, ta2, scb22))
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)
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
562 def _tanf(self, txi): # called by .Ellipsoid.auxAuthalic
563 '''(INTERNAL) Function M{tan-phi from tan-xi}.
564 '''
565 tol = _tol(_TOL, txi)
567 e2 = self.datum.ellipsoid.e2
568 qx = self._qx
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))
585 def _txif(self, ta): # in .Ellipsoid.auxAuthalic
586 '''(INTERNAL) Function M{tan-xi from tan-phi}.
587 '''
588 E = self.datum.ellipsoid
590 ca2 = _1_x21(ta)
591 sa = sqrt(ca2) * fabs(ta) # enforce odd parity
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
605class AlbersEqualArea(_AlbersBase):
606 '''An Albers equal-area (authalic) projection with a single standard parallel.
608 @see: L{AlbersEqualArea2} and L{AlbersEqualArea4}.
609 '''
610 def __init__(self, lat, k0=1, datum=_WGS84, name=NN):
611 '''New L{AlbersEqualArea} projection.
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}).
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)
626class AlbersEqualArea2(_AlbersBase):
627 '''An Albers equal-area (authalic) projection with two standard parallels.
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.
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}).
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)
649class AlbersEqualArea4(_AlbersBase):
650 '''An Albers equal-area (authalic) projection specified by the C{sin}
651 and C{cos} of both standard parallels.
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.
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}).
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
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)
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)
687 def __init__(self, datum=_WGS84, name=NN):
688 '''New L{AlbersEqualAreaCylindrical} projection.
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)
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)
703 def __init__(self, datum=_WGS84, name=NN):
704 '''New L{AlbersEqualAreaNorth} projection.
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)
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)
719 def __init__(self, datum=_WGS84, name=NN):
720 '''New L{AlbersEqualAreaSouth} projection.
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)
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)
742__all__ += _ALL_DOCS(_AlbersBase)
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.