Coverage for pygeodesy/ltp.py: 95%

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

4u'''I{Local Tangent Plane} (LTP) and I{local} cartesian coordinates.

6I{Local cartesian} and I{local tangent plane} classes L{LocalCartesian}, approximations L{ChLVa}

7and L{ChLVe} and L{Ltp}, L{ChLV}, L{LocalError}, L{Attitude} and L{Frustum}.

9@see: U{Local tangent plane coordinates<https://WikiPedia.org/wiki/Local_tangent_plane_coordinates>}

10 and class L{LocalCartesian}, transcoded from I{Charles Karney}'s C++ classU{LocalCartesian

11 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1LocalCartesian.html>}.

12'''

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

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

16from pygeodesy.basics import _args_kwds_names, _isin, map1, map2, _xinstanceof, \

17 _xsubclassof, typename # .datums

18from pygeodesy.constants import EPS, INT0, _umod_360, _0_0, _0_01, _0_5, _1_0, \

19 _2_0, _8_0, _60_0, _90_0, _100_0, _180_0, _3600_0, \

20 _N_1_0 # PYCHOK used!

21from pygeodesy.datums import Datums, _WGS84

22# from pygeodesy.dms import parseDMS # from .units

23from pygeodesy.ecef import _EcefBase, EcefKarney, Ecef9Tuple, _llhn4, _xyzn4

24from pygeodesy.errors import _NotImplementedError, _ValueError, _xattr, \

25 _xkwds, _xkwds_get, _xkwds_pop2

26from pygeodesy.fmath import fabs, fdot, Fdot_, fdot_, Fhorner

27from pygeodesy.fsums import _floor, fsumf_

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

29from pygeodesy.interns import _0_, _COMMASPACE_, _DOT_, _ecef_, _height_, \

30 _invalid_, _lat0_, _lon0_, _M_, _name_, _too_

31# from pygeodesy.lazily import _ALL_LAZY # from vector3d

32from pygeodesy.ltpTuples import Attitude4Tuple, ChLVEN2Tuple, ChLV9Tuple, \

33 ChLVYX2Tuple, Footprint5Tuple, Local9Tuple, \

34 ChLVyx2Tuple, _XyzLocals4, _XyzLocals5, Xyz4Tuple

35from pygeodesy.named import _name__, _name2__, _NamedBase, notOverloaded

36from pygeodesy.namedTuples import LatLon3Tuple, LatLon4Tuple, Vector3Tuple, \

37 Bounds4Tuple # PYCHOK used!

38from pygeodesy.props import Property, Property_RO, property_doc_, \

39 property_ROver, _update_all

40from pygeodesy.streprs import Fmt, strs, unstr

41from pygeodesy.units import Bearing, Degrees, _isHeight, Meter, parseDMS

42from pygeodesy.utily import cotd, _loneg, sincos2d, sincos2d_, tand, tand_, \

43 wrap180, wrap360

44from pygeodesy.vector3d import _ALL_LAZY, Vector3d

46# from math import fabs, floor as _floor # from .fmath, .fsums

48__all__ = _ALL_LAZY.ltp

49__version__ = '26.05.24'

51_GRS80 = Datums.GRS80

52_height0_ = _height_ + _0_

53_narrow_ = 'narrow'

54_wide_ = 'wide'

55# del Datums

58class Attitude(_NamedBase):

59 '''The pose of a plane or camera in space.

60 '''

61 _alt = Meter( alt =_0_0)

62 _roll = Degrees(roll=_0_0)

63 _tilt = Degrees(tilt=_0_0)

64 _yaw = Bearing(yaw =_0_0)

66 def __init__(self, alt_attitude=INT0, tilt=INT0, yaw=INT0, roll=INT0, **name):

67 '''New L{Attitude}.

69 @kwarg alt_attitude: Altitude (C{meter}) above earth or a previous attitude

70 (L{Attitude} or L{Attitude4Tuple}) with the C{B{alt}itude},

71 B{C{tilt}}, B{C{yaw}} and B{C{roll}}.

72 @kwarg tilt: Pitch, elevation from horizontal (C{degrees180}), negative down

73 (clockwise rotation along and around the x- or East axis), iff

74 B{C{alt_attitude}} is C{meter}, ignored otherwise.

75 @kwarg yaw: Bearing, heading (compass C{degrees360}), clockwise from North

76 (counter-clockwise rotation along and around the z- or Up axis)

77 iff B{C{alt_attitude}} is C{meter}, ignored otherwise.

78 @kwarg roll: Roll, bank (C{degrees180}), positive to the right and down

79 (clockwise rotation along and around the y- or North axis), iff

80 B{C{alt_attitude}} is C{meter}, ignored otherwise.

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

83 @raise AttitudeError: Invalid B{C{alt_attitude}}, B{C{tilt}}, B{C{yaw}} or

84 B{C{roll}}.

86 @see: U{Principal axes<https://WikiPedia.org/wiki/Aircraft_principal_axes>} and

87 U{Yaw, pitch, and roll rotations<http://MSL.CS.UIUC.edu/planning/node102.html>}.

88 '''

89 if _isHeight(alt_attitude):

90 t = Attitude4Tuple(alt_attitude, tilt, yaw, roll)

91 else:

92 try:

93 t = alt_attitude.atyr

94 except AttributeError:

95 raise AttitudeError(alt=alt_attitude, tilt=tilt, yaw=yaw, rol=roll)

96 for n, v in t.items():

97 if v:

98 setattr(self, n, v)

99 n = _name__(name, _or_nameof=t)

100 if n:

101 self.name = n

102

103 @property_doc_(' altitude above earth in C{meter}.')

104 def alt(self):

105 return self._alt

106

107 @alt.setter # PYCHOK setter!

108 def alt(self, alt): # PYCHOK no cover

109 a = Meter(alt=alt, Error=AttitudeError)

110 if self._alt != a:

111 _update_all(self)

112 self._alt = a

113

114 altitude = alt

115

116 @Property_RO

117 def atyr(self):

118 '''Return this attitude's alt[itude], tilt, yaw and roll as an L{Attitude4Tuple}.

119 '''

120 return Attitude4Tuple(self.alt, self.tilt, self.yaw, self.roll, name=self.name)

121

122 @Property_RO

123 def matrix(self):

124 '''Get the 3x3 rotation matrix C{R(yaw)·R(tilt)·R(roll)}, aka I{ZYX} (C{float}, row-order).

125

126 @see: Matrix M of case 10 in U{Appendix A

127 <https://ntrs.NASA.gov/api/citations/19770019231/downloads/19770019231.pdf>}.

128 '''

129 # to follow the definitions of rotation angles alpha, beta and gamma:

130 # negate yaw since yaw is counter-clockwise around the z-axis, swap

131 # tilt and roll since tilt is around the x- and roll around the y-axis

132 sa, ca, sb, cb, sg, cg = sincos2d_(-self.yaw, self.roll, self.tilt)

133 return ((ca * cb, fdot_(ca, sb * sg, -sa, cg), fdot_(ca, sb * cg, sa, sg)),

134 (sa * cb, fdot_(sa, sb * sg, ca, cg), fdot_(sa, sb * cg, -ca, sg)),

135 ( -sb, cb * sg, cb * cg))

136

137 @property_doc_(' roll/bank in C{degrees180}, positive to the right and down.')

138 def roll(self):

139 return self._roll

140

141 @roll.setter # PYCHOK setter!

142 def roll(self, roll):

143 r = Degrees(roll=roll, wrap=wrap180, Error=AttitudeError)

144 if self._roll != r:

145 _update_all(self)

146 self._roll = r

147

148 bank = roll

149

150 def rotate(self, x_xyz, y=None, z=None, Vector=None, **name_Vector_kwds):

