Coverage for pygeodesy/ellipsoidalGeodSolve.py: 100%

40 statements  

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

1 

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

3 

4u'''Exact ellipsoidal geodesy, intended I{for testing purposes only}. 

5 

6Ellipsoidal geodetic (lat-/longitude) L{LatLon} and geocentric 

7(ECEF) L{Cartesian} classes and functions L{areaOf}, L{intersections2}, 

8L{isclockwise}, L{nearestOn} and L{perimeterOf} based on module 

9L{geodsolve}, a wrapper invoking I{Karney}'s U{GeodSolve 

10<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>} utility. 

11''' 

12 

13# from pygeodesy.datums import _WGS84 # from .ellipsoidalBase 

14from pygeodesy.ellipsoidalBase import CartesianEllipsoidalBase, \ 

15 _nearestOn, _WGS84 

16from pygeodesy.ellipsoidalBaseDI import LatLonEllipsoidalBaseDI, _TOL_M, \ 

17 _intersection3, _intersections2 

18# from pygeodesy.errors import _xkwds # from .karney 

19from pygeodesy.karney import fabs, _polygon, Property_RO, _xkwds 

20from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _ALL_OTHER 

21from pygeodesy.points import _areaError, ispolar # PYCHOK exported 

22# from pygeodesy.props import Property_RO # from .karney 

23 

24# from math import fabs # from .karney 

25 

26__all__ = _ALL_LAZY.ellipsoidalGeodSolve 

27__version__ = '23.04.11' 

28 

29 

30class Cartesian(CartesianEllipsoidalBase): 

31 '''Extended to convert exact L{Cartesian} to exact L{LatLon} points. 

32 ''' 

33 

34 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon, datum=None 

35 '''Convert this cartesian point to an exact geodetic point. 

36 

37 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword 

38 arguments as C{datum}. Use C{B{LatLon}=..., 

39 B{datum}=...} to override this L{LatLon} 

40 class or specify C{B{LatLon}=None}. 

41 

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

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

44 with C{C} and C{M} if available. 

45 

46 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument. 

47 ''' 

48 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum) 

49 return CartesianEllipsoidalBase.toLatLon(self, **kwds) 

50 

51 

52class LatLon(LatLonEllipsoidalBaseDI): 

53 '''An ellipsoidal L{LatLon} like L{ellipsoidalKarney.LatLon} but using (exact) 

54 geodesic I{wrapper} L{GeodesicSolve} to compute the geodesic distance, 

55 initial and final bearing (azimuths) between two given points or the 

56 destination point given a start point and an (initial) bearing. 

57 ''' 

58 

59 @Property_RO 

60 def Equidistant(self): 

61 '''Get the prefered azimuthal equidistant projection I{class} (L{EquidistantGeodSolve}). 

62 ''' 

63 return _MODS.azimuthal.EquidistantGeodSolve 

64 

65 @Property_RO 

66 def geodesicx(self): 

67 '''Get this C{LatLon}'s (exact) geodesic (L{GeodesicSolve}). 

68 ''' 

69 return self.datum.ellipsoid.geodsolve 

70 

71 geodesic = geodesicx # for C{._Direct} and C{._Inverse} 

72 

73 def toCartesian(self, **Cartesian_datum_kwds): # PYCHOK Cartesian=Cartesian, datum=None 

74 '''Convert this point to exact cartesian (ECEF) coordinates. 

75 

76 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}} 

77 and other keyword arguments, ignored if C{B{Cartesian} is None}. 

78 Use C{B{Cartesian}=...} to override this L{Cartesian} class 

79 or set C{B{Cartesian} is None}. 

80 

81 @return: The cartesian (ECEF) coordinates (L{Cartesian}) or if 

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

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

84 available. 

85 

86 @raise TypeError: Invalid B{C{Cartesian}}, B{C{datum}} or other 

87 B{C{Cartesian_datum_kwds}}. 

88 ''' 

89 kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian, datum=self.datum) 

90 return LatLonEllipsoidalBaseDI.toCartesian(self, **kwds) 

