1

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

3

4 u'''(INTERNAL) Private elliposiodal and spherical C{Nvector} base classes

5 L{LatLonNvectorBase} and L{NvectorBase} and function L{sumOf}.

6

7 Pure Python implementation of C{n-vector}-based geodesy tools for ellipsoidal

8 earth models, transcoded from JavaScript originals by I{(C) Chris Veness 2005-2016}

9 and published under the same MIT Licence**, see U{Vector-based geodesy

10<https://www.Movable-Type.co.UK/scripts/latlong-vectors.html>}.

11'''

12

13# from pygeodesy.basics import map1 # from .namedTuples

14 from pygeodesy.constants import EPS, EPS1, EPS_2, R_M, _2_0, _N_2_0

15# from pygeodesy.datums import _spherical_datum # from .formy

16 from pygeodesy.errors import IntersectionError, _ValueError, VectorError, \

17 _xkwds, _xkwds_pop

18 from pygeodesy.fmath import fdot, fidw, hypot_ # PYCHOK fdot shared

19 from pygeodesy.fsums import Fsum, fsumf_

20 from pygeodesy.formy import _isequalTo, n_xyz2latlon, n_xyz2philam, \

21 _spherical_datum

22 from pygeodesy.interns import NN, _1_, _2_, _3_, _bearing_, _coincident_, \

23 _COMMASPACE_, _distance_, _h_, _insufficient_, \

24 _intersection_, _no_, _NorthPole_, _point_, \

25 _pole_, _SPACE_, _SouthPole_, _under

26 from pygeodesy.latlonBase import LatLonBase, _ALL_DOCS, _MODS

27# from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS # from .latlonBase

28 from pygeodesy.named import notImplemented, _xother3

29 from pygeodesy.namedTuples import Trilaterate5Tuple, Vector3Tuple, \

30 Vector4Tuple, map1

31 from pygeodesy.props import deprecated_method, Property_RO, property_doc_, \

32 property_RO, _update_all

33 from pygeodesy.streprs import Fmt, hstr, unstr, _xattrs

34 from pygeodesy.units import Bearing, Height, Radius_, Scalar

35 from pygeodesy.utily import sincos2d, _unrollon, _unrollon3

36 from pygeodesy.vector3d import Vector3d, _xyzhdn3

37

38 from math import fabs, sqrt

39

40__all__ = (_NorthPole_, _SouthPole_) # constants

41__version__ = '23.10.24'

42

43

44 class NvectorBase(Vector3d): # XXX kept private

45 '''Base class for ellipsoidal and spherical C{Nvector}s.

46 '''

47 _datum = None # L{Datum}, overriden

48 _h = Height(h=0) # height (C{meter})

49 _H = NN # height prefix (C{str}), '↑' in JS version

50

51 def __init__(self, x_xyz, y=None, z=None, h=0, ll=None, datum=None, name=NN):

52 '''New n-vector normal to the earth's surface.

53

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

55 (C{Nvector}, L{Vector3d}, L{Vector3Tuple} or

56 L{Vector4Tuple}).

57 @kwarg y: Y component of vector (C{scalar}), ignored if B{C{x_xyz}}

58 is not C{scalar}, otherwise same units as B{C{x_xyz}}.

59 @kwarg z: Z component of vector (C{scalar}), ignored if B{C{x_xyz}}

60 is not C{scalar}, otherwise same units as B{C{x_xyz}}.

61 @kwarg h: Optional height above surface (C{meter}).

62 @kwarg ll: Optional, original latlon (C{LatLon}).

63 @kwarg datum: Optional, I{pass-thru} datum (L{Datum}).

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

65

66 @raise TypeError: Non-scalar B{C{x}}, B{C{y}} or B{C{z}}

67 coordinate or B{C{x}} not an C{Nvector},

68 L{Vector3Tuple} or L{Vector4Tuple} or

69 invalid B{C{datum}}.

70

71 @example:

72

73 >>> from pygeodesy.sphericalNvector import Nvector

74 >>> v = Nvector(0.5, 0.5, 0.7071, 1)

75 >>> v.toLatLon() # 45.0°N, 045.0°E, +1.00m

76 '''

77 h, d, n = _xyzhdn3(x_xyz, h, datum, ll)

78 Vector3d.__init__(self, x_xyz, y=y, z=z, ll=ll, name=name or n)

79 if h:

80 self.h = h

81 if d is not None:

