Coverage for pygeodesy/lcc.py: 97%

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 Conformal Conic to Geographic Transformation Formulae

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

11U{Lambert Conformal Conic Projection<https://MathWorld.Wolfram.com/LambertConformalConicProjection.html>}, U{Lambert Conic Conformal (1SP)

12<https://EPSG.io/9801-method>}, U{Lambert Conic Conformal (2SP)<https://EPSG.io/9802-method>} and

13John 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.83333333, par2=51.16666667, E0=649328, N0=665262, k0=1, SP=2, datum=Datum(name='GRS80', ellipsoid=Ellipsoids.GRS80, transform=Transforms.WGS84))

16@var Conics.Be50Lb: Conic(name='Be50Lb', lat0=90, lon0=4.36748667, par1=49.83333333, par2=51.16666667, E0=150000, N0=5400000, k0=1, SP=2, datum=Datum(name='ED50', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.ED50))

17@var Conics.Be72Lb: Conic(name='Be72Lb', lat0=90, lon0=4.36748667, par1=51.16666723, par2=49.8333339, E0=150000.013, N0=5400088.438, k0=1, SP=2, datum=Datum(name='ED50', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.ED50))

18@var Conics.Be72NLb: Conic(name='Be72NLb', lat0=50.79936222, lon0=4.36748667, par1=51.16666723, par2=49.8333339, E0=150000.013, N0=165372.956, k0=1, SP=2, datum=Datum(name='ED50', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.ED50))

19@var Conics.Be72RLb: Conic(name='Be72RLb', lat0=90, lon0=4.35693972, par1=49.83333333, par2=51.16666667, E0=150000.01256, N0=5400088.4378, k0=1, SP=2, datum=Datum(name='ED50', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.ED50))

20@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))

21@var Conics.MaNLb: Conic(name='MaNLb', lat0=33.3, lon0=-5.4, par1=33.3, E0=500000, N0=300000, k0=0.99962577, SP=1, datum=Datum(name='NTF', ellipsoid=Ellipsoids.Clarke1880IGN, transform=Transforms.NTF))

22@var Conics.MaSLb: Conic(name='MaSLb', lat0=29.7, lon0=-5.4, par1=29.7, E0=500000, N0=300000, k0=0.9996156, SP=1, datum=Datum(name='NTF', ellipsoid=Ellipsoids.Clarke1880IGN, transform=Transforms.NTF))

23@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))

24@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))

25@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))

26@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))

27'''

28# make sure int/int division yields float quotient, see .basics

29from __future__ import division as _; del _ # noqa: E702 ;

31from pygeodesy.basics import copysign0, _isin, typename, \

32 _xinstanceof, _xsubclassof

33from pygeodesy.constants import EPS, EPS02, PI_2, _float as _F, \

34 _0_0, _0_5, _1_0, _2_0, _90_0

35from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB

36from pygeodesy.datums import Datums, _ellipsoidal_datum

37from pygeodesy.errors import _IsnotError, _ValueError

38from pygeodesy.fmath import hypot, _ALL_LAZY

39# from pygeodesy.internals import typename # from .basics

40from pygeodesy.interns import NN, _COMMASPACE_, _DMAIN_, _ED50_, _ellipsoidal_, \

41 _GRS80_, _k0_, _lat0_, _lon0_, _m_, _NTF_, _SPACE_, \

42 _WGS84_, _C_ # PYCHOK used!

43# from pygeodesy.lazily import _ALL_LAZY # from .fmath

44from pygeodesy.named import _lazyNamedEnumItem as _lazy, _name2__, _NamedBase, \

45 _NamedEnum, _NamedEnumItem, _xnamed

46from pygeodesy.namedTuples import EasNor3Tuple, LatLonDatum3Tuple, \

47 LatLon2Tuple, _LL4Tuple, PhiLam2Tuple

48from pygeodesy.props import deprecated_method, Property, Property_RO, _update_all

49from pygeodesy.streprs import Fmt, _fstrENH2, _xzipairs

50from pygeodesy.units import Easting, Height, _heigHt, Lamd, Northing, Phid, \

51 Scalar_

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

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

56__all__ = _ALL_LAZY.lcc

57__version__ = '26.07.26'

59_E0_ = 'E0'

60_N0_ = 'N0'

61_par1_ = 'par1'

62_par2_ = 'par2'

63_SP_ = 'SP'

66class Conic(_NamedEnumItem):

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

68 '''

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

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

71 _name = NN # Conic.__name__, set below

73 _e = _0_0 # ellipsoid eccentricity (C{float})

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

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

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

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

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

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

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

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

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

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

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

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

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

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

91 k0=1, opt3=0, auth=NN, **name):

