1

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

3

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

5

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

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

8

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

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

15

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

17 _xsubclassof # .datums

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

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

20 _N_1_0 # PYCHOK used!

21# from pygeodesy.datums import _WGS84 # from .ecef

22 from pygeodesy.ecef import _EcefBase, EcefKarney, Ecef9Tuple, _llhn4, \

23 _xyzn4, _WGS84

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

25 _xkwds, _xkwds_get, _xkwds_pop2

26 from pygeodesy.fmath import fabs, fdot, Fhorner

27 from pygeodesy.fsums import _floor, _Fsumf_, fsumf_, fsum1f_

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

29 _invalid_, _lat0_, _lon0_, _name_, _too_

30# from pygeodesy.lazily import _ALL_LAZY # from vector3d

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

32 ChLVYX2Tuple, Footprint5Tuple, Local9Tuple, \

33 ChLVyx2Tuple, _XyzLocals4, _XyzLocals5, Xyz4Tuple

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

35 from pygeodesy.namedTuples import LatLon3Tuple, LatLon4Tuple, Vector3Tuple

36 from pygeodesy.props import Property, Property_RO, property_doc_, property_RO, \

37 _update_all

38 from pygeodesy.streprs import Fmt, strs, unstr

39 from pygeodesy.units import Bearing, Degrees, _isHeight, Meter

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

41 wrap180, wrap360

42 from pygeodesy.vector3d import _ALL_LAZY, Vector3d

43

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

45

46__all__ = _ALL_LAZY.ltp

47__version__ = '24.06.11'

48

49_height0_ = _height_ + _0_

50_narrow_ = 'narrow'

51_wide_ = 'wide'

52

53

54 class Attitude(_NamedBase):

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

56 '''

57 _alt = Meter( alt =_0_0)

58 _roll = Degrees(roll=_0_0)

59 _tilt = Degrees(tilt=_0_0)

60 _yaw = Bearing(yaw =_0_0)

61

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

63 '''New L{Attitude}.

64

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

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

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

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

69 (clockwise rotation along and around the x- or East axis).

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

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

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

73 (clockwise rotation along and around the y- or North axis).

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

75

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

77 B{C{roll}}.

78

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

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

81 '''

82 if _isHeight(alt_attitude):

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

84 else:

85 try:

86 t = alt_attitude.atyr

87 except AttributeError:

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

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

90 if v:

91 setattr(self, n, v)

92 n = _name__(name, _or_nameof=t)

93 if n:

94 self.name = n

95

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

97 def alt(self):

98 return self._alt

99

100 @alt.setter # PYCHOK setter!

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

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

103 if self._alt != a:

104 _update_all(self)

105 self._alt = a

106

107 altitude = alt

108

109 @Property_RO

110 def atyr(self):

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

112 '''

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

114

115 @Property_RO

116 def matrix(self):

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

118

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

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

121 '''

122 _f = fsum1f_

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

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

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

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

127 return ((ca * cb, _f(ca * sb * sg, -sa * cg), _f(ca * sb * cg, sa * sg)),

128 (sa * cb, _f(sa * sb * sg, ca * cg), _f(sa * sb * cg, -ca * sg)),

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

130

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

132 def roll(self):

133 return self._roll

134

135 @roll.setter # PYCHOK setter!

136 def roll(self, roll):

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

138 if self._roll != r:

139 _update_all(self)

140 self._roll = r

141

142 bank = roll

143

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

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

146

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

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

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

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

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

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

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

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

155 C{B{Vector} is None}.

156

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

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

159

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

161

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

163

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

165 '''

166 try:

167 try:

168 xyz = map2(float, x_xyz.xyz)

169 except AttributeError:

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

171 except (TypeError, ValueError) as x:

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

173

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

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

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

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

178

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

180 def tilt(self):

181 return self._tilt

182

183 @tilt.setter # PYCHOK setter!

184 def tilt(self, tilt):

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

186 if self._tilt != t:

187 _update_all(self)

188 self._tilt = t

189

190 elevation = pitch = tilt

191

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

193 '''Format this attitude as string.

194

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

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

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

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

199

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

201 '''

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

203

204 @Property_RO

205 def tyr3d(self):

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

207

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