82 self._datum = _spherical_datum(d, name=self.name) # pass-thru

83

84 @Property_RO

85 def datum(self):

86 '''Get the I{pass-thru} datum (C{Datum}) or C{None}.

87 '''

88 return self._datum

89

90 @Property_RO

91 def Ecef(self):

92 '''Get the ECEF I{class} (L{EcefKarney}), I{lazily}.

93 '''

94 return _MODS.ecef.EcefKarney # default

95

96 @property_RO

97 def ellipsoidalNvector(self):

98 '''Get the C{Nvector type} iff ellipsoidal, overloaded in L{ellipsoidalNvector.Nvector}.

99 '''

100 return False

101

102 @property_doc_(''' the height above surface (C{meter}).''')

103 def h(self):

104 '''Get the height above surface (C{meter}).

105 '''

106 return self._h

107

108 @h.setter # PYCHOK setter!

109 def h(self, h):

110 '''Set the height above surface (C{meter}).

111

112 @raise TypeError: If B{C{h}} invalid.

113

114 @raise VectorError: If B{C{h}} invalid.

115 '''

116 h = Height(h=h, Error=VectorError)

117 if self._h != h:

118 _update_all(self)

119 self._h = h

120

121 @property_doc_(''' the height prefix (C{str}).''')

122 def H(self):

123 '''Get the height prefix (C{str}).

124 '''

125 return self._H

126

127 @H.setter # PYCHOK setter!

128 def H(self, H):

129 '''Set the height prefix (C{str}).

130 '''

131 self._H = str(H) if H else NN

132

133 def hStr(self, prec=-2, m=NN):

134 '''Return a string for the height B{C{h}}.

135

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

137 @kwarg m: Optional unit of the height (C{str}).

138

139 @see: Function L{pygeodesy.hstr}.

140 '''

141 return NN(self.H, hstr(self.h, prec=prec, m=m))

142

143 @Property_RO

144 def isEllipsoidal(self):

145 '''Check whether this n-vector is ellipsoidal (C{bool} or C{None} if unknown).

146 '''

147 return self.datum.isEllipsoidal if self.datum else None

148

149 @Property_RO

150 def isSpherical(self):

151 '''Check whether this n-vector is spherical (C{bool} or C{None} if unknown).

152 '''

153 return self.datum.isSpherical if self.datum else None

154

155 @Property_RO

156 def lam(self):

157 '''Get the (geodetic) longitude in C{radians} (C{float}).

158 '''

159 return self.philam.lam

160

161 @Property_RO

162 def lat(self):

163 '''Get the (geodetic) latitude in C{degrees} (C{float}).

164 '''

165 return self.latlon.lat

166

167 @Property_RO

168 def latlon(self):

169 '''Get the (geodetic) lat-, longitude in C{degrees} (L{LatLon2Tuple}C{(lat, lon)}).

170 '''

171 return n_xyz2latlon(self.x, self.y, self.z, name=self.name)

172

173 @Property_RO

174 def latlonheight(self):

175 '''Get the (geodetic) lat-, longitude in C{degrees} and height (L{LatLon3Tuple}C{(lat, lon, height)}).

176 '''

177 return self.latlon.to3Tuple(self.h)

178

179 @Property_RO

180 def latlonheightdatum(self):

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

182 '''

183 return self.latlonheight.to4Tuple(self.datum)

184

185 @Property_RO

186 def lon(self):

187 '''Get the (geodetic) longitude in C{degrees} (C{float}).

188 '''

189 return self.latlon.lon

190

191 @Property_RO

192 def phi(self):

193 '''Get the (geodetic) latitude in C{radians} (C{float}).

194 '''

195 return self.philam.phi

196

197 @Property_RO

198 def philam(self):

199 '''Get the (geodetic) lat-, longitude in C{radians} (L{PhiLam2Tuple}C{(phi, lam)}).

200 '''

201 return n_xyz2philam(self.x, self.y, self.z, name=self.name)

202

203 @Property_RO

204 def philamheight(self):

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

206 '''

207 return self.philam.to3Tuple(self.h)

208

209 @Property_RO

210 def philamheightdatum(self):

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

212 '''

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

214

215 @property_RO

216 def sphericalNvector(self):

217 '''Get the C{Nvector type} iff spherical, overloaded in L{sphericalNvector.Nvector}.

