1

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

3

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

5

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

7 and position class L{Lcc}.

8

9 See U{LCC<https://WikiPedia.org/wiki/Lambert_conformal_conic_projection>}, U{Lambert

10 Conformal Conic to Geographic Transformation Formulae

11<https://www.Linz.govt.NZ/data/geodetic-system/coordinate-conversion/projection-conversions/lambert-conformal-conic-geographic>},

12 U{Lambert Conformal Conic Projection<https://MathWorld.Wolfram.com/LambertConformalConicProjection.html>}

13 and 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

25 from __future__ import division as _; del _ # PYCHOK semicolon

26

27 from pygeodesy.basics import copysign0, _xinstanceof, _xsubclassof

28 from pygeodesy.constants import EPS, EPS02, PI_2, _float as _F, _0_0, _0_5, \

29 _1_0, _2_0, _90_0

30 from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB

31 from pygeodesy.datums import Datums, _ellipsoidal_datum

32 from pygeodesy.errors import _IsnotError, _ValueError

33 from pygeodesy.fmath import _ALL_LAZY, hypot

34 from 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

38 from pygeodesy.named import _lazyNamedEnumItem as _lazy, _NamedBase, \

39 _NamedEnum, _NamedEnumItem, nameof, _xnamed

40 from pygeodesy.namedTuples import EasNor3Tuple, LatLonDatum3Tuple, \

41 LatLon2Tuple, _LL4Tuple, PhiLam2Tuple

42 from pygeodesy.props import deprecated_method, Property, Property_RO, \

43 _update_all

44 from pygeodesy.streprs import Fmt, _fstrENH2, _xzipairs

45 from pygeodesy.units import Easting, Height, _heigHt, Lam_, Northing, \

46 Phi_, Scalar_

47 from pygeodesy.utily import atan1, degrees90, degrees180, sincos2, tanPI_2_2

48

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

50

51__all__ = _ALL_LAZY.lcc

52__version__ = '23.12.03'

53

54_E0_ = 'E0'

55_N0_ = 'N0'

56_par1_ = 'par1'

57_par2_ = 'par2'

58_SP_ = 'SP'

59

60

61 class 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

67

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})

73

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})

79

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})

84

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.

88

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}).

98

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

100

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

102

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}}.

106 '''

107 if latlon0 is not None:

108 _xinstanceof(_LLEB, latlon0=latlon0)

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

110

111 self._par1 = Phi_(par1=par1)

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

113

114 if k0 != 1:

115 self._k0 = Scalar_(k0=k0)

116 if E0:

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

118 if N0:

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

120 if opt3:

121 self._opt3 = Lam_(opt3=opt3)

122

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

124 self._register(Conics, name)

125 elif name:

126 self.name = name

127 if auth:

128 self._auth = str(auth)

129

130 @Property_RO

131 def auth(self):

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

133 '''

134 return self._auth

135

136 @deprecated_method

137 def convertDatum(self, datum):

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

139 return self.toDatum(datum)

140

141 @Property_RO

142 def datum(self):

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

144 '''

145 return self._datum

146

147 @Property_RO

148 def E0(self):

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

150 '''

151 return self._E0

152

153 @Property_RO

154 def k0(self):

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

156 '''

157 return self._k0

158

159 @Property_RO

160 def lat0(self):

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

162 '''

163 return degrees90(self._phi0)

164

165 @Property_RO

166 def latlon0(self):

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

168 '''

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

170

171 @Property_RO

172 def lam0(self):

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

174 '''

175 return self._lam0

176

177 @Property_RO

178 def lon0(self):

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

180 '''

181 return degrees180(self._lam0)

182

183 @Property_RO

184 def N0(self):

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

186 '''

187 return self._N0

188

189 @Property_RO

190 def name2(self):

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

192 '''

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

194

195 @Property_RO

196 def opt3(self):

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

198 '''

199 return degrees180(self._opt3)

200

201 @Property_RO

202 def par1(self):

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

204 '''

205 return degrees90(self._par1)

206

207 @Property_RO

208 def par2(self):

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

210 '''

211 return degrees90(self._par2)

212

213 @Property_RO

214 def phi0(self):

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

216 '''

217 return self._phi0

218

219 @Property_RO

220 def philam0(self):

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

222 '''

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

224

225 @Property_RO

226 def SP(self):

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

228 '''

229 return self._SP

230

231 def toDatum(self, datum):

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

233

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

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

236

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

238

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

240 '''

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

242 E = d.ellipsoid

243 if not E.isEllipsoidal:

244 raise _IsnotError(_ellipsoidal_, datum=datum)

245

246 c = self

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

248

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

250 self._dup2(c)

251 c._datum = d

252 c._e = E.e

253

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

255 m1 = c._mdef(c._phi0)

256 t1 = c._tdef(c._phi0)

257 t0 = t1

258 k = 1 # _1_0

259 n = sin(c._phi0)

260 sp = 1

261 else:

262 m1 = c._mdef(c._par1)

263 m2 = c._mdef(c._par2)

264 t1 = c._tdef(c._par1)

265 t2 = c._tdef(c._par2)

266 t0 = c._tdef(c._phi0)

267 k = c._k0

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

269 / (log(t1) - log(t2))

270 sp = 2

271

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

273

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

275 c._n = n

276 c._1_n = _1_0 / n