209 '''

210 def _r2d(r):

211 return fsumf_(_N_1_0, *r)

212

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

214

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

216 def yaw(self):

217 return self._yaw

218

219 @yaw.setter # PYCHOK setter!

220 def yaw(self, yaw):

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

222 if self._yaw != y:

223 _update_all(self)

224 self._yaw = y

225

226 bearing = heading = yaw

227

228

229 class AttitudeError(_ValueError):

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

231 '''

232 pass

233

234

235 class Frustum(_NamedBase):

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

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

238

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

240 '''

241 _h_2 = _0_0 # half hfov in degrees

242 _ltp = None # local tangent plane

243 _tan_h_2 = _0_0 # tan(_h_2)

244 _v_2 = _0_0 # half vfov in degrees

245

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

247 '''New L{Frustum}.

248

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

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

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

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

253

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

255 '''

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

257 self._v_2 = _fov_2(vfov=vfov)

258

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

260

261 if ltp:

262 self._ltp = _xLtp(ltp)

263 if name:

264 self.name # PYCHOK effect

265

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

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

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

269

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

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

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

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

274 (clockwise rotation along and around the x- or East axis).

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

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

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

278 (clockwise rotation along and around the y- or North axis).

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

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

281 frustum's C{ltp}.

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

283

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

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

286

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

288

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

290

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

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

293 over the horizon.

294

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

296 '''

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

298 # left and right corners, or swapped

299 if r < EPS: # no roll

300 r = a * tan_h_2

301 l = -r # PYCHOK l is ell

302 else: # roll

303 r, l = tand_(r - h_2, r + h_2, roll_hfov=r) # PYCHOK l is ell

304 r *= -a # negate right positive

305 l *= -a # PYCHOK l is ell

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

307 return (l, y), (r, y)

308

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

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

311 s, c = sincos2d(b)

312 _f = fsum1f_

313 for x, y in xy5:

314 yield Xyz4Tuple(_f(x * c, y * s),

315 _f(y * c, -x * s), z, ltp)

316

317 try:

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

319 except AttributeError:

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

321

322 a = Meter(altitude=a)

323 if a < EPS: # too low

324 raise _ValueError(altitude=a)

325 if z: # PYCHOK no cover

326 z = Meter(z=z)

327 a -= z

328 if a < EPS: # z above a

329 raise _ValueError(altitude_z=a)

330 else:

331 z = _0_0

332

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

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

335 if not EPS < e < _180_0:

336 raise _ValueError(tilt=t)

337 if e > _90_0:

338 e = _loneg(e)

339 b = _umod_360(b + _180_0)

340

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

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

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

344 if fabs(y) < EPS:

345 y = _0_0

346

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

348 # center and corners, clockwise from upperleft, rolled

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

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

351 # turn center and corners by yaw, clockwise

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

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

354

355 @Property_RO

356 def hfov(self):

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

358 '''

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

360

361 @Property_RO

362 def ltp(self):

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

364 '''

365 return self._ltp

366

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

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

369

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

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

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

373

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

375 '''

376 t = self.hfov, self.vfov

377 if self.ltp:

378 t += self.ltp,

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

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

381

382 @Property_RO

383 def vfov(self):

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

385 '''

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

387

388

389 class LocalError(_ValueError):

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

391 '''

392 pass

393

394

395 class LocalCartesian(_NamedBase):

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

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

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

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

400 classGeographicLib_1_1LocalCartesian.html>}.

401

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

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

404

405 The conversions all take place via geocentric coordinates using a

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

407

408 @see: Class L{Ltp}.

409 '''

410 _ecef = EcefKarney(_WGS84)

411 _Ecef = EcefKarney

412 _lon00 = INT0 # self.lon0

413 _t0 = None # origin (..., lat0, lon0, height0, ...) L{Ecef9Tuple}

414 _9Tuple = Local9Tuple

415

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

417 '''New L{LocalCartesian} converter.

418

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

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

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

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

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

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

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

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

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

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

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

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

431 for further details.

432

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

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

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

436

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

438

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

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

441 '''

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

443

444 def __eq__(self, other):

445 '''Compare this and an other instance.

446

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

448

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

450 '''

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

452 other.ecef == self.ecef and

453 other._t0 == self._t0)

454

455 @Property_RO

456 def datum(self):

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

458 '''

459 return self.ecef.datum

460

461 @Property_RO

462 def ecef(self):

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

464 '''

465 return self._ecef

466