218 '''

219 return False

220

221 @deprecated_method

222 def to2ab(self): # PYCHOK no cover

223 '''DEPRECATED, use property L{philam}.

224

225 @return: A L{PhiLam2Tuple}C{(phi, lam)}.

226 '''

227 return self.philam

228

229 @deprecated_method

230 def to3abh(self, height=None): # PYCHOK no cover

231 '''DEPRECATED, use property L{philamheight} or C{philam.to3Tuple(B{height})}.

232

233 @kwarg height: Optional height, overriding this

234 n-vector's height (C{meter}).

235

236 @return: A L{PhiLam3Tuple}C{(phi, lam, height)}.

237

238 @raise ValueError: Invalid B{C{height}}.

239 '''

240 return self.philamheight if height in (None, self.h) else \

241 self.philam.to3Tuple(height)

242

243 def toCartesian(self, h=None, Cartesian=None, datum=None, **Cartesian_kwds):

244 '''Convert this n-vector to C{Nvector}-based cartesian (ECEF) coordinates.

245

246 @kwarg h: Optional height, overriding this n-vector's height (C{meter}).

247 @kwarg Cartesian: Optional class to return the (ECEF) coordinates

248 (C{Cartesian}).

249 @kwarg datum: Optional datum (C{Datum}), overriding this datum.

250 @kwarg Cartesian_kwds: Optional, additional B{C{Cartesian}} keyword

251 arguments, ignored if C{B{Cartesian} is None}.

252

253 @return: The cartesian (ECEF) coordinates (B{C{Cartesian}}) or

254 if C{B{Cartesian} is None}, an L{Ecef9Tuple}C{(x, y, z,

255 lat, lon, height, C, M, datum)} with C{C} and C{M} if

256 available.

257

258 @raise TypeError: Invalid B{C{Cartesian}} or B{C{Cartesian_kwds}}

259 argument.

260

261 @raise ValueError: Invalid B{C{h}}.

262

263 @example:

264

265 >>> v = Nvector(0.5, 0.5, 0.7071)

266 >>> c = v.toCartesian() # [3194434, 3194434, 4487327]

267 >>> p = c.toLatLon() # 45.0°N, 45.0°E

268 '''

269 D = _spherical_datum(datum or self.datum, name=self.name)

270 E = D.ellipsoid

271 h = self.h if h is None else Height(h)

272

273 x, y, z = self.x, self.y, self.z

274 # Kenneth Gade eqn 22

275 n = E.b / hypot_(x * E.a_b, y * E.a_b, z)

276 r = h + n * E.a2_b2

277

278 x *= r

279 y *= r

280 z *= h + n

281

282 if Cartesian is None:

283 r = self.Ecef(D).reverse(x, y, z, M=True)

284 else:

285 kwds = _xkwds(Cartesian_kwds, datum=D) # h=0

286 r = Cartesian(x, y, z, **kwds)

287 return self._xnamed(r)

288

289 @deprecated_method

290 def to2ll(self): # PYCHOK no cover

291 '''DEPRECATED, use property L{latlon}.

292

293 @return: A L{LatLon2Tuple}C{(lat, lon)}.

294 '''

295 return self.latlon

296

297 @deprecated_method

298 def to3llh(self, height=None): # PYCHOK no cover

299 '''DEPRECATED, use property C{latlonheight} or C{latlon.to3Tuple(B{height})}.

300

301 @kwarg height: Optional height, overriding this

302 n-vector's height (C{meter}).

303

304 @return: A L{LatLon3Tuple}C{(lat, lon, height)}.

305

306 @raise ValueError: Invalid B{C{height}}.

307 '''

308 return self.latlonheight if height in (None, self.h) else \

309 self.latlon.to3Tuple(height)

310

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

312 '''Convert this n-vector to an C{Nvector}-based geodetic point.

313

314 @kwarg height: Optional height, overriding this n-vector's

315 height (C{meter}).

316 @kwarg LatLon: Optional class to return the geodetic point

317 (C{LatLon}) or C{None}.

318 @kwarg datum: Optional, spherical datum (C{Datum}).

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

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

321

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

323 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M,

324 datum)} with C{C} and C{M} if available.

325

326 @raise TypeError: Invalid B{C{LatLon}} or B{C{LatLon_kwds}}

327 argument.

328

329 @raise ValueError: Invalid B{C{height}}.

330

331 @example:

332

333 >>> v = Nvector(0.5, 0.5, 0.7071)

