Coverage for pyguymer3/geo/great_circle.py: 63%
86 statements
« prev ^ index » next coverage.py v7.9.2, created at 2025-07-08 18:47 +0000
1#!/usr/bin/env python3
3# Define function ...
4def great_circle(
5 lon1,
6 lat1,
7 lon2,
8 lat2,
9 /,
10 *,
11 debug = __debug__,
12 eps = 1.0e-12,
13 maxdist = None,
14 nIter = 100,
15 npoint = None,
16 prefix = ".",
17 ramLimit = 1073741824,
18):
19 """Calculate the great circle that connects two coordinates.
21 This function reads in two coordinates (in degrees) on the surface of the
22 Earth and calculates the great circle that connects them, correctly handling
23 crossing over the anti-meridian (should it occur).
25 Parameters
26 ----------
27 lon1 : float
28 the longitude of the first coordinate (in degrees)
29 lat1 : float
30 the latitude of the first coordinate (in degrees)
31 lon2 : float
32 the longitude of the second coordinate (in degrees)
33 lat2 : float
34 the latitude of the second coordinate (in degrees)
35 debug : bool, optional
36 print debug messages (in metres)
37 eps : float, optional
38 the tolerance of the Vincenty formula iterations
39 maxdist : float, optional
40 the maximum distance between points along the great circle
41 nIter : int, optional
42 the maximum number of iterations (particularly the Vincenty formula)
43 npoint : int, optional
44 the number of points along the great circle
45 prefix : str, optional
46 change the name of the output debugging CSVs
47 ramLimit : int, optional
48 the maximum RAM usage of each "large" array (in bytes)
50 Returns
51 -------
52 line : shapely.geometry.linestring.LineString, shapely.geometry.multilinestring.MultiLineString
53 the great circle
55 Notes
56 -----
57 If neither maxdist nor npoint are provided then npoint will default to 5.
58 For large great circles, prettier results are obtained by using maxdist
59 rather than npoint.
61 Copyright 2017 Thomas Guymer [1]_
63 References
64 ----------
65 .. [1] PyGuymer3, https://github.com/Guymer/PyGuymer3
66 """
68 # Import special modules ...
69 try:
70 import numpy
71 except:
72 raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None
73 try:
74 import shapely
75 import shapely.geometry
76 except:
77 raise Exception("\"shapely\" is not installed; run \"pip install --user Shapely\"") from None
79 # Import sub-functions ...
80 from .calc_dist_between_two_locs import calc_dist_between_two_locs
81 from .check import check
82 from .find_middle_of_great_circle import find_middle_of_great_circle
83 from .find_point_on_great_circle import find_point_on_great_circle
84 from ..interpolate import interpolate
85 from ..consts import CIRCUMFERENCE_OF_EARTH
87 # **************************************************************************
89 # Set default value ...
90 if maxdist is None and npoint is None:
91 npoint = 5 # [#]
93 # Check inputs ...
94 if isinstance(maxdist, float):
95 assert maxdist >= 10.0, f"the maximum distance is too small ({maxdist:,.1f}m < {10.0:,.1f}m)"
96 assert maxdist <= 0.5 * CIRCUMFERENCE_OF_EARTH, f"the maximum distance is too large ({maxdist:,.1f}m > {0.5 * CIRCUMFERENCE_OF_EARTH:,.1f}m)"
97 elif isinstance(npoint, int):
98 assert npoint >= 3, f"the number of points is too small ({npoint:,d} < {3:,d})"
99 else:
100 raise Exception("\"maxdist\" is not a float and \"npoint\" is not an integer") from None
102 # **************************************************************************
104 # Skip if the start- and end-points are the same ...
105 if lon1 == lon2 and lat1 == lat2:
106 return None
108 # Check if the user wants to make the great circle using either "maxdist" or
109 # "npoints" ...
110 if isinstance(maxdist, float):
111 # Create a poorly resolved great circle ...
112 points = [
113 (lon1, lat1),
114 find_point_on_great_circle(0.2, lon1, lat1, lon2, lat2),
115 find_point_on_great_circle(0.4, lon1, lat1, lon2, lat2),
116 find_point_on_great_circle(0.6, lon1, lat1, lon2, lat2),
117 (lon2, lat2),
118 ] # [°], [°]
120 # Initialize flag ...
121 addedPoint = True
123 # Commence infinite loop ...
124 while addedPoint:
125 # Reset the flag ..
126 addedPoint = False
128 # Loop over points ...
129 for ipoint in range(1, len(points)):
130 # Calculate the distance between this point and the previous one ...
131 dist, _, _ = calc_dist_between_two_locs(
132 points[ipoint - 1][0],
133 points[ipoint - 1][1],
134 points[ipoint][0],
135 points[ipoint][1],
136 eps = eps,
137 nIter = nIter,
138 ) # [m]
140 # Check if the distance is too large ...
141 if dist > maxdist:
142 # Replace the list of points with one which has a new
143 # intermediate point inserted ...
144 points = points[:ipoint] + [
145 find_middle_of_great_circle(
146 points[ipoint - 1][0],
147 points[ipoint - 1][1],
148 points[ipoint][0],
149 points[ipoint][1],
150 )
151 ] + points[ipoint:] # [°], [°]
153 # Set the flag ...
154 addedPoint = True
156 # Stop looping over points and start the survey again ...
157 break
159 # Convert list to array and clean up ...
160 circle = numpy.array(points, dtype = numpy.float64) # [°]
161 del points
162 elif isinstance(npoint, int):
163 # Check array size ...
164 if npoint * 2 * 8 > ramLimit:
165 raise Exception(f"\"circle\" is going to be {npoint * 2 * 8:,d} bytes, which is larger than {ramLimit:,d} bytes") from None
167 # Initialize array ...
168 circle = numpy.zeros((npoint, 2), dtype = numpy.float64) # [°]
170 # Loop over points ...
171 for ipoint in range(npoint):
172 # Set point ...
173 circle[ipoint, :] = find_point_on_great_circle(float(ipoint) / float(npoint - 1), lon1, lat1, lon2, lat2) # [°]
174 else:
175 raise Exception("\"maxdist\" is not a float and \"npoint\" is not an integer") from None
177 # **************************************************************************
179 # Check if the great circle crosses the anti-meridian W to E ...
180 if lon2 > lon1 and (circle[1, 0] < lon1 or circle[-2, 0] > lon2):
181 if debug:
182 print("INFO: The great circle crosses the anti-meridian (point #2 is E of point #1 but the great circle goes W).")
184 # Find where it crosses the anti-meridian ...
185 i = numpy.diff(circle[:, 0]).argmax()
187 # Calculate the first intersection and convert to a LineString ...
188 x = -180.0 # [°]
189 y = interpolate(circle[i, 0], circle[i + 1, 0] - 360.0, circle[i, 1], circle[i + 1, 1], x) # [°]
190 line1 = shapely.geometry.linestring.LineString(numpy.append(circle[:i + 1, :], [[x, y]], axis = 0))
191 if debug:
192 check(line1, prefix = prefix)
194 # Calculate the second intersection and convert to a LineString ...
195 x = +180.0 # [°]
196 y = interpolate(circle[i, 0] + 360.0, circle[i + 1, 0], circle[i, 1], circle[i + 1, 1], x) # [°]
197 line2 = shapely.geometry.linestring.LineString(numpy.append([[x, y]], circle[i + 1:, :], axis = 0))
198 if debug:
199 check(line2, prefix = prefix)
201 # Convert to a MultiLineString ...
202 multiline = shapely.geometry.multilinestring.MultiLineString([line1, line2])
203 if debug:
204 check(multiline, prefix = prefix)
206 # Return answer ...
207 return multiline
209 # Check if the great circle crosses the anti-meridian E to W ...
210 if lon2 < lon1 and (circle[1, 0] > lon1 or circle[-2, 0] < lon2):
211 if debug:
212 print("INFO: The great circle crosses the anti-meridian (point #2 is W of point #1 but the great circle goes E).")
214 # Find where it crosses the anti-meridian ...
215 i = numpy.diff(circle[:, 0]).argmin()
217 # Calculate the first intersection and convert to a LineString ...
218 x = +180.0 # [°]
219 y = interpolate(circle[i, 0], circle[i + 1, 0] + 360.0, circle[i, 1], circle[i + 1, 1], x) # [°]
220 line1 = shapely.geometry.linestring.LineString(numpy.append(circle[:i + 1, :], [[x, y]], axis = 0))
221 if debug:
222 check(line1, prefix = prefix)
224 # Calculate the second intersection and convert to a LineString ...
225 x = -180.0 # [°]
226 y = interpolate(circle[i, 0] - 360.0, circle[i + 1, 0], circle[i, 1], circle[i + 1, 1], x) # [°]
227 line2 = shapely.geometry.linestring.LineString(numpy.append([[x, y]], circle[i + 1:, :], axis = 0))
228 if debug:
229 check(line2, prefix = prefix)
231 # Convert to a MultiLineString ...
232 multiline = shapely.geometry.multilinestring.MultiLineString([line1, line2])
233 if debug:
234 check(multiline, prefix = prefix)
236 # Return answer ...
237 return multiline
239 # **************************************************************************
241 # Convert to a LineString ...
242 line = shapely.geometry.linestring.LineString(circle)
243 if debug:
244 check(line, prefix = prefix)
246 # Return answer ...
247 return line