467 def _ecef2local(self, ecef, Xyz, name_Xyz_kwds):

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

469

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

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

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

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

474 C{B{Xyz} is None}.

475

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

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

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

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

480 always.

481

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

483 '''

484 _xinstanceof(Ecef9Tuple, ecef=ecef)

485 ltp = self

486 if ecef.datum != ltp.datum:

487 ecef = ecef.toDatum(ltp.datum)

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

489 x, y, z = self.M.rotate(ecef.xyz, *ltp._t0_xyz)

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

491 ltp, ecef, None, name=n)

492 if Xyz:

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

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

495 return r

496

497 @Property_RO

498 def ellipsoid(self):

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

500 '''

501 return self.ecef.datum.ellipsoid

502

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

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

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

506

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

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

509 @kwarg lon: Optional C{scalar} (geodetic) longitude for C{scalar}

510 B{C{latlonh}} (C{degrees}).

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

512 to and above (or below) the ellipsoid's surface.

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

514 iff available (C{bool}).

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

516

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

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

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

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

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

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

523

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

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

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

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

528 '''

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

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

531 x, y, z = self.M.rotate(t.xyz, *self._t0_xyz)

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

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

534

535 @Property_RO

536 def height0(self):

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

538 '''

539 return self._t0.height

540

541 @Property_RO

542 def lat0(self):

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

544 '''

545 return self._t0.lat

546

547 @Property_RO

548 def latlonheight0(self):

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

550 '''

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

552

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

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

555

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

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

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

559

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

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

562 optionally including rotation matrix C{M} otherwise C{None}.

563 '''

564 _xinstanceof(*_XyzLocals5, local=local)

565 t = self.M.unrotate(local.xyz, *self._t0_xyz)

566 if nine:

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

568 return t

569

570 @Property_RO

571 def lon0(self):

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

573 '''

574 return self._t0.lon

575

576 @Property

577 def lon00(self):

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

579 '''

580 return self._lon00

581

582 @lon00.setter # PYCHOK setter!

583 def lon00(self, lon00):

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

585 '''

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

587 self._lon00 = Degrees(lon00=lon00)

588

589 @Property_RO

590 def M(self):

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

592 '''

593 return self._t0.M

594

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

596 '''Reset this converter, see L{LocalCartesian.__init__} for more details.

597 '''

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

599 if isinstance(latlonh0, LocalCartesian):

600 if self._t0:

601 _update_all(self)

602 self._ecef = latlonh0.ecef

603 self._lon00 = latlonh0.lon00

604 self._t0 = latlonh0._t0

605 n = _name__(name, _or_nameof=latlonh0)

606 else:

607 n = _name__(name, _or_nameof=self)

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

609 Error=LocalError, name=n)

610 if ecef: # PYCHOK no cover

611 _xinstanceof(self._Ecef, ecef=ecef)

612 _update_all(self)

613 self._ecef = ecef

614 elif self._t0:

615 _update_all(self)

616 self._t0 = self.ecef._forward(lat0, lon0, height0, n, M=True)

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

618 if n:

619 self.rename(n)

620

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

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

623

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

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

626 @kwarg y: Local C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).

627 @kwarg z: Local C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).

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

629 available (C{bool}).

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

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

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

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

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

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

636

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

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

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

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

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

642

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

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

645 '''

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

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

648 c = self.M.unrotate((x, y, z), *self._t0_xyz)

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

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

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

652

653 @Property_RO

654 def _t0_xyz(self):

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

656 '''

657 return self._t0.xyz

658

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

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

661

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

663

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

665 '''

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

667

668

669 class Ltp(LocalCartesian):

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

671 (re-)configurable ECEF converter.

672 '''

673 _Ecef = _EcefBase

674

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

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

677

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

679 L{EcefFarrell22}, L{EcefSudano}, L{EcefVeness} or

680 L{EcefYou} I{instance}), overriding the default

681 L{EcefKarney}C{(datum=Datums.WGS84)} for C{scalar}.

682

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

684

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

686 '''

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

688 ecef=ecef, **lon00_name)

689

690 @Property

691 def ecef(self):

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

693 '''

694 return self._ecef

695

696 @ecef.setter # PYCHOK setter!

697 def ecef(self, ecef):

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

699

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

701 '''

702 _xinstanceof(_EcefBase, ecef=ecef)

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

704 self.reset(self._t0)

705 self._ecef = ecef