277 c._r0 = c._rdef(t0)

278 c._SP = sp

279

280 return c

281

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

283 '''Return this conic as a string.

284

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

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

287 this conic's name.

288

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

290 '''

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

292 _E0_, _N0_, _k0_, _SP_]

293 if self._SP == 1:

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

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

296

297 def _dup2(self, c):

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

299

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

301 '''

302 _update_all(c)

303

304 c._auth = self._auth

305 c._datum = self._datum

306

307 c._e = self._e

308 c._E0 = self._E0

309 c._k0 = self._k0

310 c._N0 = self._N0

311 c._SP = self._SP

312

313 c._par1 = self._par1

314 c._par2 = self._par2

315 c._phi0 = self._phi0

316 c._lam0 = self._lam0

317 c._opt3 = self._opt3

318

319 c._aF = self._aF

320 c._n = self._n

321 c._1_n = self._1_n

322 c._r0 = self._r0

323

324 def _mdef(self, a):

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

326 '''

327 s, c = sincos2(a)

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

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

330

331 def _pdef(self, a):

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

333 '''

334 s = self._e * sin(a)

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

336

337 def _rdef(self, t):

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

339 '''

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

341

342 def _tdef(self, a):

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

344 '''

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

346

347 def _xdef(self, t_x):

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

349 '''

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

351

352

353 Conic._name = Conic.__name__

354

355

356 class Conics(_NamedEnum):

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

358 to accommodate the L{_LazyNamedEnumItem} properties.

359 '''

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

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

362 '''

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

364

365 Conics = Conics(Conic) # PYCHOK singleton

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

386)

387

388

389 class LCCError(_ValueError):

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

391 '''

392 pass

393

394

395 class Lcc(_NamedBase):

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

397 '''

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

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

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

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

402

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

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

405

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

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

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

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

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

411

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

413

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

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

416

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

418 '''

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

420 self.conic = conic

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

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

423 if h:

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

425 if name:

426 self.name = name

427

428 @Property

429 def conic(self):

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

431 '''

432 return self._conic

433

434 @conic.setter # PYCHOK setter!

435 def conic(self, conic):

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

437

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

439 '''

440 _xinstanceof(Conic, conic=conic)

441 if conic != self._conic:

442 _update_all(self)

443 self._conic = conic

444

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

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

447#

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

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

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

451# or C{conic} projection, respectively.

452# '''

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

454# return (e, n), kwds

455#

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

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

458# name=name or self.name)

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

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

461

462 @Property_RO

463 def easting(self):

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

465 '''

466 return self._easting

467

468 @Property_RO

469 def height(self):

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

471 '''

472 return self._height

473

474 @Property_RO

475 def latlon(self):

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

477 '''

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

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

480

481 @Property_RO

482 def latlonheight(self):

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

484 '''

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

486

487 @Property_RO

488 def latlonheightdatum(self):

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

490 '''

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

492

493 @Property_RO

494 def northing(self):

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

496 '''

497 return self._northing

498

499 @Property_RO

500 def philam(self):

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

502 '''

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

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

505

506 @Property_RO

507 def philamheight(self):

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

509 '''

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

511

512 @Property_RO

513 def philamheightdatum(self):

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

515 '''

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

517

518 @deprecated_method

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

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

521

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

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

524

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

526

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

528 '''

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

530 r = LatLonDatum3Tuple(self.latlon.lat,

531 self.latlon.lon,

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

533 else:

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

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

536 return r

537

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

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

540

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

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

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

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

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

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

547 the default height (C{meter}).

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

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

550

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

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

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

554

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

556 not ellipsoidal or not valid.

557 '''

558 if LatLon:

559 _xsubclassof(_LLEB, LatLon=LatLon)

560

561 c = self.conic

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

563 c = c.toDatum(datum)

564

565 e = self.easting - c._E0

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

567

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

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

570

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

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

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

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

575 break

576 lat = degrees90(x)

577 lon = degrees180((atan1(e, n) + c._opt3) * c._1_n + c._lam0)

578

579 h = _heigHt(self, height)

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

581 inst=self, name=self.name)

582

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

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

585

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

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

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

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

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

591

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

593 (C{str}).

594 '''

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

596 if C:

597 t += self.conic.name2,

598 T += _C_,

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

600

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

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

603

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

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

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

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

608

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

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

611 non-zero (C{str}).

612 '''

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

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

615

616

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

618 **Lcc_kwds):

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

620

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

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

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

624 default height (C{meter}).

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

626 (L{Lcc}).

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

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

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

630

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

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

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

634

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

636 '''

637 _xinstanceof(_LLEB, latlon=latlon)

638

639 a, b = latlon.philam

640 c = conic.toDatum(latlon.datum)

641

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

643 st, ct = sincos2(t)

644

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

646 e = c._E0 + r * st

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

648

649 h = _heigHt(latlon, height)

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

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

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

653

654

655 if __name__ == '__main__':

656

657 from pygeodesy.interns import _NL_, _NLATvar_

658 from pygeodesy.lazily import printf

659

660 # __doc__ of this file, force all into registery

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

662 printf(_NLATvar_.join(t.split(_NL_)))

663

664# **) MIT License

665#

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

667#

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

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

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

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

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

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

674#

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

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

677#

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

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

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

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

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

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

684# OTHER DEALINGS IN THE SOFTWARE.