Coverage for pygeodesy/lcc.py: 96%

280 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-03-31 10:52 -0400

1 

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

3 

4u'''Lambert Conformal Conic (LCC) projection. 

5 

6Lambert conformal conic projection for 1- or 2-Standard Parallels classes L{Conic}, L{Conics} registry, L{LCCError} 

7and position class L{Lcc}. 

8 

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. 

14 

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 

26 

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, Lam_, Northing, Phi_, Scalar_ 

46from pygeodesy.utily import degrees90, degrees180, sincos2, tanPI_2_2 

47 

48from math import atan, fabs, log, radians, sin, sqrt 

49 

50__all__ = _ALL_LAZY.lcc 

51__version__ = '23.03.19' 

52 

53_E0_ = 'E0' 

54_N0_ = 'N0' 

55_par1_ = 'par1' 

56_par2_ = 'par2' 

57_SP_ = 'SP' 

58 

59 

60class Conic(_NamedEnumItem): 

61 '''Lambert conformal conic projection (1- or 2-SP). 

62 ''' 

63 _auth = NN # authorization (C{str}) 

64 _datum = None # datum (L{Datum}) 

65 _name = NN # Conic.__name__, set below 

66 

67 _e = _0_0 # ellipsoid excentricity (C{float}) 

68 _E0 = _0_0 # false easting (C{float}) 

69 _k0 = _1_0 # scale factor (C{float}) 

70 _N0 = _0_0 # false northing (C{float}) 

71 _SP = 0 # 1- or 2-SP (C{int}) 

72 

73 _opt3 = _0_0 # optional, longitude (C{radians}) 

74 _par1 = _0_0 # 1st std parallel (C{radians}) 

75 _par2 = _0_0 # 2nd std parallel (C{radians}) 

76 _phi0 = _0_0 # origin lat (C{radians}) 

77 _lam0 = _0_0 # origin lon (C{radians}) 

78 

79 _aF = _0_0 # precomputed F (C{float}) 

80 _n = _0_0 # precomputed n (C{float}) 

81 _1_n = _0_0 # precomputed 1 / n (C{float}) 

82 _r0 = _0_0 # precomputed rho0 (C{float}) 

83 

84 def __init__(self, latlon0, par1, par2=None, E0=0, N0=0, 

85 k0=1, opt3=0, name=NN, auth=NN): 

86 '''New Lambert conformal conic projection. 

87 

88 @arg latlon0: Origin with (ellipsoidal) datum (C{LatLon}). 

89 @arg par1: First standard parallel (C{degrees90}). 

90 @kwarg par2: Optional, second standard parallel (C{degrees90}). 

91 @kwarg E0: Optional, false easting (C{meter}). 

92 @kwarg N0: Optional, false northing (C{meter}). 

93 @kwarg k0: Optional scale factor (C{scalar}). 

94 @kwarg opt3: Optional meridian (C{degrees180}). 

95 @kwarg name: Optional name of the conic (C{str}). 

96 @kwarg auth: Optional authentication authority (C{str}). 

97 

98 @return: A Lambert projection (L{Conic}). 

99 

100 @raise TypeError: Non-ellipsoidal B{C{latlon0}}. 

101 

102 @raise ValueError: Invalid B{C{par1}}, B{C{par2}}, 

103 B{C{E0}}, B{C{N0}}, B{C{k0}} 

104 or B{C{opt3}}. 

105 

106 @example: 

107 

108 >>> from pygeodesy import Conic, Datums, ellipsoidalNvector 

109 >>> ll0 = ellipsoidalNvector.LatLon(23, -96, datum=Datums.NAD27) 

110 >>> Snyder = Conic(ll0, 33, 45, E0=0, N0=0, name='Snyder') 

111 ''' 

112 if latlon0 is not None: 

113 _xinstanceof(_LLEB, latlon0=latlon0) 

114 self._phi0, self._lam0 = latlon0.philam 

115 

116 self._par1 = Phi_(par1=par1) 

117 self._par2 = self._par1 if par2 is None else Phi_(par2=par2) 

118 

119 if k0 != 1: 

120 self._k0 = Scalar_(k0=k0) 

121 if E0: 

122 self._E0 = Northing(E0=E0, falsed=True) 

123 if N0: 

124 self._N0 = Easting(N0=N0, falsed=True) 

125 if opt3: 

126 self._opt3 = Lam_(opt3=opt3) 

127 

128 self.toDatum(latlon0.datum)._dup2(self) 

129 self._register(Conics, name) 

130 elif name: 

131 self.name = name 

132 if auth: 

133 self._auth = str(auth) 

134 

135 @Property_RO 

136 def auth(self): 

137 '''Get the authentication authority (C{str}). 