151 '''Transform a (local) cartesian by this attitude's matrix.

152

153 @arg x_xyz: X component of vector (C{scalar}) or (3-D) vector (C{Cartesian},

154 L{Vector3d} or L{Vector3Tuple}).

155 @kwarg y: Y component of vector (C{scalar}), same units as C{scalar} B{C{x}},

156 ignored otherwise.

157 @kwarg z: Z component of vector (C{scalar}), same units as C{sclar} B{C{x}},

158 ignored otherwise.

159 @kwarg Vector: Class to return transformed point (C{Cartesian}, L{Vector3d}

160 or C{Vector3Tuple}) or C{None}.

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

162 additional B{C{Vector}} keyword arguments, ignored if C{B{Vector}

163 is None}.

164

165 @return: A named B{C{Vector}} instance or if C{B{Vector} is None},

166 a named L{Vector3Tuple}C{(x, y, z)}.

167

168 @raise AttitudeError: Invalid B{C{x_xyz}}, B{C{y}} or B{C{z}}.

169

170 @raise TypeError: Invalid B{C{Vector}} or B{C{name_Vector_kwds}} item.

171

172 @see: U{Yaw, pitch, and roll rotations<http://MSL.CS.UIUC.edu/planning/node102.html>}.

173 '''

174 try:

175 try:

176 xyz = map2(float, x_xyz.xyz3)

177 except AttributeError:

178 xyz = map1(float, x_xyz, y, z)

179 except (TypeError, ValueError) as x:

180 raise AttitudeError(x_xyz=x_xyz, y=y, z=z, cause=x)

181

182 x, y, z = (fdot(r, *xyz) for r in self.matrix)

183 n, kwds = _name2__(name_Vector_kwds, _or_nameof=self)

184 return Vector3Tuple(x, y, z, name=n) if Vector is None else \

185 Vector(x, y, z, name=n, **kwds)

186

187 @property_doc_(' tilt/pitch/elevation from horizontal in C{degrees180}, negative down.')

188 def tilt(self):

189 return self._tilt

190

191 @tilt.setter # PYCHOK setter!

192 def tilt(self, tilt):

193 t = Degrees(tilt=tilt, wrap=wrap180, Error=AttitudeError)

194 if self._tilt != t:

195 _update_all(self)

196 self._tilt = t

197

198 elevation = pitch = tilt

199

200 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature

201 '''Format this attitude as string.

202

203 @kwarg prec: The C{float} precision, number of decimal digits (0..9).

204 Trailing zero decimals are stripped for B{C{prec}} values

205 of 1 and above, but kept for negative B{C{prec}} values.

206 @kwarg sep: Separator to join (C{str}).

207

208 @return: This attitude (C{str}).

209 '''

210 return self.atyr.toStr(prec=prec, sep=sep)

211

212 @Property_RO

213 def tyr3d(self):

214 '''Get this attitude's (3-D) directional vector (L{Vector3d}).

215

216 @see: U{Yaw, pitch, and roll rotations<http://MSL.CS.UIUC.edu/planning/node102.html>}.

217 '''

218 def _r2d(r):

219 return fsumf_(_N_1_0, *r)

220

221 return Vector3d(*map(_r2d, self.matrix), name__=tyr3d)

222

223 @property_doc_(' yaw/bearing/heading in compass C{degrees360}, clockwise from North.')

224 def yaw(self):

225 return self._yaw

226

227 @yaw.setter # PYCHOK setter!

228 def yaw(self, yaw):

229 y = Bearing(yaw=yaw, Error=AttitudeError)

230 if self._yaw != y:

231 _update_all(self)

232 self._yaw = y

233

234 bearing = heading = yaw # azimuth

235

236

237class AttitudeError(_ValueError):

238 '''An L{Attitude} or L{Attitude4Tuple} issue.

239 '''

240 pass

241

242

243class Frustum(_NamedBase):

244 '''A rectangular pyramid, typically representing a camera's I{field-of-view}

245 (fov) and the intersection with (or projection to) a I{local tangent plane}.

246

247 @see: U{Viewing frustum<https://WikiPedia.org/wiki/Viewing_frustum>}.

248 '''

249 _h_2 = _0_0 # half hfov in degrees

250 _ltp = None # local tangent plane

251 _tan_h_2 = _0_0 # tan(_h_2)

252 _v_2 = _0_0 # half vfov in degrees

253

254 def __init__(self, hfov, vfov, ltp=None, **name):

255 '''New L{Frustum}.

256

257 @arg hfov: Horizontal field-of-view (C{degrees180}).

258 @arg vfov: Vertical field-of-view (C{degrees180}).

259 @kwarg ltp: Optional I{local tangent plane} (L{Ltp}).

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

261

262 @raise LocalError: Invalid B{C{hfov}} or B{C{vfov}}.

263 '''

264 self._h_2 = h = _fov_2(hfov=hfov)

265 self._v_2 = _fov_2(vfov=vfov)

266

267 self._tan_h_2 = tand(h, hfov_2=h)

268

269 if ltp:

270 self._ltp = _xLtp(ltp)

271 if name:

272 self.name # PYCHOK effect

273

274 def footprint5(self, alt_attitude, tilt=0, yaw=0, roll=0, z=_0_0, ltp=None, **name): # MCCABE 15

275 '''Compute the center and corners of the intersection with (or projection

276 to) the I{local tangent plane} (LTP).

277

278 @arg alt_attitude: An altitude (C{meter}) above I{local tangent plane} or

279 an attitude (L{Attitude} or L{Attitude4Tuple}) with the

280 C{B{alt}itude}, B{C{tilt}}, B{C{yaw}} and B{C{roll}}.

281 @kwarg tilt: Pitch, elevation from horizontal (C{degrees}), negative down

282 (clockwise rotation along and around the x- or East axis) iff

283 B{C{alt_attitude}} is C{meter}, ignored otherwise.

284 @kwarg yaw: Bearing, heading (compass C{degrees}), clockwise from North

285 (counter-clockwise rotation along and around the z- or Up axis)

286 iff B{C{alt_attitude}} is C{meter}, ignored otherwise.

287 @kwarg roll: Roll, bank (C{degrees}), positive to the right and down

288 (clockwise rotation along and around the y- or North axis) iff

289 B{C{alt_attitude}} is C{meter}, ignored otherwise.

290 @kwarg z: Optional height of the footprint (C{meter}) above I{local tangent plane}.

291 @kwarg ltp: The I{local tangent plane} (L{Ltp}), overriding this

292 frustum's C{ltp}.

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

294

295 @return: A L{Footprint5Tuple}C{(center, upperleft, upperight, loweright,

296 lowerleft)} with the C{center} and 4 corners, each an L{Xyz4Tuple}.

297

298 @raise TypeError: Invalid B{C{ltp}}.

299

300 @raise UnitError: Invalid B{C{altitude}}, B{C{tilt}}, B{C{roll}} or B{C{z}}.

301

302 @raise ValueError: If B{C{altitude}} too low, B{C{z}} too high or B{C{tilt}}

303 or B{C{roll}} -including B{C{vfov}} respectively B{C{hfov}}-

304 over the horizon.

305

306 @see: U{Principal axes<https://WikiPedia.org/wiki/Aircraft_principal_axes>}.

307 '''

308 def _xy2(a, e, h_2, tan_h_2, r):

309 # left and right corners, or swapped

310 if r < EPS: # no roll

311 r = a * tan_h_2

312 l = -r # noqa: E741 l is ell

313 else: # roll

314 r, l = tand_(r - h_2, r + h_2, roll_hfov=r) # noqa: E741 l is ell

315 r *= -a # negate right positive

316 l *= -a # noqa: E741 l is ell

317 y = a * cotd(e, tilt_vfov=e)

318 return (l, y), (r, y)

319

320 def _xyz5(b, xy5, z, ltp):

321 # rotate (x, y)'s by bearing, clockwise

322 sc = sincos2d(b)

323 for x, y in xy5:

324 yield Xyz4Tuple(fdot(sc, x, y),

325 fdot(sc, -x, y), z, ltp)

326

327 try:

328 a, t, y, r = alt_attitude.atyr

329 except AttributeError:

330 a, t, y, r = alt_attitude, tilt, yaw, roll

331

332 a = Meter(altitude=a)

333 if a < EPS: # too low

334 raise _ValueError(altitude=a)

335 if z: # PYCHOK no cover

336 z = Meter(z=z)

