Coverage for pyguymer3/geo/find_middle_of_locsSrc/find_middle_of_locs_geodesicBox.py: 63%
81 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_geodesicBox(
5 lons,
6 lats,
7 /,
8 *,
9 conv = 1.0e3,
10 debug = __debug__,
11 eps = 1.0e-12,
12 iRefine = 0,
13 midLat = None,
14 midLon = None,
15 nIter = 100,
16 nRefine = 1,
17 pad = 10.0e3,
18):
19 """Find the middle of some locations such that: a) the Geodesic distance to
20 the most Northern point is the same as the Geodesic distance to the most
21 Southern point; and b) the Geodesic distance to the most Eastern point is
22 the same as the Geodesic distance to the most Western point.
23 """
25 # Import special modules ...
26 try:
27 import numpy
28 except:
29 raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None
31 # Import sub-functions ...
32 from .find_middle_of_locs_euclideanBox import find_middle_of_locs_euclideanBox
33 from ..calc_dist_between_two_locs import calc_dist_between_two_locs
34 from ..calc_loc_from_loc_and_bearing_and_dist import calc_loc_from_loc_and_bearing_and_dist
35 from ..max_dist import max_dist
37 # **************************************************************************
39 # Check arguments ...
40 assert isinstance(lons, numpy.ndarray), "\"lons\" is not a NumPy array"
41 assert isinstance(lats, numpy.ndarray), "\"lats\" is not a NumPy array"
43 # **************************************************************************
45 # Calculate the middle of the Euclidean bounding box to use as a decent
46 # starting guess (if the user has not provided one) ...
47 if midLon is None or midLat is None:
48 midLon, midLat, _ = find_middle_of_locs_euclideanBox(
49 lons,
50 lats,
51 debug = debug,
52 pad = -1.0,
53 ) # [°], [°]
55 # Find the maximum Geodesic distance from the middle to any location ...
56 maxDist = max_dist(
57 lons,
58 lats,
59 midLon,
60 midLat,
61 eps = eps,
62 nIter = nIter,
63 space = "GeodesicSpace",
64 ) # [m]
66 if debug:
67 print(f"INFO: The initial middle is ({midLon:.6f}°, {midLat:.6f}°) and the initial maximum Geodesic distance is {0.001 * maxDist:,.1f} km.")
69 # Loop over iterations ...
70 for iIter in range(nIter):
71 # Find the Geodesic distances to the western-most and eastern-most
72 # locations ...
73 maxWestDist = 0.0 # [m]
74 maxEastDist = 0.0 # [m]
75 for iLoc in range(lons.size):
76 if lons[iLoc] < midLon: # NOTE: Location is west of the middle.
77 maxWestDist = max(
78 maxWestDist,
79 calc_dist_between_two_locs(
80 midLon,
81 midLat,
82 lons[iLoc],
83 midLat,
84 eps = eps,
85 nIter = nIter,
86 )[0]
87 ) # [m]
88 elif lons[iLoc] > midLon: # NOTE: Location is east of the middle.
89 maxEastDist = max(
90 maxEastDist,
91 calc_dist_between_two_locs(
92 midLon,
93 midLat,
94 lons[iLoc],
95 midLat,
96 eps = eps,
97 nIter = nIter,
98 )[0]
99 ) # [m]
100 else:
101 pass
103 # Find the Geodesic distances to the southern-most and northern-most
104 # locations ...
105 maxSouthDist = 0.0 # [m]
106 maxNorthDist = 0.0 # [m]
107 for iLoc in range(lons.size):
108 if lats[iLoc] < midLat: # NOTE: Location is south of the middle.
109 maxSouthDist = max(
110 maxSouthDist,
111 calc_dist_between_two_locs(
112 midLon,
113 midLat,
114 midLon,
115 lats[iLoc],
116 eps = eps,
117 nIter = nIter,
118 )[0]
119 ) # [m]
120 elif lats[iLoc] > midLat: # NOTE: Location is north of the middle.
121 maxNorthDist = max(
122 maxNorthDist,
123 calc_dist_between_two_locs(
124 midLon,
125 midLat,
126 midLon,
127 lats[iLoc],
128 eps = eps,
129 nIter = nIter,
130 )[0]
131 ) # [m]
132 else:
133 pass
135 if debug:
136 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: {0.001 * maxWestDist:,.1f} km west ← middle → {0.001 * maxEastDist:,.1f} km east.")
137 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: {0.001 * maxSouthDist:,.1f} km south ↓ middle ↑ {0.001 * maxNorthDist:,.1f} km north.")
139 # Initialize flag ...
140 moved = False
142 # Check if the middle needs moving west/east ...
143 if 0.5 * (maxWestDist - maxEastDist) > conv:
144 if debug:
145 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: Moving middle {0.001 * 0.5 * (maxWestDist - maxEastDist):,.1f} km towards the west ...")
146 midLon, midLat, _ = calc_loc_from_loc_and_bearing_and_dist(
147 midLon,
148 midLat,
149 270.0,
150 0.5 * (maxWestDist - maxEastDist),
151 eps = eps,
152 nIter = nIter,
153 ) # [°], [°]
154 if debug:
155 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: The middle is now ({midLon:.6f}°, {midLat:.6f}°).")
156 moved = True
157 elif 0.5 * (maxEastDist - maxWestDist) > conv:
158 if debug:
159 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: Moving middle {0.001 * 0.5 * (maxEastDist - maxWestDist):,.1f} km towards the east ...")
160 midLon, midLat, _ = calc_loc_from_loc_and_bearing_and_dist(
161 midLon,
162 midLat,
163 90.0,
164 0.5 * (maxEastDist - maxWestDist),
165 eps = eps,
166 nIter = nIter,
167 ) # [°], [°]
168 if debug:
169 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: The middle is now ({midLon:.6f}°, {midLat:.6f}°).")
170 moved = True
171 else:
172 pass
174 # Check if the middle needs moving south/north ...
175 if 0.5 * (maxSouthDist - maxNorthDist) > conv:
176 if debug:
177 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: Moving middle {0.001 * 0.5 * (maxSouthDist - maxNorthDist):,.1f} km towards the south ...")
178 midLon, midLat, _ = calc_loc_from_loc_and_bearing_and_dist(
179 midLon,
180 midLat,
181 180.0,
182 0.5 * (maxSouthDist - maxNorthDist),
183 eps = eps,
184 nIter = nIter,
185 ) # [°], [°]
186 if debug:
187 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: The middle is now ({midLon:.6f}°, {midLat:.6f}°).")
188 moved = True
189 elif 0.5 * (maxNorthDist - maxSouthDist) > conv:
190 if debug:
191 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: Moving middle {0.001 * 0.5 * (maxNorthDist - maxSouthDist):,.1f} km towards the north ...")
192 midLon, midLat, _ = calc_loc_from_loc_and_bearing_and_dist(
193 midLon,
194 midLat,
195 0.0,
196 0.5 * (maxNorthDist - maxSouthDist),
197 eps = eps,
198 nIter = nIter,
199 ) # [°], [°]
200 if debug:
201 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: The middle is now ({midLon:.6f}°, {midLat:.6f}°).")
202 moved = True
203 else:
204 pass
206 # Check if we can stop iterating ...
207 if not moved:
208 break
210 # Stop if the end of the loop has been reached but the answer has not
211 # converged ...
212 if iIter == nIter - 1:
213 raise Exception(f"failed to converge; the middle is currently ({midLon:.6f}°, {midLat:.6f}°); nIter = {nIter:,d}") from None
215 # Find the maximum Geodesic distance from the middle to any location ...
216 maxDist = max_dist(
217 lons,
218 lats,
219 midLon,
220 midLat,
221 eps = eps,
222 nIter = nIter,
223 space = "GeodesicSpace",
224 ) # [m]
226 if debug:
227 print(f"INFO: Maximum Geodesic distance is {0.001 * maxDist:,.1f} km.")
229 # **************************************************************************
231 # Check if a padding needs to be added ...
232 if pad > 0.0:
233 # Add padding ...
234 maxDist += pad # [m]
236 if debug:
237 print(f"INFO: Maximum (padded) Geodesic distance is {0.001 * maxDist:,.1f} km.")
239 # Return answer ...
240 if iRefine == nRefine - 1:
241 return midLon, midLat, maxDist
242 return find_middle_of_locs_geodesicBox(
243 lons,
244 lats,
245 conv = 0.5 * conv,
246 debug = debug,
247 eps = eps,
248 iRefine = iRefine + 1,
249 midLat = midLat,
250 midLon = midLon,
251 nIter = nIter,
252 nRefine = nRefine,
253 pad = pad,
254 )