138 ''' 

139 return self._auth 

140 

141 @deprecated_method 

142 def convertDatum(self, datum): 

143 '''DEPRECATED, use method L{Conic.toDatum}.''' 

144 return self.toDatum(datum) 

145 

146 @Property_RO 

147 def datum(self): 

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

149 ''' 

150 return self._datum 

151 

152 @Property_RO 

153 def E0(self): 

154 '''Get the false easting (C{meter}). 

155 ''' 

156 return self._E0 

157 

158 @Property_RO 

159 def k0(self): 

160 '''Get scale factor (C{float}). 

161 ''' 

162 return self._k0 

163 

164 @Property_RO 

165 def lat0(self): 

166 '''Get the origin latitude (C{degrees90}). 

167 ''' 

168 return degrees90(self._phi0) 

169 

170 @Property_RO 

171 def latlon0(self): 

172 '''Get the central origin (L{LatLon2Tuple}C{(lat, lon)}). 

173 ''' 

174 return LatLon2Tuple(self.lat0, self.lon0, name=self.name) 

175 

176 @Property_RO 

177 def lam0(self): 

178 '''Get the central meridian (C{radians}). 

179 ''' 

180 return self._lam0 

181 

182 @Property_RO 

183 def lon0(self): 

184 '''Get the central meridian (C{degrees180}). 

185 ''' 

186 return degrees180(self._lam0) 

187 

188 @Property_RO 

189 def N0(self): 

190 '''Get the false northing (C{meter}). 

191 ''' 

192 return self._N0 

193 

194 @Property_RO 

195 def name2(self): 

196 '''Get the conic and datum names as "conic.datum" (C{str}). 

197 ''' 

198 return self._DOT_(self.datum.name) 

199 

200 @Property_RO 

201 def opt3(self): 

202 '''Get the optional meridian (C{degrees180}). 

203 ''' 

204 return degrees180(self._opt3) 

205 

206 @Property_RO 

207 def par1(self): 

208 '''Get the 1st standard parallel (C{degrees90}). 

209 ''' 

210 return degrees90(self._par1) 

211 

212 @Property_RO 

213 def par2(self): 

214 '''Get the 2nd standard parallel (C{degrees90}). 

215 ''' 

216 return degrees90(self._par2) 

217 

218 @Property_RO 

219 def phi0(self): 

220 '''Get the origin latitude (C{radians}). 

221 ''' 

222 return self._phi0 

223 

224 @Property_RO 

225 def philam0(self): 

226 '''Get the central origin (L{PhiLam2Tuple}C{(phi, lam)}). 

227 ''' 

228 return PhiLam2Tuple(self.phi0, self.lam0, name=self.name) 

229 

230 @Property_RO 

231 def SP(self): 

232 '''Get the number of standard parallels (C{int}). 

233 ''' 

234 return self._SP 

235 

236 def toDatum(self, datum): 

237 '''Convert this conic to the given datum. 

238 

239 @arg datum: Ellipsoidal datum to use (L{Datum}, L{Ellipsoid}, 

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

241 

242 @return: Converted conic, unregistered (L{Conic}). 

243 

244 @raise TypeError: Non-ellipsoidal B{C{datum}}. 

245 ''' 

246 d = _ellipsoidal_datum(datum, name=self.name) 

247 E = d.ellipsoid 

248 if not E.isEllipsoidal: 

249 raise _IsnotError(_ellipsoidal_, datum=datum) 

250 

251 c = self 

252 if c._e != E.e or c._datum != d: 

253 

254 c = Conic(None, 0, name=self._name) 

255 self._dup2(c) 

256 c._datum = d 

257 c._e = E.e 

258 

259 if fabs(c._par1 - c._par2) < EPS: 

260 m1 = c._mdef(c._phi0) 

261 t1 = c._tdef(c._phi0) 

262 t0 = t1 

263 k = 1 # _1_0 

264 n = sin(c._phi0) 

265 sp = 1 

266 else: 

267 m1 = c._mdef(c._par1) 

268 m2 = c._mdef(c._par2) 

269 t1 = c._tdef(c._par1) 

270 t2 = c._tdef(c._par2) 

271 t0 = c._tdef(c._phi0) 

272 k = c._k0 