334 >>> p = v.toLatLon() # 45.0°N, 45.0°E

335 '''

336 d = _spherical_datum(datum or self.datum, name=self.name)

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

338 # use self.Cartesian(Cartesian=None) for better accuracy of the height

339 # than self.Ecef(d).forward(self.lat, self.lon, height=h, M=True)

340 if LatLon is None:

341 r = self.toCartesian(h=h, Cartesian=None, datum=d)

342 else:

343 kwds = _xkwds(LatLon_kwds, height=h, datum=d)

344 r = self._xnamed(LatLon(self.lat, self.lon, **kwds))

345 return r

346

347 def toStr(self, prec=5, fmt=Fmt.PAREN, sep=_COMMASPACE_): # PYCHOK expected

348 '''Return a string representation of this n-vector.

349

350 Height component is only included if non-zero.

351

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

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

354 @kwarg sep: Optional separator between components (C{str}).

355

356 @return: Comma-separated C{"(x, y, z [, h])"} enclosed in

357 B{C{fmt}} brackets (C{str}).

358

359 @example:

360

361 >>> Nvector(0.5, 0.5, 0.7071).toStr() # (0.5, 0.5, 0.7071)

362 >>> Nvector(0.5, 0.5, 0.7071, 1).toStr(-3) # (0.500, 0.500, 0.707, +1.00)

363 '''

364 t = Vector3d.toStr(self, prec=prec, fmt=NN, sep=sep)

365 if self.h:

366 t = sep.join((t, self.hStr()))

367 return (fmt % (t,)) if fmt else t

368

369 def toVector3d(self, norm=True):

370 '''Convert this n-vector to a 3-D vector, I{ignoring

371 the height}.

372

373 @kwarg norm: Normalize the 3-D vector (C{bool}).

374

375 @return: The (normalized) vector (L{Vector3d}).

376 '''

377 v = Vector3d.unit(self) if norm else self

378 return Vector3d(v.x, v.y, v.z, name=self.name)

379

380 @deprecated_method

381 def to4xyzh(self, h=None): # PYCHOK no cover

382 '''DEPRECATED, use property L{xyzh} or C{xyz.to4Tuple(B{h})}.

383 '''

384 return self.xyzh if h in (None, self.h) else Vector4Tuple(

385 self.x, self.y, self.z, h, name=self.name)

386

387 def unit(self, ll=None):

388 '''Normalize this n-vector to unit length.

389

390 @kwarg ll: Optional, original latlon (C{LatLon}).

391

392 @return: Normalized vector (C{Nvector}).

393 '''

394 return _xattrs(Vector3d.unit(self, ll=ll), _under(_h_))

395

396 @Property_RO

397 def xyzh(self):

398 '''Get this n-vector's components (L{Vector4Tuple}C{(x, y, z, h)})

399 '''

400 return self.xyz.to4Tuple(self.h)

401

402

403 NorthPole = NvectorBase(0, 0, +1, name=_NorthPole_) # North pole (C{Nvector})

404 SouthPole = NvectorBase(0, 0, -1, name=_SouthPole_) # South pole (C{Nvector})

405

406

407 class _N_vector_(NvectorBase):

408 '''(INTERNAL) Minimal, low-overhead C{n-vector}.

409 '''

410 def __init__(self, x, y, z, h=0, name=NN):

411 self._x, self._y, self._z = x, y, z

412 if h:

413 self._h = h

414 if name:

415 self.name = name

416

417

418 class LatLonNvectorBase(LatLonBase):

419 '''(INTERNAL) Base class for n-vector-based ellipsoidal

420 and spherical C{LatLon} classes.

421 '''

422

423 def _update(self, updated, *attrs, **setters): # PYCHOK _Nv=None

424 '''(INTERNAL) Zap cached attributes if updated.

425

426 @see: C{ellipsoidalNvector.LatLon} and C{sphericalNvector.LatLon}

427 for the special case of B{C{_Nv}}.

428 '''

429 if updated:

430 _Nv = _xkwds_pop(setters, _Nv=None)

431 if _Nv is not None:

432 if _Nv._fromll is not None:

433 _Nv._fromll = None

434 self._Nv = None

435 LatLonBase._update(self, updated, *attrs, **setters)

436

437# def distanceTo(self, other, **kwds): # PYCHOK no cover

438# '''I{Must be overloaded}.'''

439# _MODS.named.notOverloaded(self, other, **kwds)

440

441 def intersections2(self, radius1, other, radius2, **kwds): # PYCHOK expected

442 '''B{Not implemented}, throws a C{NotImplementedError} always.'''

443 notImplemented(self, radius1, other, radius2, **kwds)

444

445 def others(self, *other, **name_other_up):

446 '''Refined class comparison.