337 a -= z

338 if a < EPS: # z above a

339 raise _ValueError(altitude_z=a)

340 else:

341 z = _0_0

342

343 b = Degrees(yaw=y, wrap=wrap360) # bearing

344 e = -Degrees(tilt=t, wrap=wrap180) # elevation, pitch

345 if not EPS < e < _180_0:

346 raise _ValueError(tilt=t)

347 if e > _90_0:

348 e = _loneg(e)

349 b = _umod_360(b + _180_0)

350

351 r = Degrees(roll=r, wrap=wrap180) # roll center

352 x = (-a * tand(r, roll=r)) if r else _0_0

353 y = a * cotd(e, tilt=t) # ground range

354 if fabs(y) < EPS:

355 y = _0_0

356

357 v, h, t = self._v_2, self._h_2, self._tan_h_2

358 # center and corners, clockwise from upperleft, rolled

359 xy5 = ((x, y),) + _xy2(a, e - v, h, t, r) \

360 + _xy2(a, e + v, -h, -t, r) # swapped

361 # turn center and corners by yaw, clockwise

362 p = self.ltp if ltp is None else ltp # None OK

363 return Footprint5Tuple(_xyz5(b, xy5, z, p), **name) # *_xyz5

364

365 @Property_RO

366 def hfov(self):

367 '''Get the horizontal C{fov} (C{degrees}).

368 '''

369 return Degrees(hfov=self._h_2 * _2_0)

370

371 @Property_RO

372 def ltp(self):

373 '''Get the I{local tangent plane} (L{Ltp}) or C{None}.

374 '''

375 return self._ltp

376

377 def toStr(self, prec=3, fmt=Fmt.F, sep=_COMMASPACE_): # PYCHOK signature

378 '''Convert this frustum to a "hfov, vfov, ltp" string.

379

380 @kwarg prec: Number of (decimal) digits, unstripped (0..8 or C{None}).

381 @kwarg fmt: Optional, C{float} format (C{letter}).

382 @kwarg sep: Separator to join (C{str}).

383

384 @return: Frustum in the specified form (C{str}).

385 '''

386 t = self.hfov, self.vfov

387 if self.ltp:

388 t += self.ltp,

389 t = strs(t, prec=prec, fmt=fmt)

390 return sep.join(t) if sep else t

391

392 @Property_RO

393 def vfov(self):

394 '''Get the vertical C{fov} (C{degrees}).

395 '''

396 return Degrees(vfov=self._v_2 * _2_0)

397

398

399class LocalError(_ValueError):

400 '''A L{LocalCartesian} or L{Ltp} related issue.

401 '''

402 pass

403

404

405class LocalCartesian(_NamedBase):

406 '''Conversion between geodetic C{(lat, lon, height)} and I{local

407 cartesian} C{(x, y, z)} coordinates with I{geodetic} origin

408 C{(lat0, lon0, height0)}, transcoded from I{Karney}'s C++ class

409 U{LocalCartesian<https://GeographicLib.SourceForge.io/C++/doc/

410 classGeographicLib_1_1LocalCartesian.html>}.

411

412 The C{z} axis is normal to the ellipsoid, the C{y} axis points due

413 North. The plane C{z = -height0} is tangent to the ellipsoid.

414

415 The conversions all take place via geocentric coordinates using a

416 geocentric L{EcefKarney}, by default the WGS84 datum/ellipsoid.

417

418 @see: Class L{Ltp}.

419 '''

420 _Ecef = EcefKarney

421 _ecef = EcefKarney(_WGS84)

422 _lon00 = INT0 # self.lon0

423 _9t0 = None # origin (..., lat0, lon0, height0, ...) L{Ecef9Tuple}

424 _9Tuple = Local9Tuple

425

426 def __init__(self, latlonh0=INT0, lon0=INT0, height0=INT0, ecef=None, **lon00_name):

427 '''New L{LocalCartesian} converter.

428

429 @kwarg latlonh0: The (geodetic) origin (C{LatLon}, L{LatLon4Tuple}, L{Ltp}

430 L{LocalCartesian} or L{Ecef9Tuple}) or the C{scalar}

431 latitude of the (goedetic) origin (C{degrees}).

432 @kwarg lon0: Longitude of the (goedetic) origin (C{degrees}), required if

433 B{C{latlonh0}} is C{scalar}, ignored otherwise.

434 @kwarg height0: Optional height (C{meter}, conventionally) at the (goedetic)

435 origin perpendicular to and above (or below) the ellipsoid's

436 surface, like B{C{lon0}}.

437 @kwarg ecef: An ECEF converter (L{EcefKarney} I{only}), like B{C{lon0}}.

438 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and keyword argument

439 C{B{lon00}=B{lon0}} for the arbitrary I{polar} longitude

440 (C{degrees}), see method C{reverse} and property C{lon00}

441 for further details.

442

443 @raise LocalError: If B{C{latlonh0}} not C{LatLon}, L{LatLon4Tuple}, L{Ltp},

444 L{LocalCartesian} or L{Ecef9Tuple} or B{C{latlonh0}},

445 B{C{lon0}}, B{C{height0}} or B{C{lon00}} invalid.

446

447 @raise TypeError: Invalid B{C{ecef}} or not L{EcefKarney}.

448

449 @note: If BC{latlonh0} is an L{Ltp} or L{LocalCartesian}, only C{lat0}, C{lon0},

450 C{height0} and I{polar} C{lon00} are copied, I{not} the ECEF converter.

451 '''

452 self.reset(latlonh0, lon0=lon0, height0=height0, ecef=ecef, **lon00_name)

453

454 def __eq__(self, other):

455 '''Compare this and an other instance.

456

457 @arg other: The other ellipsoid (L{LocalCartesian} or L{Ltp}).

458

459 @return: C{True} if equal, C{False} otherwise.

460 '''

461 return other is self or (isinstance(other, self.__class__) and

462 other.ecef == self.ecef and

463 other._9t0 == self._9t0)

464

465 @Property_RO

466 def datum(self):

467 '''Get the ECEF converter's datum (L{Datum}).

468 '''

469 return self.ecef.datum

470

471 @Property_RO

472 def ecef(self):

473 '''Get the ECEF converter (L{EcefKarney}).

474 '''

475 return self._ecef

476

477 def _ecef2local(self, ecef, Xyz, name_Xyz_kwds): # in _EcefLocal._Ltp_ecef2local

478 '''(INTERNAL) Convert geocentric/geodetic to local, like I{forward}.

479

480 @arg ecef: Geocentric (and geodetic) (L{Ecef9Tuple}).

481 @arg Xyz: An L{XyzLocal}, L{Aer}, L{Enu} or L{Ned} I{class} or C{None}.

482 @arg name_Xyz_kwds: Optional C{B{name}=NN} (C{str}) and optionally,

483 additional B{C{Xyz}} keyword arguments, ignored if C{B{Xyz}

484 is None}.

485

486 @return: An C{B{Xyz}(x, y, z, ltp, **B{name_Xyz_kwds}} instance or

487 if C{B{Xyz} is None}, a L{Local9Tuple}C{(x, y, z, lat, lon,

488 height, ltp, ecef, M)} with this C{ltp}, B{C{ecef}}

489 (L{Ecef9Tuple}) converted to this C{datum} and C{M=None},

490 always.

491

492 @raise TypeError: Invalid B{C{Xyz}} or B{C{name_Xyz_kwds}} item.

493 '''

494 _xinstanceof(Ecef9Tuple, ecef=ecef)

495 if ecef.datum != self.datum:

496 ecef = ecef.toDatum(self.datum)

497 n, kwds = _name2__(name_Xyz_kwds, _or_nameof=ecef)

498 x, y, z = self.M.rotate(ecef.xyz, *self._9t0_xyz)

499 r = Local9Tuple(x, y, z, ecef.lat, ecef.lon, ecef.height,

500 self, ecef, None, name=n)

501 if Xyz:

502 _xsubclassof(*_XyzLocals4, Xyz=Xyz) # Vector3d

503 r = r.toXyz(Xyz=Xyz, name=n, **kwds)