273 n = (log(m1) - log(m2)) \ 

274 / (log(t1) - log(t2)) 

275 sp = 2 

276 

277 F = m1 / (n * pow(t1, n)) 

278 

279 c._aF = k * E.a * F 

280 c._n = n 

281 c._1_n = _1_0 / n 

282 c._r0 = c._rdef(t0) 

283 c._SP = sp 

284 

285 return c 

286 

287 def toStr(self, prec=8, name=NN, **unused): # PYCHOK expected 

288 '''Return this conic as a string. 

289 

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

291 @kwarg name: Override name (C{str}) or C{None} to exclude 

292 this conic's name. 

293 

294 @return: Conic attributes (C{str}). 

295 ''' 

296 a = [name, prec, _lat0_, _lon0_, _par1_, _par2_, 

297 _E0_, _N0_, _k0_, _SP_] 

298 if self._SP == 1: 

299 _ = a.pop(a.index(_par2_)) 

300 return self._instr(datum=self.datum, *a) 

301 

302 def _dup2(self, c): 

303 '''(INTERNAL) Copy this conic to C{c}. 

304 

305 @arg c: Duplicate (L{Conic}). 

306 ''' 

307 _update_all(c) 

308 

309 c._auth = self._auth 

310 c._datum = self._datum 

311 

312 c._e = self._e 

313 c._E0 = self._E0 

314 c._k0 = self._k0 

315 c._N0 = self._N0 

316 c._SP = self._SP 

317 

318 c._par1 = self._par1 

319 c._par2 = self._par2 

320 c._phi0 = self._phi0 

321 c._lam0 = self._lam0 

322 c._opt3 = self._opt3 

323 

324 c._aF = self._aF 

325 c._n = self._n 

326 c._1_n = self._1_n 

327 c._r0 = self._r0 

328 

329 def _mdef(self, a): 

330 '''(INTERNAL) Compute m(a). 

331 ''' 

332 s, c = sincos2(a) 

333 s = _1_0 - (s * self._e)**2 

334 return (c / sqrt(s)) if s > EPS02 else _0_0 

335 

336 def _pdef(self, a): 

337 '''(INTERNAL) Compute p(a). 

338 ''' 

339 s = self._e * sin(a) 

340 return pow((_1_0 - s) / (_1_0 + s), self._e * _0_5) 

341 

342 def _rdef(self, t): 

343 '''(INTERNAL) Compute r(t). 

344 ''' 

345 return self._aF * pow(t, self._n) 

346 

347 def _tdef(self, a): 

348 '''(INTERNAL) Compute t(lat). 

349 ''' 

350 return max(_0_0, tanPI_2_2(-a) / self._pdef(a)) 

351 

352 def _xdef(self, t_x): 

353 '''(INTERNAL) Compute x(t_x). 

354 ''' 

355 return PI_2 - _2_0 * atan(t_x) # XXX + self._phi0 

356 

357 

358Conic._name = Conic.__name__ 

359 

360 

361class Conics(_NamedEnum): 

362 '''(INTERNAL) L{Conic} registry, I{must} be a sub-class 

363 to accommodate the L{_LazyNamedEnumItem} properties. 

364 ''' 

365 def _Lazy(self, lat, lon, datum_name, *args, **kwds): 

366 '''(INTERNAL) Instantiate the L{Conic}. 

367 ''' 

368 return Conic(_LLEB(lat, lon, datum=Datums.get(datum_name)), *args, **kwds) 

369 

370Conics = Conics(Conic) # PYCHOK singleton 

371'''Some pre-defined L{Conic}s, all I{lazily} instantiated.''' 