447

448 @arg other: The other instance (C{LatLonNvectorBase}).

449 @kwarg name_other_up: Overriding C{name=other} and C{up=1}

450 keyword arguments.

451

452 @return: The B{C{other}} if compatible.

453

454 @raise TypeError: Incompatible B{C{other}} C{type}.

455 '''

456 if other:

457 other0 = other[0]

458 if isinstance(other0, (self.__class__, LatLonNvectorBase)): # XXX NvectorBase?

459 return other0

460

461 other, name, up = _xother3(self, other, **name_other_up)

462 if not isinstance(other, (self.__class__, LatLonNvectorBase)): # XXX NvectorBase?

463 LatLonBase.others(self, other, name=name, up=up + 1)

464 return other

465

466 def toNvector(self, Nvector=NvectorBase, **Nvector_kwds): # PYCHOK signature

467 '''Convert this point to C{Nvector} components, I{including height}.

468

469 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}} keyword

470 arguments, ignored if C{B{Nvector} is None}.

471

472 @return: An B{C{Nvector}} or a L{Vector4Tuple}C{(x, y, z, h)} if

473 B{C{Nvector}} is C{None}.

474

475 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}

476 argument.

477 '''

478 return LatLonBase.toNvector(self, Nvector=Nvector, **Nvector_kwds)

479

480 def triangulate(self, bearing1, other, bearing2, height=None, wrap=False):

481 '''Locate a point given this and an other point and a bearing

482 at this and the other point.

483

484 @arg bearing1: Bearing at this point (compass C{degrees360}).

485 @arg other: The other point (C{LatLon}).

486 @arg bearing2: Bearing at the other point (compass C{degrees360}).

487 @kwarg height: Optional height at the triangulated point,

488 overriding the mean height (C{meter}).

489 @kwarg wrap: If C{True}, use this and the B{C{other}} point

490 I{normalized} (C{bool}).

491

492 @return: Triangulated point (C{LatLon}).

493

494 @raise TypeError: Invalid B{C{other}} point.

495

496 @raise Valuerror: Points coincide.

497

498 @example:

499

500 >>> p = LatLon("47°18.228'N","002°34.326'W") # Basse Castouillet

501 >>> q = LatLon("47°18.664'N","002°31.717'W") # Basse Hergo

502 >>> t = p.triangulate(7, q, 295) # 47.323667°N, 002.568501°W'

503 '''

504 return _triangulate(self, bearing1, self.others(other), bearing2,

505 height=height, wrap=wrap, LatLon=self.classof)

506

507 def trilaterate(self, distance1, point2, distance2, point3, distance3,

508 radius=R_M, height=None, useZ=False, wrap=False):

509 '''Locate a point at given distances from this and two other points.

510

511 @arg distance1: Distance to this point (C{meter}, same units

512 as B{C{radius}}).

513 @arg point2: Second reference point (C{LatLon}).

514 @arg distance2: Distance to point2 (C{meter}, same units as

515 B{C{radius}}).

516 @arg point3: Third reference point (C{LatLon}).

517 @arg distance3: Distance to point3 (C{meter}, same units as

518 B{C{radius}}).

519 @kwarg radius: Mean earth radius (C{meter}).

520 @kwarg height: Optional height at trilaterated point, overriding

521 the mean height (C{meter}, same units as B{C{radius}}).

522 @kwarg useZ: Include Z component iff non-NaN, non-zero (C{bool}).

523 @kwarg wrap: If C{True}, use this, B{C{point2}} and B{C{point3}}

524 I{normalized} (C{bool}).

525

526 @return: Trilaterated point (C{LatLon}).

527

528 @raise IntersectionError: No intersection, trilateration failed.

529

530 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.

531

532 @raise ValueError: Some B{C{points}} coincide or invalid B{C{distance1}},

533 B{C{distance2}}, B{C{distance3}} or B{C{radius}}.

534

535 @see: U{Trilateration<https://WikiPedia.org/wiki/Trilateration>},

536 Veness' JavaScript U{Trilateration<https://www.Movable-Type.co.UK/

537 scripts/latlong-vectors.html>} and method C{LatLon.trilaterate5}