504 return r

505

506 @Property_RO

507 def ellipsoid(self):

508 '''Get the ECEF converter's ellipsoid (L{Ellipsoid}).

509 '''

510 return self.ecef.datum.ellipsoid

511

512 def forward(self, latlonh, lon=None, height=0, M=False, **name):

513 '''Convert I{geodetic} C{(lat, lon, height)} to I{local} cartesian

514 C{(x, y, z)}.

515

516 @arg latlonh: Either a C{LatLon}, L{Ltp}, L{Ecef9Tuple} or C{scalar}

517 (geodetic) latitude (C{degrees}).

518 @kwarg lon: Optional C{scalar} (geodetic) longitude (C{degrees}) iff

519 B{C{latlonh}} is C{scalar}, ignored otherwise.

520 @kwarg height: Optional height (C{meter}, conventionally) perpendicular

521 to and above (or below) the ellipsoid's surface, iff

522 B{C{latlonh}} is C{scalar}, ignored otherwise.

523 @kwarg M: Optionally, return the I{concatenated} rotation L{EcefMatrix},

524 iff available (C{bool}).

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

526

527 @return: A L{Local9Tuple}C{(x, y, z, lat, lon, height, ltp, ecef, M)}

528 with I{local} C{x}, C{y}, C{z}, I{geodetic} C{(lat}, C{lon},

529 C{height}, this C{ltp}, C{ecef} (L{Ecef9Tuple}) with

530 I{geocentric} C{x}, C{y}, C{z} (and I{geodetic} C{lat},

531 C{lon}, C{height}) and the I{concatenated} rotation matrix

532 C{M} (L{EcefMatrix}) if requested.

533

534 @raise LocalError: If B{C{latlonh}} not C{scalar}, C{LatLon}, L{Ltp},

535 L{Ecef9Tuple} or invalid or if B{C{lon}} not

536 C{scalar} for C{scalar} B{C{latlonh}} or invalid

537 or if B{C{height}} invalid.

538 '''

539 lat, lon, h, n = _llhn4(latlonh, lon, height, Error=LocalError, **name)

540 t = self.ecef._forward(lat, lon, h, n, M=M)

541 x, y, z = self.M.rotate(t.xyz, *self._9t0_xyz)

542 m = self.M.multiply(t.M) if M else None

543 return self._9Tuple(x, y, z, lat, lon, h, self, t, m, name=n or self.name)

544

545 @Property_RO

546 def height0(self):

547 '''Get the origin's height (C{meter}).

548 '''

549 return self._9t0.height

550

551 @Property_RO

552 def lat0(self):

553 '''Get the origin's latitude (C{degrees}).

554 '''

555 return self._9t0.lat

556

557 @Property_RO

558 def latlonheight0(self):

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

560 '''

561 return LatLon3Tuple(self.lat0, self.lon0, self.height0, name=self.name)

562

563 def _local2ecef(self, local, nine=False, M=False):

564 '''(INTERNAL) Convert I{local} to geocentric/geodetic, like I{.reverse}.

565

566 @arg local: Local (L{XyzLocal}, L{Enu}, L{Ned}, L{Aer} or L{Local9Tuple}).

567 @kwarg nine: If C{True}, return a 9-, otherwise a 3-tuple (C{bool}).

568 @kwarg M: Include the rotation matrix (C{bool}).

569

570 @return: A I{geocentric} 3-tuple C{(x, y, z)} or if C{B{nine}=True}, an

571 L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with

572 rotation matrix C{M} (L{EcefMatrix}) if requested.

573 '''

574 _xinstanceof(*_XyzLocals5, local=local)

575 t = self.M.unrotate(local.xyz, *self._9t0_xyz)

576 if nine:

577 t = self.ecef.reverse(*t, M=M)

578 return t

579

580 @Property_RO

581 def lon0(self):

582 '''Get the origin's longitude (C{degrees}).

583 '''

584 return self._9t0.lon

585

586 @Property

587 def lon00(self):

588 '''Get the arbitrary, I{polar} longitude (C{degrees}).

589 '''

590 return self._lon00

591

592 @lon00.setter # PYCHOK setter!

593 def lon00(self, lon00):

594 '''Set the arbitrary, I{polar} longitude (C{degrees}).

595 '''

596 # lon00 <https://GitHub.com/mrJean1/PyGeodesy/issues/77>

597 self._lon00 = Degrees(lon00=lon00)

598

599 @Property_RO

600 def M(self):

601 '''Get the rotation matrix (C{EcefMatrix}).

602 '''

603 return self._9t0.M

604

605 def reset(self, latlonh0=INT0, lon0=INT0, height0=INT0, ecef=None, **lon00_name):

606 '''Reset this converter, see L{LocalCartesian.__init__} for further details.

607 '''

608 _, name = _xkwds_pop2(lon00_name, lon00=None) # PYCHOK get **name

609 if isinstance(latlonh0, LocalCartesian):

610 if self._9t0:

611 _update_all(self)

612 self._ecef = latlonh0.ecef

613 self._lon00 = latlonh0.lon00

614 self._9t0 = latlonh0._9t0

615 n = _name__(name, _or_nameof=latlonh0)

616 else:

617 n = _name__(name, _or_nameof=self)

618 lat0, lon0, height0, n = _llhn4(latlonh0, lon0, height0, suffix=_0_,

619 Error=LocalError, name=n)

620 if ecef: # PYCHOK no cover

621 _xinstanceof(self._Ecef, ecef=ecef)

622 _update_all(self)

623 self._ecef = ecef

624 elif self._9t0:

625 _update_all(self)

626 self._9t0 = self.ecef._forward(lat0, lon0, height0, n, M=True)

627 self.lon00 = _xattr(latlonh0, lon00=_xkwds_get(lon00_name, lon00=lon0))

628 if n:

629 self.rename(n)

630

631 def reverse(self, xyz, y=None, z=None, M=False, **lon00_name):

632 '''Convert I{local} C{(x, y, z)} to I{geodetic} C{(lat, lon, height)}.

633

634 @arg xyz: A I{local} (L{XyzLocal}, L{Enu}, L{Ned}, L{Aer}, L{Local9Tuple}) or

635 local C{x} coordinate (C{scalar}).

636 @kwarg y: Local C{y} coordinate (C{meter}), iff B{C{xyz}} is C{scalar},

637 ignored otherwise.

638 @kwarg z: Local C{z} coordinate (C{meter}), iff B{C{xyz}} is C{scalar},

639 ignored otherwise.

640 @kwarg M: Optionally, return the I{concatenated} rotation L{EcefMatrix}, iff

641 available (C{bool}).

642 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and keyword argument

643 C{B{lon00}=B{lon0}} for the arbitrary I{polar} longitude

644 (C{degrees}), overriding see the property C{B{lon00}=B{lon0}}

645 value. The I{polar} longitude (C{degrees}) is returned with

646 I{polar} latitudes C{abs(B{lat0}) == 90} for local C{B{x}=0}

647 and C{B{y}=0} locations.

648

649 @return: A L{Local9Tuple}C{(x, y, z, lat, lon, height, ltp, ecef, M)} with

650 I{local} C{x}, C{y}, C{z}, I{geodetic} C{lat}, C{lon}, C{height},

651 this C{ltp}, an C{ecef} (L{Ecef9Tuple}) with the I{geocentric} C{x},

652 C{y}, C{z} (and I{geodetic} C{lat}, C{lon}, C{height}) and the

653 I{concatenated} rotation matrix C{M} (L{EcefMatrix}) if requested.

654

655 @raise LocalError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}

656 not C{scalar} for C{scalar} B{C{xyz}}.

657 '''

658 lon00, name =_xkwds_pop2(lon00_name, lon00=self.lon00)

659 x, y, z, n = _xyzn4(xyz, y, z, _XyzLocals5, Error=LocalError, name=name)

660 c = self.M.unrotate((x, y, z), *self._9t0_xyz)

661 t = self.ecef.reverse(*c, M=M, lon00=lon00)

662 m = self.M.multiply(t.M) if M else None

663 return self._9Tuple(x, y, z, t.lat, t.lon, t.height, self, t, m, name=n or self.name)

664

665 @Property_RO

666 def _9t0_xyz(self):

667 '''(INTERNAL) Get C{(x0, y0, z0)} as L{Vector3Tuple}.

