Coverage for pygeodesy/nvectorBase.py: 95%

249 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-12-02 13:46 -0500

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 

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

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, fsumf_ 

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

21 _spherical_datum 

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

23 _COMMASPACE_, _distance_, _h_, _insufficient_, \ 

24 _intersection_, _no_, _NorthPole_, _point_, \ 

25 _pole_, _SPACE_, _SouthPole_, _under 

26from pygeodesy.latlonBase import LatLonBase, _ALL_DOCS, _ALL_LAZY, _MODS 

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

28from pygeodesy.named import notImplemented, _xother3 

29from pygeodesy.namedTuples import Trilaterate5Tuple, Vector3Tuple, \ 

30 Vector4Tuple, map1 

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

32 property_RO, _update_all 

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

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

35from pygeodesy.utily import sincos2d, _unrollon, _unrollon3 

36from pygeodesy.vector3d import Vector3d, _xyzhdn3 

37 

38from math import fabs, sqrt 

39 

40__all__ = _ALL_LAZY.nvectorBase 

41__version__ = '23.11.29' 

42 

43 

44class 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, once}. 

93 ''' 

94 NvectorBase.Ecef = E = _MODS.ecef.EcefKarney # overwrite property_RO 

95 return E 

96 

97 @property_RO 

98 def ellipsoidalNvector(self): 

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

100 ''' 

101 return False 

102 

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

104 def h(self): 

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

106 ''' 

107 return self._h 

108 

109 @h.setter # PYCHOK setter! 

110 def h(self, h): 

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

112 

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

114 

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

116 ''' 

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

118 if self._h != h: 

119 _update_all(self) 

120 self._h = h 

121 

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

123 def H(self): 

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

125 ''' 

126 return self._H 

127 

128 @H.setter # PYCHOK setter! 

129 def H(self, H): 

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

131 ''' 

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

133 

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

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

136 

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

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

139 

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

141 ''' 

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

143 

144 @Property_RO 

145 def isEllipsoidal(self): 

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

147 ''' 

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

149 

150 @Property_RO 

151 def isSpherical(self): 

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

153 ''' 

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

155 

156 @Property_RO 

157 def lam(self): 

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

159 ''' 

160 return self.philam.lam 

161 

162 @Property_RO 

163 def lat(self): 

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

165 ''' 

166 return self.latlon.lat 

167 

168 @Property_RO 

169 def latlon(self): 

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

171 ''' 

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

173 

174 @Property_RO 

175 def latlonheight(self): 

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

177 ''' 

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

179 

180 @Property_RO 

181 def latlonheightdatum(self): 

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

183 ''' 

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

185 

186 @Property_RO 

187 def lon(self): 

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

189 ''' 

190 return self.latlon.lon 

191 

192 @Property_RO 

193 def phi(self): 

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

195 ''' 

196 return self.philam.phi 

197 

198 @Property_RO 

199 def philam(self): 

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

201 ''' 

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

203 

204 @Property_RO 

205 def philamheight(self): 

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

207 ''' 

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

209 

210 @Property_RO 

211 def philamheightdatum(self): 

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

213 ''' 

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

215 

216 @property_RO 

217 def sphericalNvector(self): 

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

219 ''' 

220 return False 

221 

222 @deprecated_method 

223 def to2ab(self): # PYCHOK no cover 

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

225 

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

227 ''' 

228 return self.philam 

229 

230 @deprecated_method 

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

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

233 

234 @kwarg height: Optional height, overriding this 

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

236 

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

238 

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

240 ''' 

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

242 self.philam.to3Tuple(height) 

243 

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

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

246 

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

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

249 (C{Cartesian}). 

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

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

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

253 

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

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

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

257 available. 

258 

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

260 argument. 

261 

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

263 

264 @example: 

265 

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

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

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

269 ''' 

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

271 E = D.ellipsoid 

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

273 

274 x, y, z = self.x, self.y, self.z 

275 # Kenneth Gade eqn 22 

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

277 r = h + n * E.a2_b2 

278 

279 x *= r 

280 y *= r 

281 z *= h + n 

282 

283 if Cartesian is None: 

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

285 else: 

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

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

288 return self._xnamed(r) 

289 

290 @deprecated_method 

291 def to2ll(self): # PYCHOK no cover 

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

293 

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

295 ''' 

296 return self.latlon 

297 

298 @deprecated_method 

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

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

301 

302 @kwarg height: Optional height, overriding this 

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

304 

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

306 

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

308 ''' 

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

310 self.latlon.to3Tuple(height) 

311 

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

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

314 

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

316 height (C{meter}). 

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

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

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

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

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

322 

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

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

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

326 

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

328 argument. 

329 

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

331 

332 @example: 

333 

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

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

336 ''' 

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

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

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

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

341 if LatLon is None: 

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

343 else: 

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

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

