Coverage for pygeodesy/nvectorBase.py: 97%

233 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-21 13:14 -0400

1 

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

3 

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

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

6 

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

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

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

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

11''' 

12 

13from pygeodesy.basics import len2, map1 

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

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

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

17 _xkwds, _xkwds_pop 

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

19from pygeodesy.fsums import fsum, fsum_ 

20from pygeodesy.formy import n_xyz2latlon, n_xyz2philam, _spherical_datum 

21from pygeodesy.interns import MISSING, NN, _1_, _2_, _3_, _bearing_, \ 

22 _coincident_, _COMMASPACE_, _distance_, _h_, \ 

23 _intersection_, _no_, _NorthPole_, _points_, \ 

24 _pole_, _SPACE_, _SouthPole_, _UNDER 

25from pygeodesy.latlonBase import LatLonBase, _ALL_DOCS, _MODS 

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

27from pygeodesy.named import notImplemented, _xother3 

28from pygeodesy.namedTuples import Trilaterate5Tuple, Vector3Tuple, \ 

29 Vector4Tuple 

30from pygeodesy.props import deprecated_method, property_doc_, \ 

31 Property_RO, _update_all 

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

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

34# from pygeodesy.utily import sincos2d # from vector3d 

35from pygeodesy.vector3d import Vector3d, sumOf as _sumOf, sincos2d, _xyzhdn3 

36 

37from math import fabs, sqrt # atan2, cos, sin 

38 

39__all__ = (_NorthPole_, _SouthPole_) # constants 

40__version__ = '23.03.19' 

41 

42 

43class NvectorBase(Vector3d): # XXX kept private 

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

45 ''' 

46 _datum = None # L{Datum}, overriden 

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

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

49 

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

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

52 

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

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

55 L{Vector4Tuple}). 

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

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

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

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

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

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

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

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

64 

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

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

67 L{Vector3Tuple} or L{Vector4Tuple} or 

68 invalid B{C{datum}}. 

69 

70 @example: 

71 

72 >>> from pygeodesy.sphericalNvector import Nvector 

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

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

75 ''' 

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

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

78 if h: 

79 self.h = h 

80 if d is not None: 

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

82 

83 @Property_RO 

84 def datum(self): 

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

86 ''' 

87 return self._datum 

88 

89 @Property_RO 

90 def Ecef(self): 

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

92 ''' 

93 return _MODS.ecef.EcefKarney # default 

94 

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

96 def h(self): 

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

98 ''' 

99 return self._h 

100 

101 @h.setter # PYCHOK setter! 

102 def h(self, h): 

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

104 

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

106 

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

108 ''' 

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

110 if self._h != h: 

111 _update_all(self) 

112 self._h = h 

113 

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

115 def H(self): 

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

117 ''' 

118 return self._H 

119 

120 @H.setter # PYCHOK setter! 

121 def H(self, H): 

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

123 ''' 

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

125 

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

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

128 

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

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

131 

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

133 ''' 

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

135 

136 @Property_RO 

137 def isEllipsoidal(self): 

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

139 ''' 

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

141 

142 @Property_RO 

143 def isSpherical(self): 

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

145 ''' 

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

147 

148 @Property_RO 

149 def lam(self): 

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

151 ''' 

152 return self.philam.lam 

153 

154 @Property_RO 

155 def lat(self): 

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

157 ''' 

158 return self.latlon.lat 

159 

160 @Property_RO 

161 def latlon(self): 

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

163 ''' 

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

165 

166 @Property_RO 

167 def latlonheight(self): 

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

169 ''' 

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

171 

172 @Property_RO 

173 def latlonheightdatum(self): 

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

175 ''' 

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

177 

178 @Property_RO 

179 def lon(self): 

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

181 ''' 

182 return self.latlon.lon 

183 

184 @Property_RO 

185 def phi(self): 

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

187 ''' 

188 return self.philam.phi 

189 

190 @Property_RO 

191 def philam(self): 

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

193 ''' 

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

195 

196 @Property_RO 

197 def philamheight(self): 

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

199 ''' 

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

201 

202 @Property_RO 

203 def philamheightdatum(self): 

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

205 ''' 

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