668 '''

669 return self._9t0.xyz

670

671 def toStr(self, prec=9, **unused): # PYCHOK signature

672 '''Return this L{LocalCartesian} as a string.

673

674 @kwarg prec: Precision, number of (decimal) digits (0..9).

675

676 @return: This L{LocalCartesian} representation (C{str}).

677 '''

678 return self.attrs(_lat0_, _lon0_, _height0_, _M_, _ecef_, _name_, prec=prec)

679

680

681class Ltp(LocalCartesian):

682 '''A I{local tangent plan} (LTP), a sub-class of C{LocalCartesian} with

683 (re-)configurable ECEF converter.

684 '''

685 _Ecef = _EcefBase

686

687 def __init__(self, latlonh0=INT0, lon0=INT0, height0=INT0, ecef=None, **lon00_name):

688 '''New C{Ltp}, see L{LocalCartesian.__init__} for more details.

689

690 @kwarg ecef: Optional ECEF converter (L{EcefKarney}, L{EcefFarrell21},

691 L{EcefFarrell22}, L{EcefSudano}, L{EcefVeness} or L{EcefYou}

692 I{instance}), overriding the default L{EcefKarney}C{(datum=Datums.WGS84)}

693 for C{scalar} B{C{latlonh0}}.

694

695 @see: Class L{LocalCartesian<LocalCartesian.__init__>} for further details.

696

697 @raise TypeError: Invalid B{C{ecef}}.

698 '''

699 LocalCartesian.reset(self, latlonh0, lon0=lon0, height0=height0,

700 ecef=ecef, **lon00_name)

701

702 @Property

703 def ecef(self):

704 '''Get this LTP's ECEF converter (C{Ecef...} I{instance}).

705 '''

706 return self._ecef

707

708 @ecef.setter # PYCHOK setter!

709 def ecef(self, ecef):

710 '''Set this LTP's ECEF converter (C{Ecef...} I{instance}).

711

712 @raise TypeError: Invalid B{C{ecef}}.

713 '''

714 _xinstanceof(_EcefBase, ecef=ecef)

715 if self._ecef != ecef: # PYCHOK no cover

716 self.reset(self._9t0)

717 self._ecef = ecef

718

719

720class LqRD(Ltp):

721 '''A I{local tangent plan} (LTP) for conversion between I{GRS80 (ETRS89) geodetic} and

722 I{local Netherlands}' C{quasi-B{R}ijksB{D}riehoeksmeting (RD)} coordinates.

723

724 This C{quasi-RD} transformer B{does not} implement any U{RD NAP<https://www.NSGI.NL/

725 coordinatenstelsels-en-transformaties/coordinatentransformaties/rdnap-etrs89-rdnaptrans>}

726 specification and B{does not} provide I{Netherlands}' C{B{N}ormaal B{A}msterdams B{P}eil

727 (NAP)} quasi-geodetic-height.

728

729 The L{LqRD.forward} C{x} and C{y} results differ 3 meter near the center up to 600 meter

730 at the corners of the L{RD region<LqRD.region>} with C{RDx} and C{RDy} values from formal