706

707

708 class _ChLV(object):

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

710 '''

711 _03_falsing = ChLVyx2Tuple(0.6e6, 0.2e6)

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

713 _95_falsing = ChLVEN2Tuple(2.6e6, 1.2e6)

714

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

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

717 '''

718 if bool(M): # PYCHOK no cover

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

720 n = _DOT_(self.__class__.__name__, m.__name__)

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

722 t = Y_X_h_lat_lon_h + (self, self._t0, None) # PYCHOK _t0

723 return ChLV9Tuple(t, name=name)

724

725 @property_RO

726 def _enh_n_h(self):

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

728 '''

729 _ChLV._enh_n_h = t = _args_kwds_names(_ChLV.reverse)[1:4] # overwrite property_RO

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

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

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

733 return t

734

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

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

737

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

739 @kwarg lon: Optional, C{scalar} (geodetic) longitude for C{scalar} B{C{latlonh}} (C{degrees}).

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

741 (C{meter}) for C{scalar} B{C{latlonh}} and B{C{lon}}.

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

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

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

745

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

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

748 C{lon} and C{height}, this C{ChLV*} instance and C{ecef} (L{Ecef9Tuple}) at

749 I{Bern, Ch} and rotation matrix C{M}. The returned C{ltp} is this C{ChLV},

750 C{ChLVa} or C{ChLVe} instance.

751

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

753 '''

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

755

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

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

758

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

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

761 @kwarg n: Falsed I{Swiss N_LV85} or I{x_LV03} northing for C{scalar} B{C{enh_}} and

762 B{C{h_}} (C{meter}).

763 @kwarg h_: I{Swiss h'} height for C{scalar} B{C{enh_}} and B{C{n}} (C{meter}).

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

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

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

767

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

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

770 C{lon} and C{height}, this C{ChLV*} instance and C{ecef} (L{Ecef9Tuple}) at

771 I{Bern, Ch} and rotation matrix C{M}. The returned C{ltp} is this C{ChLV},

772 C{ChLVa} or C{ChLVe} instance.

773

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

775 '''

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

777

778 @staticmethod

779 def _falsing2(LV95):

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

781 '''

782 return _ChLV._95_falsing if LV95 in (True, 95) else (

783 _ChLV._03_falsing if LV95 in (False, 3) else ChLVYX2Tuple(0, 0))

784

785 @staticmethod

786 def _llh2abh_3(lat, lon, h):

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

788 '''

789 def _deg2ab(deg, sLL):

790 # convert degrees to arc-seconds

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

792 d = _floor(ds)

793 t = (ds - d) * p

794 m = _floor(t)

795 s = (t - m) * p

796 if swap:

797 d, s = s, d

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

799

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

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

802 return (s - sLL) / ChLV._s_ab

803

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

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

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

807 return a, b, h_

808

809 @staticmethod

810 def _YXh_2abh3(Y, X, h_):

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

812 '''

813 def _YX2ab(YX):

814 return YX * ChLV._ab_m

815

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

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

818 return a, b, h

819

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

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

822 '''

823 Y, X, h_, name = _xyzn4(enh_, n, h_, ChLV9Tuple,

824 _xyz_y_z_names=self._enh_n_h, **name)

825 if isinstance(enh_, ChLV9Tuple):

826 Y, X = enh_.Y, enh_.X

827 else: # isscalar(enh_)

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

829 return Y, X, h_, name

830

831

832 class ChLV(_ChLV, Ltp):

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

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

835

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

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

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

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

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

841 '''

842 _9Tuple = ChLV9Tuple

843

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

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

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

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

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

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

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

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

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

853

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

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

856

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

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

859

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

861 '''

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

863

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

865 # overloaded for the _ChLV.forward.__doc__

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

867

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

869 # overloaded for the _ChLV.reverse.__doc__

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

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

872

873 @staticmethod

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

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

876

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

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

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

880 C{LV03} falsing, otherwise leave unfalsed.

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

882

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

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

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

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

887 '''

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

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

890

891 @staticmethod

892 def isLV03(e, n):

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

894

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

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

897

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

899 LV03}, projection C{False} otherwise.

900 '''

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

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

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

904

905 @staticmethod

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

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

908

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

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

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

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

913

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

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

916 '''

917 if ChLV.isLV03(e, n):

918 return False

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

920 return True