207 

208 @deprecated_method 

209 def to2ab(self): # PYCHOK no cover 

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

211 

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

213 ''' 

214 return self.philam 

215 

216 @deprecated_method 

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

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

219 

220 @kwarg height: Optional height, overriding this 

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

222 

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

224 

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

226 ''' 

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

228 self.philam.to3Tuple(height) 

229 

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

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

232 

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

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

235 (C{Cartesian}). 

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

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

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

239 

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

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

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

243 available. 

244 

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

246 argument. 

247 

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

249 

250 @example: 

251 

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

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

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

255 ''' 

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

257 E = d.ellipsoid 

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

259 

260 x, y, z = self.x, self.y, self.z 

261 # Kenneth Gade eqn 22 

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

263 r = h + n * E.a2_b2 

264 

265 x *= r 

266 y *= r 

267 z *= n + h 

268 

269 if Cartesian is None: 

270 r = self.Ecef(d).reverse(x, y, z, M=True) 

271 else: 

272 kwds = _xkwds(Cartesian_kwds, datum=d) # h=0 

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

274 return self._xnamed(r) 

275 

276 @deprecated_method 

277 def to2ll(self): # PYCHOK no cover 

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

279 

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

281 ''' 

282 return self.latlon 

283 

284 @deprecated_method 

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

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

287 

288 @kwarg height: Optional height, overriding this 

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

290 

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

292 

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

294 ''' 

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

296 self.latlon.to3Tuple(height) 

297 

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

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

300 

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

302 height (C{meter}). 

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

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

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

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

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

308 

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

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

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

312 

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

314 argument. 

315 

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

317 

318 @example: 

319 

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

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

322 ''' 

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

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

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

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

327 if LatLon is None: 

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

329 else: 

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

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

332 return r 

333 

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

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

336 

337 Height component is only included if non-zero. 

338 

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

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

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

342 

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

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

345 

346 @example: 

347 

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

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

350 ''' 

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

352 if self.h: 

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

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

355 

356 def toVector3d(self, norm=True): 

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

358 the height}. 

359 

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

361 

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

363 ''' 

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

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

366 

367 @deprecated_method 

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

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

370 ''' 

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

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

373 

374 def unit(self, ll=None): 

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

376 

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

378 

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

380 ''' 

381 return _xattrs(Vector3d.unit(self, ll=ll), _UNDER(_h_)) 

382 

383 @Property_RO 

384 def xyzh(self): 

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

386 ''' 

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

388 

389 

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

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

392 

393 

394class _N_vector_(NvectorBase): 

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

396 ''' 

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

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

399 if h: 

400 self._h = h 

401 if name: 

402 self.name = name 

403 

404 

405class LatLonNvectorBase(LatLonBase): 

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

407 and spherical C{LatLon} classes. 

408 ''' 

409 

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

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

412 

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

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

415 ''' 

416 if updated: 

417 _Nv = _xkwds_pop(setters, _Nv=None) 

418 if _Nv is not None: 

419 if _Nv._fromll is not None: 

420 _Nv._fromll = None 

421 self._Nv = None 

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

423 

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

425# '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}. 

426# ''' 

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

428 

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

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

431 ''' 

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

433 

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

435 '''Refined class comparison. 

436 

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

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

439 keyword arguments. 

440 

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

442 

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

444 ''' 

445 if other: 

446 other0 = other[0] 

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

448 return other0 

449 

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

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

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

453 return other 

454 

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

456 '''Convert this point to C{Nvector} components, I{including 

457 height}. 

458 

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

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

461 

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

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

464 

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

466 argument. 

467 ''' 

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

469 

470 def triangulate(self, bearing1, other, bearing2, height=None): 

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

472 at this and the other point. 

473 

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

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

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

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

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

479 

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

481 

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

483 

484 @raise Valuerror: Points coincide. 

485 

486 @example: 

487 

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

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

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

491 ''' 

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

493 height=height, LatLon=self.classof) 

494 

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

496 radius=R_M, height=None, useZ=False): 

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

498 

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

500 as B{C{radius}}). 

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

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

503 B{C{radius}}). 

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

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