346 return r 

347 

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

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

350 

351 Height component is only included if non-zero. 

352 

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

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

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

356 

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

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

359 

360 @example: 

361 

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

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

364 ''' 

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

366 if self.h: 

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

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

369 

370 def toVector3d(self, norm=True): 

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

372 the height}. 

373 

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

375 

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

377 ''' 

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

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

380 

381 @deprecated_method 

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

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

384 ''' 

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

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

387 

388 def unit(self, ll=None): 

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

390 

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

392 

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

394 ''' 

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

396 

397 @Property_RO 

398 def xyzh(self): 

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

400 ''' 

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

402 

403 

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

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

406 

407 

408class _N_vector_(NvectorBase): 

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

410 ''' 

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

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

413 if h: 

414 self._h = h 

415 if name: 

416 self.name = name 

417 

418 

419class LatLonNvectorBase(LatLonBase): 

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

421 and spherical C{LatLon} classes. 

422 ''' 

423 

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

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

426 

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

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

429 ''' 

430 if updated: 

431 _Nv = _xkwds_pop(setters, _Nv=None) 

432 if _Nv is not None: 

433 if _Nv._fromll is not None: 

434 _Nv._fromll = None 

435 self._Nv = None 

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

437 

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

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

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

441 

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

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

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

445 

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

447 '''Refined class comparison. 

448 

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

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

451 keyword arguments. 

452 

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

454 

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

456 ''' 

457 if other: 

458 other0 = other[0] 

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

460 return other0 

461 

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

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

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

465 return other 

466 

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

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

469 

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

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

472 

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

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

475 

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

477 argument. 

478 ''' 

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

480 

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

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

483 at this and the other point. 

484 

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

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

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

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

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

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

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

492 

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

494 

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

496 

497 @raise Valuerror: Points coincide. 

498 

499 @example: 

500 

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

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

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

504 ''' 

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

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

507 

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

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

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

511 

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

513 as B{C{radius}}). 

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

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

516 B{C{radius}}). 

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

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

519 B{C{radius}}). 

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

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

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

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

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

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

526 

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

528 

529 @raise IntersectionError: No intersection, trilateration failed. 

530 

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

532 

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

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

535 

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

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

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

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

540 ''' 

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

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

543 radius=radius, height=height, useZ=useZ, 

544 wrap=wrap, LatLon=self.classof) 

545 

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

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

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

549 C{trilaterate} otherwise. 

550 

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

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

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

554 margin and count C{n = 1}. 

555 

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

557 (yet) supported. 

558 

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

560 ''' 

561 if area: 

562 notImplemented(self, area=area) 

563 

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

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

566 radius=radius, useZ=True, wrap=wrap, 

567 LatLon=self.classof) 

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

569 # like function C{.latlonBase._trilaterate5} 

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

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

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

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

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

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

576 

577 

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

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

580 ''' 

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

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

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

584 X += v.x 

585 Y += v.y 

586 Z += v.z 

587 H += v.h 

588 if n < 1: 

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

590 

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

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

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

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

595 

596 

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

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

599 

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

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

602 or C{None}. 

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

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

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

606 

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

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

609 

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

611 ''' 

612 try: 

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

614 except (TypeError, ValueError) as x: 

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

616 

617 

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

619 wrap=False, **LatLon_and_kwds): 

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

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

622 

623 def _gc(p, b, _i_): 

624 n = p.toNvector() 

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

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

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

628 dest = de.times(s) 

629 dnct = dn.times(c) 

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

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

632 

633 if wrap: 

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

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

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

637 

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

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

640 

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

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

643 kwds = _xkwds(LatLon_and_kwds, height=h) 

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

645 

646 

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

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

649 wrap=False, **LatLon_and_kwds): 

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

651 # three other points, see LatLon.triangulate above 

652 

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

654 for q in qs: 

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

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

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

658 

659 p1, r = point1, Radius_(radius) 

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

661 

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

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

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

665 

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

667 y = n3.minus(n1) 

668 x = n2.minus(n1) 

669 z = None 

670 

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

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

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

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

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

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

677 if fabs(j) > EPS_2: 

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

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

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

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

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

683 if z > EPS: 

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

685 if useZ: # include Z component 

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

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

688 if height is None: 

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

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

691 else: 

692 h = Height(height) 

693 kwds = _xkwds(LatLon_and_kwds, height=h) 

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

695 

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

697 t = _SPACE_(_no_, _intersection_, NN) 

698 raise IntersectionError(point1=point1, distance1=distance1, 

699 point2=point2, distance2=distance2, 

700 point3=point3, distance3=distance3, 

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

702 

703 

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

705 

706# **) MIT License 

707# 

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

709# 

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

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

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

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

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

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

716# 

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

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

719# 

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

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

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

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

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

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

726# OTHER DEALINGS IN THE SOFTWARE.