731 C{RD NAP 2018} implementations like U{pyrdnap<https://PyPI.org/project/pyrdnap>}.

732

733 The L{LqRD.forward} C{z} values represent perpendicular distances to this local tangent

734 plane (LTP). Other heights in L{LqRD} are I{GRS80 (ETRS89) heights} above (or below)

735 the ellipsoid's surface. B{None} are C{NAP} quasi-geodetic-heights.

736 '''

737 Amersfoort = LatLon4Tuple(parseDMS('52 9 22.178N'), # height=0.0, not .height0_ETRS!

738 parseDMS(' 5 23 15.5E'), _0_0, _GRS80, name='Amersfoort')

739# _Ecef = _EcefBase

740 _ecef = EcefKarney(_GRS80)

741 _x0 = Meter(x0=155029.8) # 155000.0 see pyrdnap -v1 -forward Amersfoort.latlon

742 _y0 = Meter(y0=463109.9) # 463000.0 see pyrdnap -v1 -forward Amersfoort.latlon

743

744 def __init__(self, latlonh0=Amersfoort, **other_Ltp_kwds):

745 '''New ECEF-based I{GRS80 (ETRS89)} L{LqRD} converter, centered at I{Amersfoort, NL}.

746

747 @kwarg latlonh0: The I{geodetic} origin and height, overriding C{Amersfoort}.

748 @kwarg other_Ltp_kwds: Optional, other L{Ltp.__init__} keyword arguments.

749

750 @see: Class L{Ltp<Ltp.__init__>} for more information.

751 '''

752 Ltp.__init__(self, latlonh0, **_xkwds(other_Ltp_kwds, ecef=None, name=LqRD.Amersfoort.name))

753

754 def forward(self, lat_latlonh, lon=None, height=0, **M_name): # PYCHOK signature

755 '''Convert I{geodetic} C{(lat, lon, height)} to I{local} C{quasi-RD (x, y, z)}.

756

757 @return: A L{Local9Tuple}C{(x, y, z, lat, lon, height, ltp, ecef, M)}.

758

759 @see: Method L{LocalCartesian.forward} for more information.

760 '''

761 r = Ltp.forward(self, lat_latlonh, lon, height, **M_name)

762 return r.dup(x=r.x + self.x0, y=r.y + self.y0, name=r.name)

763

764 @property_ROver

765 def height0_ETRS(self):

766 '''Get C{Amersfoort}'s I{GRS80 (ETRS89) height} (C{Meter}).

767 '''

768 return Meter(height0_ETRS=43.0) # see pyrdnap h0_ETRS

769

770 @property_ROver

771 def region(self):

772 '''Get the C{RD} region as L{Bounds4Tuple}C{(latS, lonW, latN, lonE)}, all C{GRS80 (ETRS89) degrees}.

773 '''

774 return Bounds4Tuple(50.0, _2_0, 56.0, _8_0).toUnits(name='RD region ') # like pyrdnap

775

776 def reverse(self, x_xyz, y=None, z=None, **M_name): # PYCHOK signature

777 '''Convert I{local} C{quasi-RD (x, y, z)} to I{geodetic} C{(lat, lon, height)}.

778

779 @return: A L{Local9Tuple}C{(x, y, z, lat, lon, height, ltp, ecef, M)}.

780

781 @see: Method L{LocalCartesian.reverse} for more information.

782 '''

783 x, y, z = x_xyz.xyz if y is z is None else map1(Meter, x_xyz, y, z)

784 r = Ltp.reverse(self, x - self.x0, y - self.y0, z, **M_name)

785 return r.dup(x=x, y=y, name=r.name)

786

787 @property_doc_(' the C{quasi-RD} false Easting (C{meter}).')

788 def x0(self):

789 return self._x0

790

791 @x0.setter # PYCHOK setter!

792 def x0(self, meter):

793 self._x0 = Meter(x0=meter)

794

795 @property_doc_(' the C{quasi-RD} false Northing (C{meter}).')

796 def y0(self):

797 return self._y0

798

799 @y0.setter # PYCHOK setter!

800 def y0(self, meter):

801 self._y0 = Meter(y0=meter)

802

803

804class _ChLV(object):

805 '''(INTERNAL) Base class for C{ChLV*} classes.

806 '''

807 _03_falsing = ChLVyx2Tuple(0.6e6, 0.2e6)

808# _92_falsing = ChLVYX2Tuple(2.0e6, 1.0e6) # _95_ - _03_

809 _95_falsing = ChLVEN2Tuple(2.6e6, 1.2e6)

810

811 def _ChLV9Tuple(self, fw, M, name, *Y_X_h_lat_lon_h):

812 '''(INTERNAL) Helper for C{ChLVa/e.forward} and C{.reverse}.

813 '''

814 if bool(M): # PYCHOK no cover

815 m = self.forward if fw else self.reverse # PYCHOK attr

816 n = _DOT_(*map1(typename, type(self), m))

817 raise _NotImplementedError(unstr(n, M=M), txt=None)

818 t = Y_X_h_lat_lon_h + (self, self._9t0, None) # PYCHOK _9t0

819 return ChLV9Tuple(t, name=name)

820

821 @property_ROver

822 def _enh_n_h(self):

823 '''(INTERNAL) Get C{ChLV*.reverse} args[1:4] names, I{once}.

824 '''

825 t = _args_kwds_names(_ChLV.reverse)[1:4]

826 # assert _args_kwds_names( ChLV.reverse)[1:4] == t

827 # assert _args_kwds_names(ChLVa.reverse)[1:4] == t

828 # assert _args_kwds_names(ChLVe.reverse)[1:4] == t

829 return t # overwrite property_ROver

830

831 def forward(self, latlonh, lon=None, height=0, M=None, **name): # PYCHOK no cover

832 '''Convert WGS84 geodetic to I{Swiss} projection coordinates. I{Must be overloaded}.

833

834 @arg latlonh: Either a C{LatLon}, L{Ltp} or C{scalar} (geodetic) latitude (C{degrees}).

835 @kwarg lon: Optional, C{scalar} (geodetic) longitude (C{degrees}) iff B{C{latlonh}} is

836 C{scalar}, ignored otherwise.

837 @kwarg height: Optional, height, vertically above (or below) the surface of the ellipsoid

838 (C{meter}) iff B{C{latlonh}} and B{C{lon}} are C{scalar}, ignored otherwise.

839 @kwarg M: If C{True}, return the I{concatenated} rotation L{EcefMatrix} iff available

840 and for C{ChLV} only, C{None} otherwise (C{bool}).

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

842

843 @return: A L{ChLV9Tuple}C{(Y, X, h_, lat, lon, height, ltp, ecef, M)} with the unfalsed

844 I{Swiss Y, X} coordinates, I{Swiss h_} height, the given I{geodetic} C{lat},

845 C{lon} and C{height}, C{ecef} (L{Ecef9Tuple}) at I{Bern, Ch}, rotation matrix

846 C{M} and C{ltp} this C{ChLV}, C{ChLVa} or C{ChLVe} instance.

847

848 @raise LocalError: Invalid or non-C{scalar} B{C{latlonh}}, B{C{lon}} or B{C{height}}.

849 '''

850 notOverloaded(self, latlonh, lon=lon, height=height, M=M, **name)

851

852 def reverse(self, enh_, n=None, h_=0, M=None, **name): # PYCHOK no cover

853 '''Convert I{Swiss} projection to WGS84 geodetic coordinates.

854

855 @arg enh_: A Swiss projection (L{ChLV9Tuple}) or the C{scalar}, falsed I{Swiss E_LV95}

856 or I{y_LV03} easting (C{meter}).

857 @kwarg n: Falsed I{Swiss N_LV85} or I{x_LV03} northing (C{meter}) iff B{C{enh_}} is

858 C{scalar}, ignored otherwise.

859 @kwarg h_: I{Swiss h'} height (C{meter}) iff B{C{enh_}} and B{C{n}} are C{scalar},

860 ignored otherwise.

861 @kwarg M: If C{True}, return the I{concatenated} rotation L{EcefMatrix} iff available

862 and for C{ChLV} only, C{None} otherwise (C{bool}).

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

864

865 @return: A L{ChLV9Tuple}C{(Y, X, h_, lat, lon, height, ltp, ecef, M)} with the unfalsed

866 I{Swiss Y, X} coordinates, I{Swiss h_} height, the given I{geodetic} C{lat},

867 C{lon} and C{height}, C{ecef} (L{Ecef9Tuple}) at I{Bern, Ch}, rotation matrix

868 C{M} and C{ltp} this C{ChLV}, C{ChLVa} or C{ChLVe} instance.

869

870 @raise LocalError: Invalid or non-C{scalar} B{C{enh_}}, B{C{n}} or B{C{h_}}.

871 '''

872 notOverloaded(self, enh_, n=n, h_=h_, M=M, **name)

873

874 @staticmethod

875 def _falsing2(LV95):

876 '''(INTERNAL) Get the C{LV95} or C{LV03} falsing.

877 '''

878 return _ChLV._95_falsing if _isin(LV95, True, 95) else (

879 _ChLV._03_falsing if _isin(LV95, False, 3) else ChLVYX2Tuple(0, 0))

880

881 @staticmethod

882 def _llh2abh_3(lat, lon, h):

883 '''(INTERNAL) Helper for C{ChLVa/e.forward}.

884 '''

885 def _deg2ab(deg, sLL):

886 # convert degrees to arc-seconds

887 def _dms(ds, p, q, swap):

888 d = _floor(ds)

889 t = (ds - d) * p

890 m = _floor(t)

891 s = (t - m) * p

892 if swap:

893 d, s = s, d

894 return d + (m + s * q) * q

895

896 s = _dms(deg, _60_0, _0_01, False) # deg2sexag

897 s = _dms( s, _100_0, _60_0, True) # sexag2asec

898 return (s - sLL) / ChLV._s_ab

899

900 a = _deg2ab(lat, ChLV._sLat) # phi', lat_aux

901 b = _deg2ab(lon, ChLV._sLon) # lam', lng_aux

902 h_ = fsumf_(h, -ChLV.Bern.height, 2.73 * b, 6.94 * a)

903 return a, b, h_

904

905 @staticmethod

906 def _YXh_2abh3(Y, X, h_):

907 '''(INTERNAL) Helper for C{ChLVa/e.reverse}.

908 '''

909 def _YX2ab(YX):

910 return YX * ChLV._ab_m

911

912 a, b = map1(_YX2ab, Y, X)

913 h = fsumf_(h_, ChLV.Bern.height, -12.6 * a, -22.64 * b)

914 return a, b, h

915

916 def _YXh_n4(self, enh_, n, h_, **name):

917 '''(INTERNAL) Helper for C{ChLV*.reverse}.

918 '''

919 Y, X, h_, name = _xyzn4(enh_, n, h_, (ChLV9Tuple,),

920 _xyz_y_z_names=self._enh_n_h, **name)

921 if isinstance(enh_, ChLV9Tuple):

922 Y, X = enh_.Y, enh_.X

923 else: # isscalar(enh_)

924 Y, X = ChLV.unfalse2(Y, X) # PYCHOK ChLVYX2Tuple

925 return Y, X, h_, name

926

927

928class ChLV(_ChLV, Ltp):

929 '''Conversion between I{WGS84 geodetic} and I{Swiss} projection coordinates using

930 L{pygeodesy.EcefKarney}'s Earth-Centered, Earth-Fixed (ECEF) methods.

931

932 @see: U{Swiss projection formulas<https://www.SwissTopo.admin.CH/en/maps-data-online/

933 calculation-services.html>}, page 7ff, U{NAVREF<https://www.SwissTopo.admin.CH/en/

934 maps-data-online/calculation-services/navref.html>}, U{REFRAME<https://www.SwissTopo.admin.CH/

935 en/maps-data-online/calculation-services/reframe.html>} and U{SwissTopo Scripts GPS WGS84

936 <-> LV03<https://GitHub.com/ValentinMinder/Swisstopo-WGS84-LV03>}.

937 '''

938 _9Tuple = ChLV9Tuple

939

940 _ab_d = 0.36 # a, b units per degree, ...

941 _ab_m = 1.0e-6 # ... per meter and ...

942 _ab_M = _1_0 # ... per 1,000 Km or 1 Mm

943 _s_d = _3600_0 # arc-seconds per degree ...

944 _s_ab = _s_d / _ab_d # ... and per a, b unit

945 _sLat = 169028.66 # Bern, Ch in ...

946 _sLon = 26782.5 # ... arc-seconds ...

947 # lat, lon, height == 46°57'08.66", 7°26'22.50", 49.55m ("new" 46°57'07.89", 7°26'22.335")

948 Bern = LatLon4Tuple(_sLat / _s_d, _sLon / _s_d, 49.55, _WGS84, name='Bern')

949

950 def __init__(self, latlonh0=Bern, **other_Ltp_kwds):

951 '''New ECEF-based I{WGS84-Swiss} L{ChLV} converter, centered at I{Bern, Ch}.

952

953 @kwarg latlonh0: The I{geodetic} origin and height, overriding C{Bern, Ch}.

954 @kwarg other_Ltp_kwds: Optional, other L{Ltp.__init__} keyword arguments.

955

956 @see: L{Ltp.__init__} for more information.

957 '''

958 Ltp.__init__(self, latlonh0, **_xkwds(other_Ltp_kwds, ecef=None, name=ChLV.Bern.name))

959

960 def forward(self, latlonh, lon=None, height=0, M=None, **name): # PYCHOK unused M

961 # overloaded for the _ChLV.forward.__doc__

962 return Ltp.forward(self, latlonh, lon=lon, height=height, M=M, **name)

963

964 def reverse(self, enh_, n=None, h_=0, M=None, **name): # PYCHOK signature

965 # overloaded for the _ChLV.reverse.__doc__

966 Y, X, h_, n = self._YXh_n4(enh_, n, h_, **name)

967 return Ltp.reverse(self, Y, X, h_, M=M, name=n)

968

969 @staticmethod

970 def false2(Y, X, LV95=True, **name):

971 '''Add the I{Swiss LV95} or I{LV03} falsing.

972

973 @arg Y: Unfalsed I{Swiss Y} easting (C{meter}).

974 @arg X: Unfalsed I{Swiss X} northing (C{meter}).

975 @kwarg LV95: If C{True}, add C{LV95} falsing, if C{False} add

976 C{LV03} falsing, otherwise leave unfalsed.

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

978

979 @return: A L{ChLVEN2Tuple}C{(E_LV95, N_LV95)} or a

980 L{ChLVyx2Tuple}C{(y_LV03, x_LV03)} with falsed B{C{Y}}

981 and B{C{X}}, otherwise a L{ChLVYX2Tuple}C{(Y, X)}

982 with B{C{Y}} and B{C{X}} as-is.

983 '''

984 e, n = t = _ChLV._falsing2(LV95)

985 return t.classof(e + Y, n + X, **name)

986

987 @staticmethod

988 def isLV03(e, n):

989 '''Is C{(B{e}, B{n})} a valid I{Swiss LV03} projection?

990

991 @arg e: Falsed (or unfalsed) I{Swiss} easting (C{meter}).

992 @arg n: Falsed (or unfalsed) I{Swiss} northing (C{meter}).

993

994 @return: C{True} if C{(B{e}, B{n})} is a valid, falsed I{Swiss

995 LV03}, projection C{False} otherwise.

996 '''

997 # @see: U{Map<https://www.SwissTopo.admin.CH/en/knowledge-facts/

998 # surveying-geodesy/reference-frames/local/lv95.html>}

999 return 400.0e3 < e < 900.0e3 and 40.0e3 < n < 400.0e3

1000

1001 @staticmethod

1002 def isLV95(e, n, raiser=True):

1003 '''Is C{(B{e}, B{n})} a valid I{Swiss LV95} or I{LV03} projection?

1004

1005 @arg e: Falsed (or unfalsed) I{Swiss} easting (C{meter}).

1006 @arg n: Falsed (or unfalsed) I{Swiss} northing (C{meter}).

1007 @kwarg raiser: If C{True}, throw a L{LocalError} if B{C{e}} and

1008 B{C{n}} are invalid I{Swiss LV95} nor I{LV03}.

1009

1010 @return: C{True} or C{False} if C{(B{e}, B{n})} is a valid I{Swiss

1011 LV95} respectively I{LV03} projection, C{None} otherwise.

1012 '''

1013 if ChLV.isLV03(e, n):

1014 return False

1015 elif ChLV.isLV03(e - 2.0e6, n - 1.0e6): # _92_falsing = _95_ - _03_

1016 return True

1017 elif raiser: # PYCHOK no cover

1018 raise LocalError(unstr(ChLV.isLV95, e=e, n=n))

1019 return None

1020

1021 @staticmethod

1022 def unfalse2(e, n, LV95=None, **name):

1023 '''Remove the I{Swiss LV95} or I{LV03} falsing.

1024

1025 @arg e: Falsed I{Swiss E_LV95} or I{y_LV03} easting (C{meter}).

1026 @arg n: Falsed I{Swiss N_LV95} or I{x_LV03} northing (C{meter}).

1027 @kwarg LV95: If C{True}, remove I{LV95} falsing, if C{False} remove

1028 I{LV03} falsing, otherwise use method C{isLV95(B{e}, B{n})}.

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

1030

1031 @return: A L{ChLVYX2Tuple}C{(Y, X)} with the unfalsed B{C{e}}

1032 respectively B{C{n}}.

1033 '''

1034 Y, X = _ChLV._falsing2(ChLV.isLV95(e, n) if LV95 is None else LV95)

1035 return ChLVYX2Tuple(e - Y, n - X, **name)

1036

1037

1038class ChLVa(_ChLV, LocalCartesian):

1039 '''Conversion between I{WGS84 geodetic} and I{Swiss} projection coordinates

1040 using the U{Approximate<https://www.SwissTopo.admin.CH/en/maps-data-online/

1041 calculation-services.html>} formulas, page 13.

1042

1043 @see: Older U{references<https://GitHub.com/alphasldiallo/Swisstopo-WGS84-LV03>}.

1044 '''

1045 def __init__(self, name=ChLV.Bern.name):

1046 '''New I{Approximate WGS84-Swiss} L{ChLVa} converter, centered at I{Bern, Ch}.

1047

1048 @kwarg name: Optional C{B{name}=Bern.name} (C{str}).

1049 '''

1050 LocalCartesian.__init__(self, latlonh0=ChLV.Bern, name=name)

1051

1052 def forward(self, latlonh, lon=None, height=0, M=None, **name):

1053 # overloaded for the _ChLV.forward.__doc__

1054 lat, lon, h, n = _llhn4(latlonh, lon, height, **name)

1055 a, b, h_ = _ChLV._llh2abh_3(lat, lon, h)

1056 a2, b2 = a**2, b**2

1057

1058 Y = fdot_(211455.93, b,

1059 -10938.51, b * a,

1060 -0.36, b * a2,

1061 -44.54, b * b2, start=72.37) # + 600_000

1062 X = fdot_(308807.95, a,

1063 3745.25, b2,

1064 76.63, a2,

1065 -194.56, b2 * a,

1066 119.79, a2 * a, start=147.07) # + 200_000

1067 return self._ChLV9Tuple(True, M, n, Y, X, h_, lat, lon, h)

1068

1069 def reverse(self, enh_, n=None, h_=0, M=None, **name): # PYCHOK signature

1070 # overloaded for the _ChLV.reverse.__doc__

1071 Y, X, h_, n = self._YXh_n4(enh_, n, h_, **name)

1072 a, b, h = _ChLV._YXh_2abh3(Y, X, h_)

1073 ab_d, a2, b2 = ChLV._ab_d, a**2, b**2

1074

1075 lat = Fdot_(3.238272, b,

1076 -0.270978, a2,

1077 -0.002528, b2,

1078 -0.0447, a2 * b,

1079 -0.014, b2 * b, start=16.9023892).fover(ab_d)

1080 lon = Fdot_(4.728982, a,

1081 0.791484, a * b,

1082 0.1306, a * b2,

1083 -0.0436, a * a2, start=2.6779094).fover(ab_d)

1084 return self._ChLV9Tuple(False, M, n, Y, X, h_, lat, lon, h)

1085

1086

1087class ChLVe(_ChLV, LocalCartesian):

1088 '''Conversion between I{WGS84 geodetic} and I{Swiss} projection coordinates

1089 using the U{Ellipsoidal approximate<https://www.SwissTopo.admin.CH/en/

1090 maps-data-online/calculation-services.html>} formulas, pp 10-11 and U{Bolliger,

1091 J.<https://eMuseum.GGGS.CH/literatur-lv/liste-Dateien/1967_Bolliger_a.pdf>}

1092 pp 148-151 (also U{GGGS<https://eMuseum.GGGS.CH/literatur-lv/liste.htm>}).

1093

1094 @note: Methods L{ChLVe.forward} and L{ChLVe.reverse} have an additional keyword

1095 argument C{B{gamma}=False} to approximate the I{meridian convergence}.

1096 If C{B{gamma}=True} a 2-tuple C{(t, gamma)} is returned with C{t} the

1097 usual result (C{ChLV9Tuple}) and C{gamma}, the I{meridian convergence}

1098 (decimal C{degrees}). To convert C{gamma} to C{grades} or C{gons}, use

1099 function L{pygeodesy.degrees2grades}.

1100

1101 @see: Older U{references<https://GitHub.com/alphasldiallo/Swisstopo-WGS84-LV03>}.

1102 '''

1103 def __init__(self, name=ChLV.Bern.name):

1104 '''New I{Approximate WGS84-Swiss} L{ChLVe} converter, centered at I{Bern, Ch}.

1105

1106 @kwarg name: Optional C{B{name}=Bern.name} (C{str}).

1107 '''

1108 LocalCartesian.__init__(self, latlonh0=ChLV.Bern, name=name)

1109

1110 def forward(self, latlonh, lon=None, height=0, M=None, gamma=False, **name): # PYCHOK gamma

1111 # overloaded for the _ChLV.forward.__doc__

1112 lat, lon, h, n = _llhn4(latlonh, lon, height, **name)

1113 a, b, h_ = _ChLV._llh2abh_3(lat, lon, h)

1114 ab_M, z, _H = ChLV._ab_M, 0, Fhorner

1115

1116 B1 = _H(a, 211428.533991, -10939.608605, -2.658213, -8.539078, -0.00345, -0.007992)

1117 B3 = _H(a, -44.232717, 4.291740, -0.309883, 0.013924)

1118 B5 = _H(a, 0.019784, -0.004277)

1119 Y = _H(b, z, B1, z, B3, z, B5).fover(ab_M) # 1,000 Km!

1120

1121 B0 = _H(a, z, 308770.746371, 75.028131, 120.435227, 0.009488, 0.070332, -0.00001)

1122 B2 = _H(a, 3745.408911, -193.792705, 4.340858, -0.376174, 0.004053)

1123 B4 = _H(a, -0.734684, 0.144466, -0.011842)

1124 X = _H(b, B0, z, B2, z, B4, z, 0.000488).fover(ab_M) # 1,000 Km!

1125

1126 t = self._ChLV9Tuple(True, M, n, Y, X, h_, lat, lon, h)

1127 if gamma:

1128 U1 = _H(a, 2255515.207166, 2642.456961, 1.284180, 2.577486, 0.001165)

1129 U3 = _H(a, -412.991934, 64.106344, -2.679566, 0.123833)

1130 U5 = _H(a, 0.204129, -0.037725)

1131 g = _H(b, z, U1, z, U3, z, U5).fover(ChLV._ab_m) # * ChLV._ab_d degrees?

1132 t = t, g

1133 return t

1134

1135 def reverse(self, enh_, n=None, h_=0, M=None, gamma=False, **name): # PYCHOK gamma

1136 # overloaded for the _ChLV.reverse.__doc__

1137 Y, X, h_, n = self._YXh_n4(enh_, n, h_, **name)

1138 a, b, h = _ChLV._YXh_2abh3(Y, X, h_)

1139 s_d, _H, z = ChLV._s_d, Fhorner, 0

1140

1141 A0 = _H(b, ChLV._sLat, 32386.4877666, -25.486822, -132.457771, 0.48747, 0.81305, -0.0069)

1142 A2 = _H(b, -2713.537919, -450.442705, -75.53194, -14.63049, -2.7604)

1143 A4 = _H(b, 24.42786, 13.20703, 4.7476)

1144 lat = _H(a, A0, z, A2, z, A4, z, -0.4249).fover(s_d)

1145

1146 A1 = _H(b, 47297.3056722, 7925.714783, 1328.129667, 255.02202, 48.17474, 9.0243)

1147 A3 = _H(b, -442.709889, -255.02202, -96.34947, -30.0808)

1148 A5 = _H(b, 9.63495, 9.0243)

1149 lon = _H(a, ChLV._sLon, A1, z, A3, z, A5).fover(s_d)

1150 # == (ChLV._sLon + a * (A1 + a**2 * (A3 + a**2 * A5))) / s_d

1151

1152 t = self._ChLV9Tuple(False, M, n, Y, X, h_, lat, lon, h)

1153 if gamma:

1154 U1 = _H(b, 106679.792202, 17876.57022, 4306.5241, 794.87772, 148.1545, 27.8725)

1155 U3 = _H(b, -1435.508, -794.8777, -296.309, -92.908)

1156 U5 = _H(b, 29.631, 27.873)

1157 g = _H(a, z, U1, z, U3, z, U5).fover(ChLV._s_ab) # degrees

1158 t = t, g

1159 return t

1160

1161

1162def _fov_2(**fov):

1163 # Half a field-of-view angle in C{degrees}.

1164 f = Degrees(Error=LocalError, **fov) * _0_5

1165 if EPS < f < _90_0:

1166 return f

1167 t = _invalid_ if f < 0 else _too_(_wide_ if f > EPS else _narrow_)

1168 raise LocalError(txt=t, **fov)

1169

1170

1171def tyr3d(tilt=INT0, yaw=INT0, roll=INT0, Vector=Vector3d, **name_Vector_kwds):

1172 '''Convert an attitude pose into a (3-D) direction vector.

1173

1174 @kwarg tilt: Pitch, elevation from horizontal (C{degrees}), negative down

1175 (clockwise rotation along and around the x-axis).

1176 @kwarg yaw: Bearing, heading (compass C{degrees360}), clockwise from North

1177 (counter-clockwise rotation along and around the z-axis).

1178 @kwarg roll: Roll, bank (C{degrees}), positive to the right and down

1179 (clockwise rotation along and around the y-axis).

1180 @kwarg Vector: Class to return the direction vector (C{Cartesian},

1181 L{Vector3d} or C{Vector3Tuple}) or C{None}.

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

1183 additional B{C{Vector}} keyword arguments, ignored if C{B{Vector}

1184 is None}.

1185

1186 @return: A named B{C{Vector}} instance or if C{B{Vector} is None},

1187 a named L{Vector3Tuple}C{(x, y, z)}.

1188

1189 @raise AttitudeError: Invalid B{C{tilt}}, B{C{yaw}} or B{C{roll}}.

1190

1191 @raise TypeError: Invalid B{C{Vector}} or B{C{name_Vector_kwds}}.

1192

1193 @see: U{Yaw, pitch, and roll rotations<http://MSL.CS.UIUC.edu/planning/node102.html>}

1194 and function L{pygeodesy.hartzell} argument C{los}, Line-Of-Sight.

1195 '''

1196 v = Attitude4Tuple(_0_0, tilt, yaw, roll).tyr3d

1197 if Vector is not type(v):

1198 n, kwds = _name2__(name_Vector_kwds, name__=tyr3d)

1199 v = Vector3Tuple(v.x, v.y, v.z, name=n) if Vector is None else \

1200 Vector(v.x, v.y, v.z, name=n, **kwds)

1201 elif name_Vector_kwds:

1202 n, _ = _name2__(name_Vector_kwds)

1203 if n:

1204 v = v.copy(name=n)

1205 return v

1206

1207

1208def _xLtp(ltp, *dflt):

1209 '''(INTERNAL) Validate B{C{ltp}} if not C{None} else B{C{dflt}}.

1210 '''

1211 if dflt and ltp is None:

1212 ltp = dflt[0]

1213 _xinstanceof(Ltp, LocalCartesian, ltp=ltp)

1214 return ltp

1215

1216# **) MIT License

1217#

1219#

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

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

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

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

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

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

1226#

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

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

1229#

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

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

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

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

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

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

1236# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pygeodesy / ltp.py: 95%

453 statements