Coverage for pygeodesy/lcc.py: 96%
281 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-28 15:52 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-28 15:52 -0400
2# -*- coding: utf-8 -*-
4u'''Lambert Conformal Conic (LCC) projection.
6Lambert conformal conic projection for 1- or 2-Standard Parallels classes L{Conic}, L{Conics} registry, L{LCCError}
7and position class L{Lcc}.
9See U{LCC<https://WikiPedia.org/wiki/Lambert_conformal_conic_projection>}, U{Lambert
10Conformal Conic to Geographic Transformation Formulae
11<https://www.Linz.govt.NZ/data/geodetic-system/coordinate-conversion/projection-conversions/lambert-conformal-conic-geographic>},
12U{Lambert Conformal Conic Projection<https://MathWorld.Wolfram.com/LambertConformalConicProjection.html>}
13and John P. Snyder U{'Map Projections - A Working Manual'<https://Pubs.USGS.gov/pp/1395/report.pdf>}, 1987, pp 107-109.
15@var Conics.Be08Lb: Conic(name='Be08Lb', lat0=50.797815, lon0=4.35921583, par1=49.8333339, par2=51.1666672, E0=649328, N0=665262, k0=1, SP=2, datum=Datum(name='GRS80', ellipsoid=Ellipsoids.GRS80, transform=Transforms.WGS84),
16@var Conics.Be72Lb: Conic(name='Be72Lb', lat0=90, lon0=4.3674867, par1=49.8333339, par2=51.1666672, E0=150000.013, N0=5400088.438, k0=1, SP=2, datum=Datum(name='NAD83', ellipsoid=Ellipsoids.GRS80, transform=Transforms.NAD83),
17@var Conics.Fr93Lb: Conic(name='Fr93Lb', lat0=46.5, lon0=3, par1=49, par2=44, E0=700000, N0=6600000, k0=1, SP=2, datum=Datum(name='WGS84', ellipsoid=Ellipsoids.WGS84, transform=Transforms.WGS84),
18@var Conics.MaNLb: Conic(name='MaNLb', lat0=33.3, lon0=-5.4, par1=31.73, par2=34.87, E0=500000, N0=300000, k0=1, SP=2, datum=Datum(name='NTF', ellipsoid=Ellipsoids.Clarke1880IGN, transform=Transforms.NTF),
19@var Conics.MxLb: Conic(name='MxLb', lat0=12, lon0=-102, par1=17.5, par2=29.5, E0=2500000, N0=0, k0=1, SP=2, datum=Datum(name='WGS84', ellipsoid=Ellipsoids.WGS84, transform=Transforms.WGS84),
20@var Conics.PyT_Lb: Conic(name='PyT_Lb', lat0=46.8, lon0=2.33722917, par1=45.8989389, par2=47.6960144, E0=600000, N0=200000, k0=1, SP=2, datum=Datum(name='NTF', ellipsoid=Ellipsoids.Clarke1880IGN, transform=Transforms.NTF),
21@var Conics.USA_Lb: Conic(name='USA_Lb', lat0=23, lon0=-96, par1=33, par2=45, E0=0, N0=0, k0=1, SP=2, datum=Datum(name='WGS84', ellipsoid=Ellipsoids.WGS84, transform=Transforms.WGS84),
22@var Conics.WRF_Lb: Conic(name='WRF_Lb', lat0=40, lon0=-97, par1=33, par2=45, E0=0, N0=0, k0=1, SP=2, datum=Datum(name='WGS84', ellipsoid=Ellipsoids.WGS84, transform=Transforms.WGS84)
23'''
24# make sure int/int division yields float quotient, see .basics
25from __future__ import division as _; del _ # PYCHOK semicolon
27from pygeodesy.basics import copysign0, _xinstanceof, _xsubclassof
28from pygeodesy.constants import EPS, EPS02, PI_2, _float as _F, _0_0, _0_5, \
29 _1_0, _2_0, _90_0
30from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB
31from pygeodesy.datums import Datums, _ellipsoidal_datum
32from pygeodesy.errors import _IsnotError, _ValueError
33from pygeodesy.fmath import _ALL_LAZY, hypot
34from pygeodesy.interns import NN, _COMMASPACE_, _ellipsoidal_, _GRS80_, _k0_, \
35 _lat0_, _lon0_, _m_, _NAD83_, _NTF_, _SPACE_, _WGS84_, \
36 _C_ # PYCHOK used!
37# from pygeodesy.lazily import _ALL_LAZY # from .fmath
38from pygeodesy.named import _lazyNamedEnumItem as _lazy, _NamedBase, \
39 _NamedEnum, _NamedEnumItem, nameof, _xnamed
40from pygeodesy.namedTuples import EasNor3Tuple, LatLonDatum3Tuple, \
41 LatLon2Tuple, _LL4Tuple, PhiLam2Tuple
42from pygeodesy.props import deprecated_method, Property, Property_RO, \
43 _update_all
44from pygeodesy.streprs import Fmt, _fstrENH2, _xzipairs
45from pygeodesy.units import Easting, Height, _heigHt, Lam_, Northing, \
46 Phi_, Scalar_
47from pygeodesy.utily import degrees90, degrees180, sincos2, tanPI_2_2
49from math import atan, fabs, log, radians, sin, sqrt
51__all__ = _ALL_LAZY.lcc
52__version__ = '23.05.26'
54_E0_ = 'E0'
55_N0_ = 'N0'
56_par1_ = 'par1'
57_par2_ = 'par2'
58_SP_ = 'SP'
61class Conic(_NamedEnumItem):
62 '''Lambert conformal conic projection (1- or 2-SP).
63 '''
64 _auth = NN # authorization (C{str})
65 _datum = None # datum (L{Datum})
66 _name = NN # Conic.__name__, set below
68 _e = _0_0 # ellipsoid excentricity (C{float})
69 _E0 = _0_0 # false easting (C{float})
70 _k0 = _1_0 # scale factor (C{float})
71 _N0 = _0_0 # false northing (C{float})
72 _SP = 0 # 1- or 2-SP (C{int})
74 _opt3 = _0_0 # optional, longitude (C{radians})
75 _par1 = _0_0 # 1st std parallel (C{radians})
76 _par2 = _0_0 # 2nd std parallel (C{radians})
77 _phi0 = _0_0 # origin lat (C{radians})
78 _lam0 = _0_0 # origin lon (C{radians})
80 _aF = _0_0 # precomputed F (C{float})
81 _n = _0_0 # precomputed n (C{float})
82 _1_n = _0_0 # precomputed 1 / n (C{float})
83 _r0 = _0_0 # precomputed rho0 (C{float})
85 def __init__(self, latlon0, par1, par2=None, E0=0, N0=0,
86 k0=1, opt3=0, name=NN, auth=NN):
87 '''New Lambert conformal conic projection.
89 @arg latlon0: Origin with (ellipsoidal) datum (C{LatLon}).
90 @arg par1: First standard parallel (C{degrees90}).
91 @kwarg par2: Optional, second standard parallel (C{degrees90}).
92 @kwarg E0: Optional, false easting (C{meter}).
93 @kwarg N0: Optional, false northing (C{meter}).
94 @kwarg k0: Optional scale factor (C{scalar}).
95 @kwarg opt3: Optional meridian (C{degrees180}).
96 @kwarg name: Optional name of the conic (C{str}).
97 @kwarg auth: Optional authentication authority (C{str}).
99 @return: A Lambert projection (L{Conic}).
101 @raise TypeError: Non-ellipsoidal B{C{latlon0}}.
103 @raise ValueError: Invalid B{C{par1}}, B{C{par2}},
104 B{C{E0}}, B{C{N0}}, B{C{k0}}
105 or B{C{opt3}}.
107 @example:
109 >>> from pygeodesy import Conic, Datums, ellipsoidalNvector
110 >>> ll0 = ellipsoidalNvector.LatLon(23, -96, datum=Datums.NAD27)
111 >>> Snyder = Conic(ll0, 33, 45, E0=0, N0=0, name='Snyder')
112 '''
113 if latlon0 is not None:
114 _xinstanceof(_LLEB, latlon0=latlon0)
115 self._phi0, self._lam0 = latlon0.philam
117 self._par1 = Phi_(par1=par1)
118 self._par2 = self._par1 if par2 is None else Phi_(par2=par2)
120 if k0 != 1:
121 self._k0 = Scalar_(k0=k0)
122 if E0:
123 self._E0 = Northing(E0=E0, falsed=True)
124 if N0:
125 self._N0 = Easting(N0=N0, falsed=True)
126 if opt3:
127 self._opt3 = Lam_(opt3=opt3)
129 self.toDatum(latlon0.datum)._dup2(self)
130 self._register(Conics, name)
131 elif name:
132 self.name = name
133 if auth:
134 self._auth = str(auth)
136 @Property_RO
137 def auth(self):
138 '''Get the authentication authority (C{str}).
139 '''
140 return self._auth
142 @deprecated_method
143 def convertDatum(self, datum):
144 '''DEPRECATED, use method L{Conic.toDatum}.'''
145 return self.toDatum(datum)
147 @Property_RO
148 def datum(self):
149 '''Get the datum (L{Datum}).
150 '''
151 return self._datum
153 @Property_RO
154 def E0(self):
155 '''Get the false easting (C{meter}).
156 '''
157 return self._E0
159 @Property_RO
160 def k0(self):
161 '''Get scale factor (C{float}).
162 '''
163 return self._k0
165 @Property_RO
166 def lat0(self):
167 '''Get the origin latitude (C{degrees90}).
168 '''
169 return degrees90(self._phi0)
171 @Property_RO
172 def latlon0(self):
173 '''Get the central origin (L{LatLon2Tuple}C{(lat, lon)}).
174 '''
175 return LatLon2Tuple(self.lat0, self.lon0, name=self.name)
177 @Property_RO
178 def lam0(self):
179 '''Get the central meridian (C{radians}).
180 '''
181 return self._lam0
183 @Property_RO
184 def lon0(self):
185 '''Get the central meridian (C{degrees180}).
186 '''
187 return degrees180(self._lam0)
189 @Property_RO
190 def N0(self):
191 '''Get the false northing (C{meter}).
192 '''
193 return self._N0
195 @Property_RO
196 def name2(self):
197 '''Get the conic and datum names as "conic.datum" (C{str}).
198 '''
199 return self._DOT_(self.datum.name)
201 @Property_RO
202 def opt3(self):
203 '''Get the optional meridian (C{degrees180}).
204 '''
205 return degrees180(self._opt3)
207 @Property_RO
208 def par1(self):
209 '''Get the 1st standard parallel (C{degrees90}).
210 '''
211 return degrees90(self._par1)
213 @Property_RO
214 def par2(self):
215 '''Get the 2nd standard parallel (C{degrees90}).
216 '''
217 return degrees90(self._par2)
219 @Property_RO
220 def phi0(self):
221 '''Get the origin latitude (C{radians}).
222 '''
223 return self._phi0
225 @Property_RO
226 def philam0(self):
227 '''Get the central origin (L{PhiLam2Tuple}C{(phi, lam)}).
228 '''
229 return PhiLam2Tuple(self.phi0, self.lam0, name=self.name)
231 @Property_RO
232 def SP(self):
233 '''Get the number of standard parallels (C{int}).
234 '''
235 return self._SP
237 def toDatum(self, datum):
238 '''Convert this conic to the given datum.
240 @arg datum: Ellipsoidal datum to use (L{Datum}, L{Ellipsoid},
241 L{Ellipsoid2} or L{a_f2Tuple}).
243 @return: Converted conic, unregistered (L{Conic}).
245 @raise TypeError: Non-ellipsoidal B{C{datum}}.
246 '''
247 d = _ellipsoidal_datum(datum, name=self.name)
248 E = d.ellipsoid
249 if not E.isEllipsoidal:
250 raise _IsnotError(_ellipsoidal_, datum=datum)
252 c = self
253 if c._e != E.e or c._datum != d:
255 c = Conic(None, 0, name=self._name)
256 self._dup2(c)
257 c._datum = d
258 c._e = E.e
260 if fabs(c._par1 - c._par2) < EPS:
261 m1 = c._mdef(c._phi0)
262 t1 = c._tdef(c._phi0)
263 t0 = t1
264 k = 1 # _1_0
265 n = sin(c._phi0)
266 sp = 1
267 else:
268 m1 = c._mdef(c._par1)
269 m2 = c._mdef(c._par2)
270 t1 = c._tdef(c._par1)
271 t2 = c._tdef(c._par2)
272 t0 = c._tdef(c._phi0)
273 k = c._k0
274 n = (log(m1) - log(m2)) \
275 / (log(t1) - log(t2))
276 sp = 2
278 F = m1 / (n * pow(t1, n))
280 c._aF = k * E.a * F
281 c._n = n
282 c._1_n = _1_0 / n
283 c._r0 = c._rdef(t0)
284 c._SP = sp
286 return c
288 def toStr(self, prec=8, name=NN, **unused): # PYCHOK expected
289 '''Return this conic as a string.
291 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
292 @kwarg name: Override name (C{str}) or C{None} to exclude
293 this conic's name.
295 @return: Conic attributes (C{str}).
296 '''
297 a = [name, prec, _lat0_, _lon0_, _par1_, _par2_,
298 _E0_, _N0_, _k0_, _SP_]
299 if self._SP == 1:
300 _ = a.pop(a.index(_par2_))
301 return self._instr(datum=self.datum, *a)
303 def _dup2(self, c):
304 '''(INTERNAL) Copy this conic to C{c}.
306 @arg c: Duplicate (L{Conic}).
307 '''
308 _update_all(c)
310 c._auth = self._auth
311 c._datum = self._datum
313 c._e = self._e
314 c._E0 = self._E0
315 c._k0 = self._k0
316 c._N0 = self._N0
317 c._SP = self._SP
319 c._par1 = self._par1
320 c._par2 = self._par2
321 c._phi0 = self._phi0
322 c._lam0 = self._lam0
323 c._opt3 = self._opt3
325 c._aF = self._aF
326 c._n = self._n
327 c._1_n = self._1_n
328 c._r0 = self._r0
330 def _mdef(self, a):
331 '''(INTERNAL) Compute m(a).
332 '''
333 s, c = sincos2(a)
334 s = _1_0 - (s * self._e)**2
335 return (c / sqrt(s)) if s > EPS02 else _0_0
337 def _pdef(self, a):
338 '''(INTERNAL) Compute p(a).
339 '''
340 s = self._e * sin(a)
341 return pow((_1_0 - s) / (_1_0 + s), self._e * _0_5)
343 def _rdef(self, t):
344 '''(INTERNAL) Compute r(t).
345 '''
346 return self._aF * pow(t, self._n)
348 def _tdef(self, a):
349 '''(INTERNAL) Compute t(lat).
350 '''
351 return max(_0_0, tanPI_2_2(-a) / self._pdef(a))
353 def _xdef(self, t_x):
354 '''(INTERNAL) Compute x(t_x).
355 '''
356 return PI_2 - _2_0 * atan(t_x) # XXX + self._phi0
359Conic._name = Conic.__name__
362class Conics(_NamedEnum):
363 '''(INTERNAL) L{Conic} registry, I{must} be a sub-class
364 to accommodate the L{_LazyNamedEnumItem} properties.
365 '''
366 def _Lazy(self, lat, lon, datum_name, *args, **kwds):
367 '''(INTERNAL) Instantiate the L{Conic}.
368 '''
369 return Conic(_LLEB(lat, lon, datum=Datums.get(datum_name)), *args, **kwds)
371Conics = Conics(Conic) # PYCHOK singleton
372'''Some pre-defined L{Conic}s, all I{lazily} instantiated.'''
373Conics._assert( # <https://SpatialReference.org/ref/sr-org/...>
374# AsLb = _lazy('AsLb', _F(-14.2666667), _F(170), _NAD27_, _0_0, _0_0,
375# E0=_F(500000), N0=_0_0, auth='EPSG:2155'), # American Samoa ... SP=1 !
376 Be08Lb = _lazy('Be08Lb', _F(50.7978150), _F(4.359215833), _GRS80_, _F(49.8333339), _F(51.1666672),
377 E0=_F(649328.0), N0=_F(665262.0), auth='EPSG:9802'), # Belgium
378 Be72Lb = _lazy('Be72Lb', _90_0, _F(4.3674867), _NAD83_, _F(49.8333339), _F(51.1666672),
379 E0=_F(150000.013), N0=_F(5400088.438), auth='EPSG:31370'), # Belgium
380 Fr93Lb = _lazy('Fr93Lb', _F(46.5), _F(3), _WGS84_, _F(49), _F(44),
381 E0=_F(700000), N0=_F(6600000), auth='EPSG:2154'), # RFG93, France
382 MaNLb = _lazy('MaNLb', _F(33.3), _F(-5.4), _NTF_, _F(31.73), _F(34.87),
383 E0=_F(500000), N0=_F(300000)), # Marocco
384 MxLb = _lazy('MxLb', _F(12), _F(-102), _WGS84_, _F(17.5), _F(29.5),
385 E0=_F(2500000), N0=_0_0, auth='EPSG:2155'), # Mexico
386 PyT_Lb = _lazy('PyT_Lb', _F(46.8), _F(2.33722917), _NTF_, _F(45.89893890000052), _F(47.69601440000037),
387 E0=_F(600000), N0=_F(200000), auth='Test'), # France?
388 USA_Lb = _lazy('USA_Lb', _F(23), _F(-96), _WGS84_, _F(33), _F(45),
389 E0=_0_0, N0=_0_0), # Conterminous, contiguous USA?
390 WRF_Lb = _lazy('WRF_Lb', _F(40), _F(-97), _WGS84_, _F(33), _F(45),
391 E0=_0_0, N0=_0_0, auth='EPSG:4326') # World
392)
395class LCCError(_ValueError):
396 '''Lambert Conformal Conic C{LCC} or other L{Lcc} issue.
397 '''
398 pass
401class Lcc(_NamedBase):
402 '''Lambert conformal conic East-/Northing location.
403 '''
404 _conic = Conics.WRF_Lb # Lambert projection (L{Conic})
405 _easting = _0_0 # Easting (C{float})
406 _height = 0 # height (C{meter})
407 _northing = _0_0 # Northing (C{float})
409 def __init__(self, e, n, h=0, conic=Conics.WRF_Lb, name=NN):
410 '''New L{Lcc} Lamber conformal conic position.
412 @arg e: Easting (C{meter}).
413 @arg n: Northing (C{meter}).
414 @kwarg h: Optional height (C{meter}).
415 @kwarg conic: Optional, the conic projection (L{Conic}).
416 @kwarg name: Optional name (C{str}).
418 @return: The Lambert location (L{Lcc}).
420 @raise LCCError: Invalid B{C{h}} or invalid or
421 negative B{C{e}} or B{C{n}}.
423 @raise TypeError: If B{C{conic}} is not L{Conic}.
425 @example:
427 >>> lb = Lcc(448251, 5411932.0001)
428 '''
429 if conic not in (None, Lcc._conic):
430 self.conic = conic
431 self._easting = Easting(e, falsed=conic.E0 > 0, Error=LCCError)
432 self._northing = Northing(n, falsed=conic.N0 > 0, Error=LCCError)
433 if h:
434 self._height = Height(h=h, Error=LCCError)
435 if name:
436 self.name = name
438 @Property
439 def conic(self):
440 '''Get the conic projection (L{Conic}).
441 '''
442 return self._conic
444 @conic.setter # PYCHOK setter!
445 def conic(self, conic):
446 '''Set the conic projection (L{Conic}).
448 @raise TypeError: Invalid B{C{conic}}.
449 '''
450 _xinstanceof(Conic, conic=conic)
451 if conic != self._conic:
452 _update_all(self)
453 self._conic = conic
455# def dup(self, name=NN, **e_n_h_conic): # PYCHOK signature
456# '''Duplicate this location with some attributes modified.
457#
458# @kwarg e_n_h_conic: Use keyword argument C{B{e}=...}, C{B{n}=...},
459# C{B{h}=...} and/or C{B{conic}=...} to override
460# the current C{easting}, C{northing} C{height}
461# or C{conic} projection, respectively.
462# '''
463# def _args_kwds(e=None, n=None, **kwds):
464# return (e, n), kwds
465#
466# kwds = _xkwds(e_n_h_conic, e=self.easting, n=self.northing,
467# h=self.height, conic=self.conic,
468# name=name or self.name)
469# args, kwds = _args_kwds(**kwds)
470# return self.__class__(*args, **kwds) # .classof
472 @Property_RO
473 def easting(self):
474 '''Get the easting (C{meter}).
475 '''
476 return self._easting
478 @Property_RO
479 def height(self):
480 '''Get the height (C{meter}).
481 '''
482 return self._height
484 @Property_RO
485 def latlon(self):
486 '''Get the lat- and longitude in C{degrees} (L{LatLon2Tuple}).
487 '''
488 ll = self.toLatLon(LatLon=None, datum=None)
489 return LatLon2Tuple(ll.lat, ll.lon, name=self.name)
491 @Property_RO
492 def latlonheight(self):
493 '''Get the lat-, longitude and height (L{LatLon3Tuple}C{(lat, lon, height)}).
494 '''
495 return self.latlon.to3Tuple(self.height)
497 @Property_RO
498 def latlonheightdatum(self):
499 '''Get the lat-, longitude in C{degrees} with height and datum (L{LatLon4Tuple}C{(lat, lon, height, datum)}).
500 '''
501 return self.latlonheight.to4Tuple(self.conic.datum)
503 @Property_RO
504 def northing(self):
505 '''Get the northing (C{meter}).
506 '''
507 return self._northing
509 @Property_RO
510 def philam(self):
511 '''Get the lat- and longitude in C{radians} (L{PhiLam2Tuple}).
512 '''
513 return PhiLam2Tuple(radians(self.latlon.lat),
514 radians(self.latlon.lon), name=self.name)
516 @Property_RO
517 def philamheight(self):
518 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).
519 '''
520 return self.philam.to3Tuple(self.height)
522 @Property_RO
523 def philamheightdatum(self):
524 '''Get the lat-, longitude in C{radians} with height and datum (L{PhiLam4Tuple}C{(phi, lam, height, datum)}).
525 '''
526 return self.philamheight.to4Tuple(self.datum)
528 @deprecated_method
529 def to3lld(self, datum=None): # PYCHOK no cover
530 '''DEPRECATED, use method C{toLatLon}.
532 @kwarg datum: Optional datum to use, otherwise use this
533 B{C{Lcc}}'s conic.datum (C{Datum}).
535 @return: A L{LatLonDatum3Tuple}C{(lat, lon, datum)}.
537 @raise TypeError: If B{C{datum}} is not ellipsoidal.
538 '''
539 if datum in (None, self.conic.datum):
540 r = LatLonDatum3Tuple(self.latlon.lat,
541 self.latlon.lon,
542 self.conic.datum, name=self.name)
543 else:
544 r = self.toLatLon(LatLon=None, datum=datum)
545 r = LatLonDatum3Tuple(r.lat, r.lon, r.datum, name=r.name)
546 return r
548 def toLatLon(self, LatLon=None, datum=None, height=None, **LatLon_kwds):
549 '''Convert this L{Lcc} to an (ellipsoidal) geodetic point.
551 @kwarg LatLon: Optional, ellipsoidal class to return the
552 geodetic point (C{LatLon}) or C{None}.
553 @kwarg datum: Optional datum to use, otherwise use this
554 B{C{Lcc}}'s conic.datum (L{Datum}, L{Ellipsoid},
555 L{Ellipsoid2} or L{a_f2Tuple}).
556 @kwarg height: Optional height for the point, overriding
557 the default height (C{meter}).
558 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
559 arguments, ignored if C{B{LatLon} is None}.
561 @return: The point (B{C{LatLon}}) or a
562 L{LatLon4Tuple}C{(lat, lon, height, datum)}
563 if B{C{LatLon}} is C{None}.
565 @raise TypeError: If B{C{LatLon}} or B{C{datum}} is
566 not ellipsoidal or not valid.
567 '''
568 if LatLon:
569 _xsubclassof(_LLEB, LatLon=LatLon)
571 c = self.conic
572 if datum not in (None, self.conic.datum):
573 c = c.toDatum(datum)
575 e = self.easting - c._E0
576 n = c._r0 - self.northing + c._N0
578 r_ = copysign0(hypot(e, n), c._n)
579 t_ = pow(r_ / c._aF, c._1_n)
581 x = c._xdef(t_) # XXX c._lam0
582 for self._iteration in range(10): # max 4 trips
583 p, x = x, c._xdef(t_ * c._pdef(x))
584 if fabs(x - p) < 1e-9: # XXX EPS too small?
585 break
586 lat = degrees90(x)
587 lon = degrees180((atan(e / n) + c._opt3) * c._1_n + c._lam0)
589 h = _heigHt(self, height)
590 return _LL4Tuple(lat, lon, h, c.datum, LatLon, LatLon_kwds,
591 inst=self, name=self.name)
593 def toRepr(self, prec=0, fmt=Fmt.SQUARE, sep=_COMMASPACE_, m=_m_, C=False, **unused): # PYCHOK expected
594 '''Return a string representation of this L{Lcc} position.
596 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
597 @kwarg fmt: Enclosing backets format (C{str}).
598 @kwarg sep: Optional separator between name:values (C{str}).
599 @kwarg m: Optional unit of the height, default meter (C{str}).
600 @kwarg C: Optionally, include name of conic and datum (C{bool}).
602 @return: This Lcc as "[E:meter, N:meter, H:m, C:Conic.Datum]"
603 (C{str}).
604 '''
605 t, T = _fstrENH2(self, prec, m)
606 if C:
607 t += self.conic.name2,
608 T += _C_,
609 return _xzipairs(T, t, sep=sep, fmt=fmt)
611 def toStr(self, prec=0, sep=_SPACE_, m=_m_): # PYCHOK expected
612 '''Return a string representation of this L{Lcc} position.
614 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
615 @kwarg sep: Optional separator to join (C{str}) or C{None}
616 to return an unjoined C{tuple} of C{str}s.
617 @kwarg m: Optional height units, default C{meter} (C{str}).
619 @return: This Lcc as I{"easting nothing"} in C{meter} plus
620 I{" height"} suffixed with B{C{m}} if height is
621 non-zero (C{str}).
622 '''
623 t, _ = _fstrENH2(self, prec, m)
624 return t if sep is None else sep.join(t)
627def toLcc(latlon, conic=Conics.WRF_Lb, height=None, Lcc=Lcc, name=NN,
628 **Lcc_kwds):
629 '''Convert an (ellipsoidal) geodetic point to a I{Lambert} location.
631 @arg latlon: Ellipsoidal point (C{LatLon}).
632 @kwarg conic: Optional Lambert projection to use (L{Conic}).
633 @kwarg height: Optional height for the point, overriding the
634 default height (C{meter}).
635 @kwarg Lcc: Optional class to return the I{Lambert} location
636 (L{Lcc}).
637 @kwarg name: Optional B{C{Lcc}} name (C{str}).
638 @kwarg Lcc_kwds: Optional, additional B{C{Lcc}} keyword
639 arguments, ignored if B{C{Lcc}} is C{None}.
641 @return: The I{Lambert} location (L{Lcc}) or an
642 L{EasNor3Tuple}C{(easting, northing, height)}
643 if C{B{Lcc} is None}.
645 @raise TypeError: If B{C{latlon}} is not ellipsoidal.
646 '''
647 _xinstanceof(_LLEB, latlon=latlon)
649 a, b = latlon.philam
650 c = conic.toDatum(latlon.datum)
652 t = c._n * (b - c._lam0) - c._opt3
653 st, ct = sincos2(t)
655 r = c._rdef(c._tdef(a))
656 e = c._E0 + r * st
657 n = c._N0 + c._r0 - r * ct
659 h = _heigHt(latlon, height)
660 r = EasNor3Tuple(e, n, h) if Lcc is None else \
661 Lcc(e, n, h=h, conic=c, **Lcc_kwds)
662 return _xnamed(r, name or nameof(latlon))
665if __name__ == '__main__':
667 from pygeodesy.interns import _NL_, _NLATvar_
668 from pygeodesy.lazily import printf
670 # __doc__ of this file, force all into registery
671 t = _NL_ + Conics.toRepr(all=True, asorted=True)
672 printf(_NLATvar_.join(t.split(_NL_)))
674# **) MIT License
675#
676# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
677#
678# Permission is hereby granted, free of charge, to any person obtaining a
679# copy of this software and associated documentation files (the "Software"),
680# to deal in the Software without restriction, including without limitation
681# the rights to use, copy, modify, merge, publish, distribute, sublicense,
682# and/or sell copies of the Software, and to permit persons to whom the
683# Software is furnished to do so, subject to the following conditions:
684#
685# The above copyright notice and this permission notice shall be included
686# in all copies or substantial portions of the Software.
687#
688# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
689# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
690# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
691# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
692# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
693# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
694# OTHER DEALINGS IN THE SOFTWARE.