538 of other, non-C{Nvector LatLon} classes.

539 '''

540 return _trilaterate(self, distance1, self.others(point2=point2), distance2,

541 self.others(point3=point3), distance3,

542 radius=radius, height=height, useZ=useZ,

543 wrap=wrap, LatLon=self.classof)

544

545 def trilaterate5(self, distance1, point2, distance2, point3, distance3, # PYCHOK signature

546 area=False, eps=EPS1, radius=R_M, wrap=False):

547 '''B{Not implemented} for C{B{area}=True} and falls back to method

548 C{trilaterate} otherwise.

549

550 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)}

551 with a single trilaterated intersection C{minPoint I{is}

552 maxPoint}, C{min I{is} max} the nearest intersection

553 margin and count C{n = 1}.

554

555 @raise NotImplementedError: Keyword argument C{B{area}=True} not

556 (yet) supported.

557

558 @see: Method L{trilaterate} for other and more details.

559 '''

560 if area:

561 notImplemented(self, area=area)

562

563 t = _trilaterate(self, distance1, self.others(point2=point2), distance2,

564 self.others(point3=point3), distance3,

565 radius=radius, useZ=True, wrap=wrap,

566 LatLon=self.classof)

567 # ... and handle B{C{eps}} and C{IntersectionError}

568 # like function C{.latlonBase._trilaterate5}

569 d = self.distanceTo(t, radius=radius, wrap=wrap) # PYCHOK distanceTo

570 d = min(fabs(distance1 - d), fabs(distance2 - d), fabs(distance3 - d))

571 if d < eps: # min is max, minPoint is maxPoint

572 return Trilaterate5Tuple(d, t, d, t, 1) # n = 1

573 t = _SPACE_(_no_(_intersection_), Fmt.PAREN(min.__name__, Fmt.f(d, prec=3)))

574 raise IntersectionError(area=area, eps=eps, radius=radius, wrap=wrap, txt=t)

575

576

577 def _nsumOf(nvs, h_None, Vector, Vector_kwds): # .sphericalNvector, .vector3d

578 '''(INTERNAL) Separated to allow callers to embellish exceptions.

579 '''

580 X, Y, Z, n = Fsum(), Fsum(), Fsum(), 0

581 H = Fsum() if h_None is None else n

582 for n, v in enumerate(nvs or ()): # one pass

583 X += v.x

584 Y += v.y

585 Z += v.z

586 H += v.h

587 if n < 1:

588 raise ValueError(_SPACE_(Fmt.PARENSPACED(len=n), _insufficient_))

589

590 x, y, z = map1(float, X, Y, Z)

591 h = H.fover(n) if h_None is None else h_None

592 return Vector3Tuple(x, y, z).to4Tuple(h) if Vector is None else \

593 Vector(x, y, z, **_xkwds(Vector_kwds, h=h))

594

595

596 def sumOf(nvectors, Vector=None, h=None, **Vector_kwds):

597 '''Return the I{vectorial} sum of two or more n-vectors.

598

599 @arg nvectors: Vectors to be added (C{Nvector}[]).

600 @kwarg Vector: Optional class for the vectorial sum (C{Nvector})

601 or C{None}.

602 @kwarg h: Optional height, overriding the mean height (C{meter}).

603 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword

604 arguments, ignored if C{B{Vector} is None}.

605

606 @return: Vectorial sum (B{C{Vector}}) or a L{Vector4Tuple}C{(x, y,

607 z, h)} if B{C{Vector}} is C{None}.

608

609 @raise VectorError: No B{C{nvectors}}.