92 '''New Lambert conformal conic projection.

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

95 @arg par1: Optional, first standard parallel (C{degrees90}).

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

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

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

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

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

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

102 @kwarg name: Optional C{B{name}=NN} for the conic (C{str}).

103

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

105

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

107

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

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

110 or B{C{opt3}}.

111 '''

112 if latlon0 is not None:

113 _xinstanceof(_LLEB, latlon0=latlon0)

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

115

116 self._par1 = self._phi0 if par1 is None else Phid(par1=par1)

117 self._par2 = self._par1 if par2 is None else Phid(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 = Lamd(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 def forward3(self, latlon, lon=None, **height_name):

159 '''Convert I{geodetic} C{(lat, lon, height)} to (conformal) conic

160 C{easting} and C{northing)}.

161

162 @arg latlon: Either a C{LatLon}, L{Lcc} or C{scalar} (geodetic)

163 latitude (C{degrees}).

164 @kwarg lon: The C{scalar} (geodetic) longitude (C{degrees}), required

165 if B{C{latlonh}} is C{scalar}, ignored otherwise.

166 @kwarg height_name: Optionally, a C{B{height}=None} overriding the

167 default height (C{meter}) and C{B{name}=NN} (C{str}).

168

169 @return: An L{EasNor3Tuple}C{(easting, northing, height)}.

170

171 @raise LCCError: If B{C{latlonh}} not C{LatLon} L{Lcc}, C{scalar}

172 or invalid or if B{C{lon}} not C{scalar} for

173 C{scalar} B{C{latlonh}} or invalid or if B{C{height}}

174 invalid.

175 '''

176 try:

177 ll = _LLEB(latlon.lat, latlon.lon, datum=latlon.datum)

178 except AttributeError:

179 ll = _LLEB(latlon, lon, datum=self.datum)

180 return toLcc(ll, conic=self, Lcc=None, **height_name)

181

182 @Property_RO

183 def k0(self):

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

185 '''

186 return self._k0

187

188 @Property_RO

189 def lat0(self):

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

191 '''

192 return degrees90(self._phi0)

193

194 @Property_RO

195 def latlon0(self):

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

197 '''

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

199

200 @Property_RO

201 def lam0(self):

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

203 '''

204 return self._lam0

205

206 @Property_RO

207 def lon0(self):

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

209 '''

210 return degrees180(self._lam0)

211

212 @Property_RO

213 def N0(self):

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

215 '''

216 return self._N0

217

218 @Property_RO

219 def name2(self):

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

221 '''

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

223

224 @Property_RO

225 def opt3(self):

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

227 '''

228 return degrees180(self._opt3)

229

230 @Property_RO

231 def par1(self):

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

233 '''

234 return degrees90(self._par1)

235

236 @Property_RO

237 def par2(self):

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

239 '''

240 return degrees90(self._par2)

241

242 @Property_RO

243 def phi0(self):

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

245 '''

246 return self._phi0

247

248 @Property_RO

249 def philam0(self):

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

251 '''

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

253

254 def reverse4(self, enh, n=None, **height_name):

255 '''Convert I{conic} C{(easting, northing, height)} to geodetic C{lat-}

256 and C{longitude}.

257

258 @arg enh: Either an L{Lcc}, L{EasNor3Tuple} or C{scalar} easting

259 (C{meter}).

260 @kwarg n: The C{scalar} northing (C{meter}), required if B{C{enh}}

261 is C{scalar}, ignored otherwise.

262 @kwarg height_name: Optionally, a C{B{height}=None} overriding the

263 default height (C{meter}) and C{B{name}=NN} (C{str}).

264

265 @return: A L{LatLon4Tuple}C{(lat, lon, height, datum)} with this

266 conic's C{datum}.

267

268 @raise LCCError: If B{C{enh}} not L{Lcc}, L{EasNor3Tuple}, C{scalar}

269 or invalid or if B{C{n}} not C{scalar} for C{scalar}

270 B{C{n}} or invalid or if B{C{height}} invalid.

271 '''

272 lcc = enh if isinstance(enh, Lcc) else (

273 Lcc(enh.easting, enh.northing, conic=self) if n is None else

274 Lcc(enh, n, conic=self))

275 return lcc.toLatLon(LatLon=None, **height_name)

276

277 @Property_RO

278 def SP(self):

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

280 '''

281 return self._SP

282

283 def toDatum(self, datum):

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

285

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

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

288

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

290

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

292 '''

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

294 E = d.ellipsoid

295 if not E.isEllipsoidal:

296 raise _IsnotError(_ellipsoidal_, datum=datum)

297

298 c = self

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

300

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

302 self._dup2(c)

303 c._datum = d

304 c._e = E.e

305

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

307 m1 = c._mdef(c._phi0)

308 t1 = c._tdef(c._phi0)

309 t0 = t1

310 k = 1 # _1_0

311 n = sin(c._phi0)

312 sp = 1

313 else:

314 m1 = c._mdef(c._par1)

315 m2 = c._mdef(c._par2)

316 t1 = c._tdef(c._par1)

317 t2 = c._tdef(c._par2)

318 t0 = c._tdef(c._phi0)

319 k = c._k0

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

321 / (log(t1) - log(t2))

322 sp = 2

323

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

325

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

327 c._n = n

328 c._1_n = _1_0 / n

329 c._r0 = c._rdef(t0)

330 c._SP = sp

331

332 return c

333

334 def toStr(self, prec=8, **name): # PYCHOK expected

335 '''Return this conic as a string.

336

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

338 @kwarg name: Overriding C{B{name}=NN} (C{str}) or C{None} to

339 exclude this conic's name.

340

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

342 '''

343 a = [_lat0_, _lon0_, _par1_, _par2_, _E0_, _N0_, _k0_, _SP_] # , _auth_]

344 if self._SP == 1:

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

346 return self._instr(datum=self.datum, prec=prec, *a, **name)

347

348 def _dup2(self, c):

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

350

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

352 '''

353 _update_all(c)

354

355 c._auth = self._auth

356 c._datum = self._datum

357

358 c._e = self._e

359 c._E0 = self._E0

360 c._k0 = self._k0

361 c._N0 = self._N0

362 c._SP = self._SP

363

364 c._par1 = self._par1

365 c._par2 = self._par2

366 c._phi0 = self._phi0

367 c._lam0 = self._lam0

368 c._opt3 = self._opt3

369

370 c._aF = self._aF

371 c._n = self._n

372 c._1_n = self._1_n

373 c._r0 = self._r0

374

375 def _mdef(self, a):

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

377 '''

378 s, c = sincos2(a)

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

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

381

382 def _pdef(self, a):

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

384 '''

385 s = self._e * sin(a)

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

387

388 def _rdef(self, t):

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

390 '''

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

392

393 def _tdef(self, a):

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

395 '''

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

397

398 def _xdef(self, t_x):

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

400 '''

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

402

403

404Conic._name = typename(Conic)

405

406

407class Conics(_NamedEnum):

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

409 to accommodate the L{_LazyNamedEnumItem} properties.

410 '''

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

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

413 '''

414 ll0 = _LLEB(lat, lon, datum=Datums.get(datum_name))

415 return Conic(ll0, *args, **kwds)

416

417Conics = Conics(Conic) # PYCHOK singleton

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

419Conics._assert( # <https://SpatialReference.org/ref/epsg/...>

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

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

422 Be08Lb = _lazy('Be08Lb', _F(50.797815), _F(4.35921583333333), _GRS80_, _F(49.8333333333333), _F(51.1666666666667),

423 E0=_F(649328.0), N0=_F(665262.0), auth='EPSG:3812'), # ETRS89 / Belgian Lambert 2008

424 Be50Lb = _lazy('Be50Lb', _90_0, _F(4.36748666666667), _ED50_, _F(49.8333333333333), _F(51.1666666666667),

425 E0=_F(150000), N0=_F(5400000), auth='EPSG:21500'), # BD50 / Belgian Lambert 1950

426 Be72Lb = _lazy('Be72Lb', _90_0, _F(4.36748666666667), _ED50_, _F(51.1666672333333), _F(49.8333339),

427 E0=_F(150000.013), N0=_F(5400088.438), auth='EPSG:31370'), # BD72 / Belgian Lambert 1972

428 Be72NLb = _lazy('Be72NLb', '50 47 57.704N', _F(4.36748666666667), _ED50_, _F(51.1666672333333), _F(49.8333339),

429 E0=_F(150000.013), N0=_F(165372.956), auth='EPSG:31370'), # BD72 / lat0 not 90

430 Be72RLb = _lazy('Be72RLb', _90_0, _F(4.35693972222222), _ED50_, _F(49.8333333333333), _F(51.1666666666667),

431 E0=_F(150000.01256), N0=_F(5400088.4378), auth='EPSG:19902'), # BD72 / "rotated" == EPSG:31300

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

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

434 MaNLb = _lazy('MaNLb', _F(33.3), _F(-5.4), _NTF_, # par1 == lat0

435 E0=_F(500000), N0=_F(300000), k0=_F(0.999625769), auth='EPSG:26191'), # Nord Maroc

436 MaSLb = _lazy('MaSLb', _F(29.7), _F(-5.4), _NTF_, # par1 == lat0

437 E0=_F(500000), N0=_F(300000), k0=_F(0.999615596), auth='EPSG:26192'), # Sud Maroc

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

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

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

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

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

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

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

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

446)

447

448

449class LCCError(_ValueError):

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

451 '''

452 pass

453

454

455class Lcc(_NamedBase):

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

457 '''

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

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

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

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

462

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

464 '''New L{Lcc} Lambert conformal conic position.

465

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

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

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

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

470 @kwarg name: Optional C{B{name}=NN} (C{str}).

471

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

473

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

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

476

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

478 '''

479 if not _isin(conic, None, Lcc._conic):

480 self.conic = conic

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

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

483 if h:

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

485 if name:

486 self.name = name

487

488 @Property

489 def conic(self):

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

491 '''

492 return self._conic

493

494 @conic.setter # PYCHOK setter!

495 def conic(self, conic):

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

497

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

499 '''

500 _xinstanceof(Conic, conic=conic)

501 if conic != self._conic:

502 _update_all(self)

503 self._conic = conic

504

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

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

507#

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

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

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

511# or C{conic} projection, respectively.

512# '''

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

514# return (e, n), kwds

515#

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

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

518# name=self._name__(name))

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

520# return type(self)(*args, **kwds) # .classof

521

522 @Property_RO

523 def easting(self):

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

525 '''

526 return self._easting

527

528 @Property_RO

529 def height(self):

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

531 '''

532 return self._height

533

534 @Property_RO

535 def latlon(self):

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

537 '''

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

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

540

541 @Property_RO

542 def latlonheight(self):

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

544 '''

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

546

547 @Property_RO

548 def latlonheightdatum(self):

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

550 '''

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

552

553 @Property_RO

554 def northing(self):

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

556 '''

557 return self._northing

558

559 @Property_RO

560 def philam(self):

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

562 '''

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

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

565

566 @Property_RO

567 def philamheight(self):

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

569 '''

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

571

572 @Property_RO

573 def philamheightdatum(self):

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

575 '''

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

577

578 @deprecated_method

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

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

581

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

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

584

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

586

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

588 '''

589 if _isin(datum, None, self.conic.datum):

590 r = LatLonDatum3Tuple(self.latlon.lat,

591 self.latlon.lon,

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

593 else:

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

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

596 return r

597

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

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

600

601 @kwarg LatLon: Optional, ellipsoidal class to return the geodetic

602 point (C{LatLon}) or C{None}.

603 @kwarg datum: Optional datum to use, otherwise use this B{C{Lcc}}'s

604 conic.datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or

605 L{a_f2Tuple}).

606 @kwarg height: Optional height for the point, overriding the default height

607 (C{meter}).

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

609 ignored if C{B{LatLon} is None}.

610

611 @return: The point (B{C{LatLon}}) or if C{B{LatLon} is None}, a

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

613

614 @raise TypeError: If B{C{LatLon}} or B{C{datum}} is not ellipsoidal or

615 not valid.

616 '''

617 if LatLon:

618 _xsubclassof(_LLEB, LatLon=LatLon)

619

620 c = self.conic

621 if not _isin(datum, None, self.conic.datum):

622 c = c.toDatum(datum)

623

624 e = self.easting - c._E0

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

626

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

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

629

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

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

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

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

634 break

635 lat = degrees90(x)

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

637

638 h = _heigHt(self, height)

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

640 inst=self, name=self.name)

641

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

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

644

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

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

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

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

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

650

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

652 (C{str}).

653 '''

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

655 if C:

656 t += self.conic.name2,

657 T += _C_,

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

659

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

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

662

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

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

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

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

667

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

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

670 non-zero (C{str}).

671 '''

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

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

674

675

676def toLcc(latlon, conic=Conics.WRF_Lb, height=None, Lcc=Lcc,

677 **name_Lcc_kwds):

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

679

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

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

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

683 default height (C{meter}).

684 @kwarg Lcc: Class to return the I{Lambert} location (L{Lcc}).

685 @kwarg name_Lcc_kwds: Optional C{B{name}=NN} (C{str}) and optionally,

686 additional B{C{Lcc}} keyword arguments, ignored if B{C{Lcc}

687 is None}.

688

689 @return: The I{Lambert} location (L{Lcc}) or if C{B{Lcc} is None},

690 an L{EasNor3Tuple}C{(easting, northing, height)}.

691

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

693 '''

694 _xinstanceof(_LLEB, latlon=latlon)

695 name, Lcc_kwds = _name2__(name_Lcc_kwds)

696

697 a, b = latlon.philam

698 c = conic.toDatum(latlon.datum)

699

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

701 st, ct = sincos2(t)

702

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

704 e = c._E0 + r * st

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

706

707 h = _heigHt(latlon, height)

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

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

710 return _xnamed(r, name) if name else r

711

712

713if __name__ == _DMAIN_:

714 # __doc__ of this file, force all into registry

715 from pygeodesy.internals import _pregistry

716 _pregistry(Conics)

717

718# **) MIT License

719#

721#

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

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

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

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

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

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

728#

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

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

731#

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

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

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

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

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

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

738# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pygeodesy / lcc.py: 97%

288 statements