91 

92 

93def areaOf(points, datum=_WGS84, wrap=True): 

94 '''Compute the area of an (ellipsoidal) polygon or composite. 

95 

96 @arg points: The polygon points (L{LatLon}[], L{BooleanFHP} or 

97 L{BooleanGH}). 

98 @kwarg datum: Optional datum (L{Datum}). 

99 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

100 

101 @return: Area (C{meter}, same as units of the B{C{datum}}'s 

102 ellipsoid axes, I{squared}). 

103 

104 @raise PointsError: Insufficient number of B{C{points}}. 

105 

106 @raise TypeError: Some B{C{points}} are not L{LatLon}. 

107 

108 @raise ValueError: Invalid C{B{wrap}=False}, unwrapped, unrolled 

109 longitudes not supported. 

110 

111 @see: Functions L{pygeodesy.areaOf}, L{ellipsoidalExact.areaOf}, 

112 L{ellipsoidalKarney.areaOf}, L{sphericalNvector.areaOf} 

113 and L{sphericalTrigonometry.areaOf}. 

114 ''' 

115 return fabs(_polygon(datum.ellipsoid.geodsolve, points, True, False, wrap)) 

116 

117 

118def intersection3(start1, end1, start2, end2, height=None, wrap=True, 

119 equidistant=None, tol=_TOL_M, LatLon=LatLon, **LatLon_kwds): 

120 '''Iteratively compute the intersection point of two lines, each defined 

121 by two (ellipsoidal) points or by an (ellipsoidal) start point and an 

122 (initial) bearing from North. 

123 

124 @arg start1: Start point of the first line (L{LatLon}). 

125 @arg end1: End point of the first line (L{LatLon}) or the initial bearing 

126 at the first point (compass C{degrees360}). 

127 @arg start2: Start point of the second line (L{LatLon}). 

128 @arg end2: End point of the second line (L{LatLon}) or the initial bearing 

129 at the second point (compass C{degrees360}). 

130 @kwarg height: Optional height at the intersection (C{meter}, conventionally) 

131 or C{None} for the mean height. 

132 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

133 @kwarg equidistant: An azimuthal equidistant projection (I{class} or function 

134 L{pygeodesy.equidistant}) or C{None} for the preferred 

135 C{B{start1}.Equidistant}. 

136 @kwarg tol: Tolerance for convergence and for skew line distance and length 

137 (C{meter}, conventionally). 

138 @kwarg LatLon: Optional class to return the intersection points (L{LatLon}) 

139 or C{None}. 

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

141 ignored if C{B{LatLon} is None}. 

142 

143 @return: An L{Intersection3Tuple}C{(point, outside1, outside2)} with C{point} 

144 a B{C{LatLon}} or if C{B{LatLon} is None}, a L{LatLon4Tuple}C{(lat, 

145 lon, height, datum)}. 

146 

147 @raise IntersectionError: Skew, colinear, parallel or otherwise 

148 non-intersecting lines or no convergence 

149 for the given B{C{tol}}. 

150 

151 @raise TypeError: Invalid or non-ellipsoidal B{C{start1}}, B{C{end1}}, 

152 B{C{start2}} or B{C{end2}} or invalid B{C{equidistant}}. 

153 

154 @note: For each line specified with an initial bearing, a pseudo-end point 

155 is computed as the C{destination} along that bearing at about 1.5 

156 times the distance from the start point to an initial gu-/estimate 

157 of the intersection point (and between 1/8 and 3/8 of the authalic 

158 earth perimeter). 

159 

160 @see: U{The B{ellipsoidal} case<https://GIS.StackExchange.com/questions/48937/ 

161 calculating-intersection-of-two-circles>} and U{Karney's paper 

162 <https://ArXiv.org/pdf/1102.1215.pdf>}, pp 20-21, section B{14. MARITIME 

163 BOUNDARIES} for more details about the iteration algorithm. 

164 ''' 