610 '''

611 try:

612 return _nsumOf(nvectors, h, Vector, Vector_kwds)

613 except (TypeError, ValueError) as x:

614 raise VectorError(nvectors=nvectors, Vector=Vector, cause=x)

615

616

617 def _triangulate(point1, bearing1, point2, bearing2, height=None,

618 wrap=False, **LatLon_and_kwds):

619 # (INTERNAL) Locate a point given two known points and initial

620 # bearings from those points, see C{LatLon.triangulate} above

621

622 def _gc(p, b, _i_):

623 n = p.toNvector()

624 de = NorthPole.cross(n, raiser=_pole_).unit() # east vector @ n

625 dn = n.cross(de) # north vector @ n

626 s, c = sincos2d(Bearing(b, name=_bearing_ + _i_))

627 dest = de.times(s)

628 dnct = dn.times(c)

629 d = dnct.plus(dest) # direction vector @ n

630 return n.cross(d) # great circle point + bearing

631

632 if wrap:

633 point2 = _unrollon(point1, point2, wrap=wrap)

634 if _isequalTo(point1, point2, eps=EPS):

635 raise _ValueError(points=point2, wrap=wrap, txt=_coincident_)

636

637 gc1 = _gc(point1, bearing1, _1_) # great circle p1 + b1

638 gc2 = _gc(point2, bearing2, _2_) # great circle p2 + b2

639

640 n = gc1.cross(gc2, raiser=_point_) # n-vector of intersection point

641 h = point1._havg(point2, h=height)

642 kwds = _xkwds(LatLon_and_kwds, height=h)

643 return n.toLatLon(**kwds) # Nvector(n.x, n.y, n.z).toLatLon(...)

644

645

646 def _trilaterate(point1, distance1, point2, distance2, point3, distance3,

647 radius=R_M, height=None, useZ=False,

648 wrap=False, **LatLon_and_kwds):

649 # (INTERNAL) Locate a point at given distances from

650 # three other points, see LatLon.triangulate above

651

652 def _nr2(p, d, r, _i_, *qs): # .toNvector and angular distance squared

653 for q in qs:

654 if _isequalTo(p, q, eps=EPS):

655 raise _ValueError(points=p, txt=_coincident_)

656 return p.toNvector(), (Scalar(d, name=_distance_ + _i_) / r)**2

657

658 p1, r = point1, Radius_(radius)

659 p2, p3, _ = _unrollon3(p1, point2, point3, wrap)

660

661 n1, r12 = _nr2(p1, distance1, r, _1_)

662 n2, r22 = _nr2(p2, distance2, r, _2_, p1)

663 n3, r32 = _nr2(p3, distance3, r, _3_, p1, p2)

664

665 # the following uses x,y coordinate system with origin at n1, x axis n1->n2

666 y = n3.minus(n1)

667 x = n2.minus(n1)

668 z = None

669

670 d = x.length # distance n1->n2

671 if d > EPS_2: # and y.length > EPS_2:

672 X = x.unit() # unit vector in x direction n1->n2

673 i = X.dot(y) # signed magnitude of x component of n1->n3

674 Y = y.minus(X.times(i)).unit() # unit vector in y direction

675 j = Y.dot(y) # signed magnitude of y component of n1->n3

676 if fabs(j) > EPS_2:

677 # courtesy of U{Carlos Freitas<https://GitHub.com/mrJean1/PyGeodesy/issues/33>}

678 x = fsumf_(r12, -r22, d**2) / (d * _2_0) # n1->intersection x- and ...

679 y = fsumf_(r12, -r32, i**2, j**2, x * i * _N_2_0) / (j * _2_0) # ... y-component

680 # courtesy of U{AleixDev<https://GitHub.com/mrJean1/PyGeodesy/issues/43>}

681 z = fsumf_(max(r12, r22, r32), -(x**2), -(y**2)) # XXX not just r12!

682 if z > EPS:

683 n = n1.plus(X.times(x)).plus(Y.times(y))

684 if useZ: # include Z component

685 Z = X.cross(Y) # unit vector perpendicular to plane

686 n = n.plus(Z.times(sqrt(z)))

687 if height is None:

688 h = fidw((point1.height, point2.height, point3.height),

689 map1(fabs, distance1, distance2, distance3))

690 else:

691 h = Height(height)

692 kwds = _xkwds(LatLon_and_kwds, height=h)

693 return n.toLatLon(**kwds) # Nvector(n.x, n.y, n.z).toLatLon(...)

694

695 # no intersection, d < EPS_2 or fabs(j) < EPS_2 or z < EPS

696 t = _SPACE_(_no_, _intersection_, NN)

697 raise IntersectionError(point1=point1, distance1=distance1,

698 point2=point2, distance2=distance2,

699 point3=point3, distance3=distance3,

700 txt=unstr(t, z=z, useZ=useZ, wrap=wrap))

701

702

703__all__ += _ALL_DOCS(LatLonNvectorBase, NvectorBase, sumOf) # classes

704

705# **) MIT License

706#

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

708#

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

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

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

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

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

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

715#

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

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

718#

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

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

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

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

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

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

725# OTHER DEALINGS IN THE SOFTWARE.