Coverage for pyguymer3/geo/find_middle_of_locsSrc/find_middle_of_locs_geodesicCircle.py: 63%
59 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 find_middle_of_locs_geodesicCircle(
5 lons,
6 lats,
7 /,
8 *,
9 angConv = 0.1,
10 conv = 1.0e3,
11 debug = __debug__,
12 eps = 1.0e-12,
13 iRefine = 0,
14 midLat = None,
15 midLon = None,
16 nAng = 9,
17 nIter = 100,
18 nRefine = 1,
19 pad = 10.0e3,
20 useSciPy = False,
21):
22 """Find the middle of some locations such that they are encompassed by the
23 smallest Geodesic circle possible.
24 """
26 # Import special modules ...
27 try:
28 import numpy
29 except:
30 raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None
31 try:
32 import scipy
33 except:
34 raise Exception("\"scipy\" is not installed; run \"pip install --user scipy\"") from None
36 # Import sub-functions ...
37 from .find_middle_of_locs_euclideanBox import find_middle_of_locs_euclideanBox
38 from ..calc_loc_from_loc_and_bearing_and_dist import calc_loc_from_loc_and_bearing_and_dist
39 from ..find_min_max_dist_bearing import find_min_max_dist_bearing
40 from ..max_dist import max_dist
41 from ...consts import RESOLUTION_OF_EARTH
43 # **************************************************************************
45 # Check arguments ...
46 assert isinstance(lons, numpy.ndarray), "\"lons\" is not a NumPy array"
47 assert isinstance(lats, numpy.ndarray), "\"lats\" is not a NumPy array"
49 # **************************************************************************
51 # Calculate the middle of the Euclidean bounding box to use as a decent
52 # starting guess (if the user has not provided one) ...
53 if midLon is None or midLat is None:
54 midLon, midLat, _ = find_middle_of_locs_euclideanBox(
55 lons,
56 lats,
57 debug = debug,
58 pad = -1.0,
59 ) # [°], [°]
61 # Find the maximum Geodesic distance from the middle to any location ...
62 maxDist = max_dist(
63 lons,
64 lats,
65 midLon,
66 midLat,
67 eps = eps,
68 nIter = nIter,
69 space = "GeodesicSpace",
70 ) # [m]
72 if debug:
73 print(f"INFO: The initial middle is ({midLon:.6f}°, {midLat:.6f}°) and the initial maximum Geodesic distance is {0.001 * maxDist:,.1f} km.")
75 # **************************************************************************
77 # Check if the input is already converged or if the user wants to use SciPy ...
78 if maxDist < conv:
79 pass
80 elif useSciPy:
81 # Use SciPy to find the minimum maximum Geodesic distance ...
82 ans = scipy.optimize.minimize(
83 lambda x: max_dist(
84 lons,
85 lats,
86 x[0],
87 x[1],
88 eps = eps,
89 nIter = nIter,
90 space = "GeodesicSpace",
91 ),
92 [midLon, midLat],
93 bounds = [
94 (-180.0, +180.0),
95 ( -90.0, +90.0),
96 ],
97 method = "L-BFGS-B",
98 options = {
99 "disp" : debug,
100 "maxiter" : nIter,
101 },
102 tol = conv / RESOLUTION_OF_EARTH,
103 )
104 if not ans.success:
105 print(lons)
106 print(lats)
107 print(ans)
108 raise Exception("failed to converge") from None
110 # Update the middle location ...
111 midLon = ans.x[0] # [°]
112 midLat = ans.x[1] # [°]
114 # Find the maximum Geodesic distance from the middle to any location ...
115 maxDist = max_dist(
116 lons,
117 lats,
118 midLon,
119 midLat,
120 eps = eps,
121 nIter = nIter,
122 space = "GeodesicSpace",
123 ) # [m]
125 if debug:
126 print(f"INFO: The middle is finally ({midLon:.6f}°, {midLat:.6f}°) and the maximum Geodesic distance is finally {0.001 * maxDist:,.1f} km.")
127 else:
128 # Loop over iterations ...
129 for iIter in range(nIter):
130 # Find the angle towards the minimum maximum Geodesic distance ...
131 minAng = find_min_max_dist_bearing(
132 midLon,
133 midLat,
134 lons,
135 lats,
136 angConv = angConv,
137 angHalfRange = 180.0,
138 debug = debug,
139 dist = conv,
140 eps = eps,
141 first = True,
142 iIter = 0,
143 nAng = nAng,
144 nIter = nIter,
145 space = "GeodesicSpace",
146 startAng = 180.0,
147 ) # [°]
149 if debug:
150 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: Moving middle {0.001 * conv:,.1f} km towards {minAng:.1f}° ...")
152 # Find the new location ...
153 newMidLon, newMidLat, _ = calc_loc_from_loc_and_bearing_and_dist(
154 midLon,
155 midLat,
156 minAng,
157 conv,
158 eps = eps,
159 nIter = nIter,
160 ) # [°], [°]
162 # Find the maximum Geodesic distance from the middle to any location ...
163 newMaxDist = max_dist(
164 lons,
165 lats,
166 newMidLon,
167 newMidLat,
168 eps = eps,
169 nIter = nIter,
170 space = "GeodesicSpace",
171 ) # [m]
173 # Stop iterating if the answer isn't getting any better ...
174 if newMaxDist > maxDist:
175 if debug:
176 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: The middle is finally ({midLon:.6f}°, {midLat:.6f}°) and the maximum Geodesic distance is finally {0.001 * maxDist:,.1f} km.")
177 break
179 # Update values ...
180 maxDist = newMaxDist # [m]
181 midLon = newMidLon # [°]
182 midLat = newMidLat # [°]
184 # Stop if the end of the loop has been reached but the answer has
185 # not converged ...
186 if iIter == nIter - 1:
187 raise Exception(f"failed to converge; the middle is currently ({midLon:.6f}°, {midLat:.6f}°); nIter = {nIter:,d}") from None
189 if debug:
190 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: The middle is now ({midLon:.6f}°, {midLat:.6f}°) and the maximum Geodesic distance is now {0.001 * maxDist:,.1f} km.")
192 # **************************************************************************
194 # Check if a padding needs to be added ...
195 if pad > 0.0:
196 # Add padding ...
197 maxDist += pad # [m]
199 if debug:
200 print(f"INFO: Maximum (padded) Geodesic distance is finally {0.001 * maxDist:,.1f} km.")
202 # Return answer ...
203 if iRefine == nRefine - 1:
204 return midLon, midLat, maxDist
205 return find_middle_of_locs_geodesicCircle(
206 lons,
207 lats,
208 angConv = angConv,
209 conv = 0.5 * conv,
210 debug = debug,
211 eps = eps,
212 iRefine = iRefine + 1,
213 midLat = midLat,
214 midLon = midLon,
215 nAng = nAng,
216 nIter = nIter,
217 nRefine = nRefine,
218 pad = pad,
219 useSciPy = useSciPy,
220 )