165 return _intersection3(start1, end1, start2, end2, height=height, wrap=wrap, 

166 equidistant=equidistant, tol=tol, LatLon=LatLon, **LatLon_kwds) 

167 

168 

169def intersections2(center1, radius1, center2, radius2, height=None, wrap=True, 

170 equidistant=None, tol=_TOL_M, LatLon=LatLon, **LatLon_kwds): 

171 '''Iteratively compute the intersection points of two circles, each defined 

172 by an (ellipsoidal) center point and a radius. 

173 

174 @arg center1: Center of the first circle (L{LatLon}). 

175 @arg radius1: Radius of the first circle (C{meter}, conventionally). 

176 @arg center2: Center of the second circle (L{LatLon}). 

177 @arg radius2: Radius of the second circle (C{meter}, same units as 

178 B{C{radius1}}). 

179 @kwarg height: Optional height for the intersection points (C{meter}, 

180 conventionally) or C{None} for the I{"radical height"} 

181 at the I{radical line} between both centers. 

182 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

183 @kwarg equidistant: An azimuthal equidistant projection (I{class} or 

184 function L{pygeodesy.equidistant}) or C{None} for 

185 the preferred C{B{center1}.Equidistant}. 

186 @kwarg tol: Convergence tolerance (C{meter}, same units as B{C{radius1}} 

187 and B{C{radius2}}). 

188 @kwarg LatLon: Optional class to return the intersection points (L{LatLon}) 

189 or C{None}. 

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

191 ignored if C{B{LatLon} is None}. 

192 

193 @return: 2-Tuple of the intersection points, each a B{C{LatLon}} instance 

194 or L{LatLon4Tuple}C{(lat, lon, height, datum)} if C{B{LatLon} is 

195 None}. For abutting circles, both points are the same instance, 

196 aka the I{radical center}. 

197 

198 @raise IntersectionError: Concentric, antipodal, invalid or non-intersecting 

199 circles or no convergence for the B{C{tol}}. 

200 

201 @raise TypeError: Invalid or non-ellipsoidal B{C{center1}} or B{C{center2}} 

202 or invalid B{C{equidistant}}. 

203 

204 @raise UnitError: Invalid B{C{radius1}}, B{C{radius2}} or B{C{height}}. 

205 

206 @see: U{The B{ellipsoidal} case<https://GIS.StackExchange.com/questions/48937/ 

207 calculating-intersection-of-two-circles>}, U{Karney's paper 

208 <https://ArXiv.org/pdf/1102.1215.pdf>}, pp 20-21, section B{14. MARITIME BOUNDARIES}, 

209 U{Circle-Circle<https://MathWorld.Wolfram.com/Circle-CircleIntersection.html>} and 

210 U{Sphere-Sphere<https://MathWorld.Wolfram.com/Sphere-SphereIntersection.html>} 

211 intersections. 

212 ''' 

213 return _intersections2(center1, radius1, center2, radius2, height=height, wrap=wrap, 

214 equidistant=equidistant, tol=tol, LatLon=LatLon, **LatLon_kwds) 

215 

216 

217def isclockwise(points, datum=_WGS84, wrap=True): 

218 '''Determine the direction of a path or polygon. 

219 

220 @arg points: The path or polygon points (C{LatLon}[]). 

221 @kwarg datum: Optional datum (L{Datum}). 

222 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

223 

224 @return: C{True} if B{C{points}} are clockwise, C{False} otherwise. 

225 

226 @raise PointsError: Insufficient number of B{C{points}}. 

227 

228 @raise TypeError: Some B{C{points}} are not C{LatLon}. 

229 

230 @raise ValueError: The B{C{points}} enclose a pole or zero 

231 area. 

232 

233 @see: L{pygeodesy.isclockwise}. 

234 ''' 

235 a = _polygon(datum.ellipsoid.geodsolve, points, True, False, wrap) 

236 if a < 0: 

237 return True 

238 elif a > 0: 

239 return False 

240 raise _areaError(points) 

241 

242 

243def nearestOn(point, point1, point2, within=True, height=None, wrap=False, 

244 equidistant=None, tol=_TOL_M, LatLon=LatLon, **LatLon_kwds): 