372Conics._assert( # <https://SpatialReference.org/ref/sr-org/...> 

373# AsLb = _lazy('AsLb', _F(-14.2666667), _F(170), _NAD27_, _0_0, _0_0, 

374# E0=_F(500000), N0=_0_0, auth='EPSG:2155'), # American Samoa ... SP=1 ! 

375 Be08Lb = _lazy('Be08Lb', _F(50.7978150), _F(4.359215833), _GRS80_, _F(49.8333339), _F(51.1666672), 

376 E0=_F(649328.0), N0=_F(665262.0), auth='EPSG:9802'), # Belgium 

377 Be72Lb = _lazy('Be72Lb', _90_0, _F(4.3674867), _NAD83_, _F(49.8333339), _F(51.1666672), 

378 E0=_F(150000.013), N0=_F(5400088.438), auth='EPSG:31370'), # Belgium 

379 Fr93Lb = _lazy('Fr93Lb', _F(46.5), _F(3), _WGS84_, _F(49), _F(44), 

380 E0=_F(700000), N0=_F(6600000), auth='EPSG:2154'), # RFG93, France 

381 MaNLb = _lazy('MaNLb', _F(33.3), _F(-5.4), _NTF_, _F(31.73), _F(34.87), 

382 E0=_F(500000), N0=_F(300000)), # Marocco 

383 MxLb = _lazy('MxLb', _F(12), _F(-102), _WGS84_, _F(17.5), _F(29.5), 

384 E0=_F(2500000), N0=_0_0, auth='EPSG:2155'), # Mexico 

385 PyT_Lb = _lazy('PyT_Lb', _F(46.8), _F(2.33722917), _NTF_, _F(45.89893890000052), _F(47.69601440000037), 

386 E0=_F(600000), N0=_F(200000), auth='Test'), # France? 

387 USA_Lb = _lazy('USA_Lb', _F(23), _F(-96), _WGS84_, _F(33), _F(45), 

388 E0=_0_0, N0=_0_0), # Conterminous, contiguous USA? 

389 WRF_Lb = _lazy('WRF_Lb', _F(40), _F(-97), _WGS84_, _F(33), _F(45), 

390 E0=_0_0, N0=_0_0, auth='EPSG:4326') # World 

391) 

392 

393 

394class LCCError(_ValueError): 

395 '''Lambert Conformal Conic C{LCC} or other L{Lcc} issue. 

396 ''' 

397 pass 

398 

399 

400class Lcc(_NamedBase): 

401 '''Lambert conformal conic East-/Northing location. 

402 ''' 

403 _conic = Conics.WRF_Lb # Lambert projection (L{Conic}) 

404 _easting = _0_0 # Easting (C{float}) 

405 _height = 0 # height (C{meter}) 

406 _northing = _0_0 # Northing (C{float}) 

407 

408 def __init__(self, e, n, h=0, conic=Conics.WRF_Lb, name=NN): 

409 '''New L{Lcc} Lamber conformal conic position. 

410 

411 @arg e: Easting (C{meter}). 

412 @arg n: Northing (C{meter}). 

413 @kwarg h: Optional height (C{meter}). 

414 @kwarg conic: Optional, the conic projection (L{Conic}). 

415 @kwarg name: Optional name (C{str}). 

416 

417 @return: The Lambert location (L{Lcc}). 

418 

419 @raise LCCError: Invalid B{C{h}} or invalid or 

420 negative B{C{e}} or B{C{n}}. 

421 

422 @raise TypeError: If B{C{conic}} is not L{Conic}. 

423 

424 @example: 

425 

426 >>> lb = Lcc(448251, 5411932.0001) 

427 ''' 

428 if conic not in (None, Lcc._conic): 

429 self.conic = conic 

430 self._easting = Easting(e, falsed=conic.E0 > 0, Error=LCCError) 

431 self._northing = Northing(n, falsed=conic.N0 > 0, Error=LCCError) 

432 if h: 

433 self._height = Height(h=h, Error=LCCError) 

434 if name: 

435 self.name = name 

436 

437 @Property 

438 def conic(self): 

439 '''Get the conic projection (L{Conic}). 

440 ''' 

441 return self._conic 

442 

443 @conic.setter # PYCHOK setter! 

444 def conic(self, conic): 

445 '''Set the conic projection (L{Conic}). 

446 

447 @raise TypeError: Invalid B{C{conic}}. 

448 ''' 

449 _xinstanceof(Conic, conic=conic) 

450 if conic != self._conic: 

451 _update_all(self) 

452 self._conic = conic 

453 

454# def dup(self, name=NN, **e_n_h_conic): # PYCHOK signature 

455# '''Duplicate this location with some attributes modified. 

456# 

457# @kwarg e_n_h_conic: Use keyword argument C{B{e}=...}, C{B{n}=...}, 

458# C{B{h}=...} and/or C{B{conic}=...} to override 

459# the current C{easting}, C{northing} C{height} 

460# or C{conic} projection, respectively. 

461# ''' 

462# def _args_kwds(e=None, n=None, **kwds): 

463# return (e, n), kwds 

464# 

465# kwds = _xkwds(e_n_h_conic, e=self.easting, n=self.northing, 

466# h=self.height, conic=self.conic, 

467# name=name or self.name) 

468# args, kwds = _args_kwds(**kwds) 

469# return self.__class__(*args, **kwds) # .classof 