506 B{C{radius}}). 

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

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

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

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

511 

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

513 

514 @raise IntersectionError: No intersection, trilateration failed. 

515 

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

517 

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

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

520 

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

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

523 scripts/latlong-vectors.html>} and method C{LatLon.trilaterate2} 

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

525 ''' 

526 return _trilaterate(self, distance1, 

527 self.others(point2=point2), distance2, 

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

529 radius=radius, height=height, useZ=useZ, 

530 LatLon=self.classof) 

531 

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

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

534 '''B{Not implemented} for C{B{area}=True} or C{B{wrap}=True} 

535 and falls back to method C{trilaterate} otherwise. 

536 

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

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

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

540 margin and count C{n = 1}. 

541 

542 @raise IntersectionError: No intersection, trilateration failed. 

543 

544 @raise NotImplementedError: Keyword argument C{B{area}=True} or 

545 C{B{wrap}=True} not (yet) supported. 

546 

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

548 

549 @raise ValueError: Coincident B{C{points}} or invalid B{C{distance1}}, 

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

551 ''' 

552 if area or wrap: 

553 notImplemented(self, area=area, wrap=wrap) 

554 

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

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

557 radius=radius, height=None, useZ=True, 

558 LatLon=self.classof) 

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

560 # like function C{.latlonBase._trilaterate5} 

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

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

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

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

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

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

567 

568 

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

570 '''Return the vectorial sum of two or more n-vectors. 

571 

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

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

574 or C{None}. 

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

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

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

578 

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

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

581 

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

583 ''' 

584 n, nvectors = len2(nvectors) 

585 if n < 1: 

586 raise VectorError(nvectors=n, txt=MISSING) 

587 

588 if h is None: 

589 h = fsum(v.h for v in nvectors) / float(n) 

590 

591 if Vector is None: 

592 r = _sumOf(nvectors, Vector=Vector3Tuple).to4Tuple(h) 

593 else: 

594 r = _sumOf(nvectors, Vector=Vector, **_xkwds(Vector_kwds, h=h)) 

595 return r 

596 

597 

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

599 **LatLon_LatLon_kwds): 

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

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

602 

603 def _gc(p, b, _i_): 

604 n = p.toNvector() 

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

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

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

608 dest = de.times(s) 

609 dnct = dn.times(c) 

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

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

612 

613 if point1.isequalTo(point2, EPS): 

614 raise _ValueError(points=point2, txt=_coincident_) 

615 

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

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

618 

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

620 

621 h = point1._havg(point2) if height is None else Height(height) 

622 kwds = _xkwds(LatLon_LatLon_kwds, height=h) 

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

624 

625 

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

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

628 **LatLon_LatLon_kwds): 

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

630 # three other points, see LatLon.triangulate above 

631 

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

633 for q in qs: 

634 if p.isequalTo(q, EPS): 

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

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

637 

638 r = Radius_(radius) 

639 

640 n1, r12 = _nd2(point1, distance1, r, _1_) 

641 n2, r22 = _nd2(point2, distance2, r, _2_, point1) 

642 n3, r32 = _nd2(point3, distance3, r, _3_, point1, point2) 

643 

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

645 y = n3.minus(n1) 

646 x = n2.minus(n1) 

647 z = None 

648 

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

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

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

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

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

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

655 if fabs(j) > EPS_2: 

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

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

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

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

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

661 if z > EPS: 

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

663 if useZ: # include Z component 

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

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

666 if height is None: 

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

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

669 else: 

670 h = Height(height) 

671 kwds = _xkwds(LatLon_LatLon_kwds, height=h) 

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

673 

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

675 t = _SPACE_(_no_, _intersection_, NN) 

676 raise IntersectionError(point1=point1, distance1=distance1, 

677 point2=point2, distance2=distance2, 

678 point3=point3, distance3=distance3, 

679 txt=unstr(t, z=z, useZ=useZ)) 

680 

681 

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

683 

684# **) MIT License 

685# 

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

687# 

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

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

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

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

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

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

694# 

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

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

697# 

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

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

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

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

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

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

704# OTHER DEALINGS IN THE SOFTWARE.