1

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

3

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

5

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

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

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

9 L{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

14 from pygeodesy.ellipsoidalBase import CartesianEllipsoidalBase, \

15 _nearestOn, _WGS84

16 from pygeodesy.ellipsoidalBaseDI import LatLonEllipsoidalBaseDI, _TOL_M, \

17 _intersection3, _intersections2

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

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

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

21 from 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

30 class 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

52 class 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

93 def 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

118 def 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

169 def 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

217 def 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

243 def 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

285 def 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.