470 

471 @Property_RO 

472 def easting(self): 

473 '''Get the easting (C{meter}). 

474 ''' 

475 return self._easting 

476 

477 @Property_RO 

478 def height(self): 

479 '''Get the height (C{meter}). 

480 ''' 

481 return self._height 

482 

483 @Property_RO 

484 def latlon(self): 

485 '''Get the lat- and longitude in C{degrees} (L{LatLon2Tuple}). 

486 ''' 

487 ll = self.toLatLon(LatLon=None, datum=None) 

488 return LatLon2Tuple(ll.lat, ll.lon, name=self.name) 

489 

490 @Property_RO 

491 def latlonheight(self): 

492 '''Get the lat-, longitude and height (L{LatLon3Tuple}C{(lat, lon, height)}). 

493 ''' 

494 return self.latlon.to3Tuple(self.height) 

495 

496 @Property_RO 

497 def latlonheightdatum(self): 

498 '''Get the lat-, longitude in C{degrees} with height and datum (L{LatLon4Tuple}C{(lat, lon, height, datum)}). 

499 ''' 

500 return self.latlonheight.to4Tuple(self.conic.datum) 

501 

502 @Property_RO 

503 def northing(self): 

504 '''Get the northing (C{meter}). 

505 ''' 

506 return self._northing 

507 

508 @Property_RO 

509 def philam(self): 

510 '''Get the lat- and longitude in C{radians} (L{PhiLam2Tuple}). 

511 ''' 

512 return PhiLam2Tuple(radians(self.latlon.lat), 

513 radians(self.latlon.lon), name=self.name) 

514 

515 @Property_RO 

516 def philamheight(self): 

517 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}). 

518 ''' 

519 return self.philam.to3Tuple(self.height) 

520 

521 @Property_RO 

522 def philamheightdatum(self): 

523 '''Get the lat-, longitude in C{radians} with height and datum (L{PhiLam4Tuple}C{(phi, lam, height, datum)}). 

524 ''' 

525 return self.philamheight.to4Tuple(self.datum) 

526 

527 @deprecated_method 

528 def to3lld(self, datum=None): # PYCHOK no cover 

529 '''DEPRECATED, use method C{toLatLon}. 

530 

531 @kwarg datum: Optional datum to use, otherwise use this 

532 B{C{Lcc}}'s conic.datum (C{Datum}). 

533 

534 @return: A L{LatLonDatum3Tuple}C{(lat, lon, datum)}. 

535 

536 @raise TypeError: If B{C{datum}} is not ellipsoidal. 

537 ''' 

538 if datum in (None, self.conic.datum): 

539 r = LatLonDatum3Tuple(self.latlon.lat, 

540 self.latlon.lon, 

541 self.conic.datum, name=self.name) 

542 else: 

543 r = self.toLatLon(LatLon=None, datum=datum) 

544 r = LatLonDatum3Tuple(r.lat, r.lon, r.datum, name=r.name) 

545 return r 

546 

547 def toLatLon(self, LatLon=None, datum=None, height=None, **LatLon_kwds): 

548 '''Convert this L{Lcc} to an (ellipsoidal) geodetic point. 

549 

550 @kwarg LatLon: Optional, ellipsoidal class to return the 

551 geodetic point (C{LatLon}) or C{None}. 

552 @kwarg datum: Optional datum to use, otherwise use this 

553 B{C{Lcc}}'s conic.datum (L{Datum}, L{Ellipsoid}, 

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

555 @kwarg height: Optional height for the point, overriding 

556 the default height (C{meter}). 

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

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

559 

560 @return: The point (B{C{LatLon}}) or a 

561 L{LatLon4Tuple}C{(lat, lon, height, datum)} 

562 if B{C{LatLon}} is C{None}. 

563 

564 @raise TypeError: If B{C{LatLon}} or B{C{datum}} is 

565 not ellipsoidal or not valid. 