921 elif raiser: # PYCHOK no cover

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

923 return None

924

925 @staticmethod

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

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

928

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

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

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

932 I{LV03} falsing, otherwise use method C{isLV95(B{e},

933 B{n})}.

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

935

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

937 respectively B{C{n}}.

938 '''

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

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

941

942

943 class ChLVa(_ChLV, LocalCartesian):

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

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

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

947

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

949 '''

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

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

952

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

954 '''

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

956

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

958 # overloaded for the _ChLV.forward.__doc__

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

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

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

962

963 Y = fsumf_( 72.37, 211455.93 * b,

964 -10938.51 * b * a,

965 -0.36 * b * a2,

966 -44.54 * b * b2) # + 600_000

967 X = fsumf_(147.07, 308807.95 * a,

968 3745.25 * b2,

969 76.63 * a2,

970 -194.56 * b2 * a,

971 119.79 * a2 * a) # + 200_000

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

973

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

975 # overloaded for the _ChLV.reverse.__doc__

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

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

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

979

980 lat = _Fsumf_(16.9023892, 3.238272 * b,

981 -0.270978 * a2,

982 -0.002528 * b2,

983 -0.0447 * a2 * b,

984 -0.014 * b2 * b).fover(ab_d)

985 lon = _Fsumf_( 2.6779094, 4.728982 * a,

986 0.791484 * a * b,

987 0.1306 * a * b2,

988 -0.0436 * a * a2).fover(ab_d)

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

990

991

992 class ChLVe(_ChLV, LocalCartesian):

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

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

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

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

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

998

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

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

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

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

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

1004 use function L{pygeodesy.degrees2grades}.

1005

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

1007 '''

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

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

1010

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

1012 '''

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

1014

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

1016 # overloaded for the _ChLV.forward.__doc__

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

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

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

1020

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

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

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

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

1025

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

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

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

1029 B6 = 0.000488

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

1031

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

1033 if gamma:

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

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

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

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

1038 t = t, g

1039 return t

1040

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

1042 # overloaded for the _ChLV.reverse.__doc__

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

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

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

1046

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

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

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

1050 A6 = -0.4249

1051 lat = _H(a, A0, z, A2, z, A4, z, A6).fover(s_d)

1052

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

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

1055 A5 = _H(b, 9.63495, 9.0243)

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

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

1058

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

1060 if gamma:

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

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

1063 U5 = _H(b, 29.631, 27.873)

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

1065 t = t, g

1066 return t

1067

1068

1069 def _fov_2(**fov):

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

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

1072 if EPS < f < _90_0:

1073 return f

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

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

1076

1077

1078 def _toLocal(inst, ltp, Xyz, Xyz_kwds):

1079 '''(INTENRAL) Helper for C{CartesianBase.toLocal} and C{latLonBase.toLocal}.

1080 '''

1081 return _xLtp(ltp, inst._Ltp)._ecef2local(inst._ecef9, Xyz, Xyz_kwds)

1082

1083

1084 def _toLtp(inst, Ecef, ecef9, name):

1085 '''(INTENRAL) Helper for C{CartesianBase.toLtp} and C{latLonBase.toLtp}.

1086 '''

1087 return inst._Ltp if Ecef in (None, inst.Ecef) and not name else \

1088 Ltp(ecef9, ecef=Ecef(inst.datum), name=inst._name__(name))

1089

1090

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

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

1093

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

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

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

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

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

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

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

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

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

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

1104 C{B{Vector} is None}.

1105

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

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

1108

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

1110

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

1112

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

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

1115 '''

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

1117 if Vector is not type(v):

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

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

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

1121 elif name_Vector_kwds:

1122 n, _ = _name2__(name_Vector_kwds)

1123 if n:

1124 v = v.copy(name=n)

1125 return v

1126

1127

1128 def _xLtp(ltp, *dflt):

1129 '''(INTERNAL) Validate B{C{ltp}} or ist B{C{dflt}}.

1130 '''

1131 if dflt and ltp is None:

1132 ltp = dflt[0]

1133 _xinstanceof(Ltp, LocalCartesian, ltp=ltp)

1134 return ltp

1135

1136# **) MIT License

1137#

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

1139#

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

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

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

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

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

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

1146#

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

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

1149#

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

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

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

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

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

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

1156# OTHER DEALINGS IN THE SOFTWARE.