245 '''Iteratively locate the closest point on the geodesic between 

246 two other (ellipsoidal) points. 

247 

248 @arg point: Reference point (C{LatLon}). 

249 @arg point1: Start point of the geodesic (C{LatLon}). 

250 @arg point2: End point of the geodesic (C{LatLon}). 

251 @kwarg within: If C{True} return the closest point I{between} 

252 B{C{point1}} and B{C{point2}}, otherwise the 

253 closest point elsewhere on the geodesic (C{bool}). 

254 @kwarg height: Optional height for the closest point (C{meter}, 

255 conventionally) or C{None} or C{False} for the 

256 interpolated height. If C{False}, the closest 

257 takes the heights of the points into account. 

258 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

259 @kwarg equidistant: An azimuthal equidistant projection (I{class} 

260 or function L{pygeodesy.equidistant}) or C{None} 

261 for the preferred C{B{point}.Equidistant}. 

262 @kwarg tol: Convergence tolerance (C{meter}). 

263 @kwarg LatLon: Optional class to return the closest point 

264 (L{LatLon}) or C{None}. 

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

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

267 

268 @return: Closest point, a B{C{LatLon}} instance or if C{B{LatLon} 

269 is None}, a L{LatLon4Tuple}C{(lat, lon, height, datum)}. 

270 

271 @raise TypeError: Invalid or non-ellipsoidal B{C{point}}, B{C{point1}} 

272 or B{C{point2}} or invalid B{C{equidistant}}. 

273 

274 @raise ValueError: No convergence for the B{C{tol}}. 

275 

276 @see: U{The B{ellipsoidal} case<https://GIS.StackExchange.com/questions/48937/ 

277 calculating-intersection-of-two-circles>} and U{Karney's paper 

278 <https://ArXiv.org/pdf/1102.1215.pdf>}, pp 20-21, section B{14. MARITIME 

279 BOUNDARIES} for more details about the iteration algorithm. 

280 ''' 

281 return _nearestOn(point, point1, point2, within=within, height=height, wrap=wrap, 

282 equidistant=equidistant, tol=tol, LatLon=LatLon, **LatLon_kwds) 

283 

284 

285def perimeterOf(points, closed=False, datum=_WGS84, wrap=True): 

286 '''Compute the perimeter of an (ellipsoidal) polygon or composite. 

287 

288 @arg points: The polygon points (L{LatLon}[], L{BooleanFHP} or 

289 L{BooleanGH}). 

290 @kwarg closed: Optionally, close the polygon (C{bool}). 

291 @kwarg datum: Optional datum (L{Datum}). 

292 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

293 

294 @return: Perimeter (C{meter}, same as units of the B{C{datum}}'s 

295 ellipsoid axes). 

296 

297 @raise PointsError: Insufficient number of B{C{points}}. 

298 

299 @raise TypeError: Some B{C{points}} are not L{LatLon}. 

300 

301 @raise ValueError: Invalid C{B{wrap}=False}, unwrapped, unrolled 

302 longitudes not supported or C{B{closed}=False} 

303 with C{B{points}} a composite. 

304 

305 @see: Functions L{pygeodesy.perimeterOf}, L{ellipsoidalExact.perimeterOf}, 

306 L{ellipsoidalKarney.perimeterOf}, L{sphericalNvector.perimeterOf} 

307 and L{sphericalTrigonometry.perimeterOf}. 

308 ''' 

309 return _polygon(datum.ellipsoid.geodsolve, points, closed, True, wrap) 

310 

311 

312__all__ += _ALL_OTHER(Cartesian, LatLon, # classes 

313 areaOf, # functions 

314 intersection3, intersections2, isclockwise, ispolar, 

315 nearestOn, perimeterOf) 

316 

317# **) MIT License 

318# 

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

320# 

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

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

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

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

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

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

327# 

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

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

330# 

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

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

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

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

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

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

337# OTHER DEALINGS IN THE SOFTWARE.