566 ''' 

567 if LatLon: 

568 _xsubclassof(_LLEB, LatLon=LatLon) 

569 

570 c = self.conic 

571 if datum not in (None, self.conic.datum): 

572 c = c.toDatum(datum) 

573 

574 e = self.easting - c._E0 

575 n = c._r0 - self.northing + c._N0 

576 

577 r_ = copysign0(hypot(e, n), c._n) 

578 t_ = pow(r_ / c._aF, c._1_n) 

579 

580 x = c._xdef(t_) # XXX c._lam0 

581 for self._iteration in range(10): # max 4 trips 

582 p, x = x, c._xdef(t_ * c._pdef(x)) 

583 if fabs(x - p) < 1e-9: # XXX EPS too small? 

584 break 

585 lat = degrees90(x) 

586 lon = degrees180((atan(e / n) + c._opt3) * c._1_n + c._lam0) 

587 

588 h = self.height if height is None else Height(height) 

589 return _LL4Tuple(lat, lon, h, c.datum, LatLon, LatLon_kwds, 

590 inst=self, name=self.name) 

591 

592 def toRepr(self, prec=0, fmt=Fmt.SQUARE, sep=_COMMASPACE_, m=_m_, C=False, **unused): # PYCHOK expected 

593 '''Return a string representation of this L{Lcc} position. 

594 

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

596 @kwarg fmt: Enclosing backets format (C{str}). 

597 @kwarg sep: Optional separator between name:values (C{str}). 

598 @kwarg m: Optional unit of the height, default meter (C{str}). 

599 @kwarg C: Optionally, include name of conic and datum (C{bool}). 

600 

601 @return: This Lcc as "[E:meter, N:meter, H:m, C:Conic.Datum]" 

602 (C{str}). 

603 ''' 

604 t, T = _fstrENH2(self, prec, m) 

605 if C: 

606 t += self.conic.name2, 

607 T += _C_, 

608 return _xzipairs(T, t, sep=sep, fmt=fmt) 

609 

610 def toStr(self, prec=0, sep=_SPACE_, m=_m_): # PYCHOK expected 

611 '''Return a string representation of this L{Lcc} position. 

612 

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

614 @kwarg sep: Optional separator to join (C{str}) or C{None} 

615 to return an unjoined C{tuple} of C{str}s. 

616 @kwarg m: Optional height units, default C{meter} (C{str}). 

617 

618 @return: This Lcc as I{"easting nothing"} in C{meter} plus 

619 I{" height"} suffixed with B{C{m}} if height is 

620 non-zero (C{str}). 

621 ''' 

622 t, _ = _fstrENH2(self, prec, m) 

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

624 

625 

626def toLcc(latlon, conic=Conics.WRF_Lb, height=None, Lcc=Lcc, name=NN, 

627 **Lcc_kwds): 

628 '''Convert an (ellipsoidal) geodetic point to a I{Lambert} location. 

629 

630 @arg latlon: Ellipsoidal point (C{LatLon}). 

631 @kwarg conic: Optional Lambert projection to use (L{Conic}). 

632 @kwarg height: Optional height for the point, overriding the 

633 default height (C{meter}). 

634 @kwarg Lcc: Optional class to return the I{Lambert} location 

635 (L{Lcc}). 

636 @kwarg name: Optional B{C{Lcc}} name (C{str}). 

637 @kwarg Lcc_kwds: Optional, additional B{C{Lcc}} keyword 

638 arguments, ignored if B{C{Lcc}} is C{None}. 

639 

640 @return: The I{Lambert} location (L{Lcc}) or an 

641 L{EasNor3Tuple}C{(easting, northing, height)} 

642 if C{B{Lcc} is None}. 

643 

644 @raise TypeError: If B{C{latlon}} is not ellipsoidal. 

645 ''' 

646 _xinstanceof(_LLEB, latlon=latlon) 

647 

648 a, b = latlon.philam 

649 c = conic.toDatum(latlon.datum) 

650 

651 t = c._n * (b - c._lam0) - c._opt3 

652 st, ct = sincos2(t) 

653 

654 r = c._rdef(c._tdef(a)) 

655 e = c._E0 + r * st 

656 n = c._N0 + c._r0 - r * ct 

657 

658 h = latlon.height if height is None else Height(height) 

659 r = EasNor3Tuple(e, n, h) if Lcc is None else \ 

660 Lcc(e, n, h=h, conic=c, **Lcc_kwds) 

661 return _xnamed(r, name or nameof(latlon)) 

662 

663 

664if __name__ == '__main__': 

665 

666 from pygeodesy.interns import _NL_, _NLATvar_ 

667 

668 # __doc__ of this file, force all into registery 

669 t = _NL_ + Conics.toRepr(all=True, asorted=True) 

670 print(_NLATvar_.join(t.split(_NL_))) 

671 

672# **) MIT License 

673# 

674# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

675# 

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

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

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

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

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

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

682# 

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

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

685# 

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

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

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

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

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

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

692# OTHER DEALINGS IN THE SOFTWARE.