Coverage for pyguymer3/geo/find_min_max_dist_bearing.py: 83%
42 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_min_max_dist_bearing(
5 midLon,
6 midLat,
7 lons,
8 lats,
9 /,
10 *,
11 angConv = 0.1,
12 angHalfRange = 180.0,
13 debug = __debug__,
14 dist = 1000.0,
15 eps = 1.0e-12,
16 first = True,
17 iIter = 0,
18 nAng = 9,
19 nIter = 100,
20 space = "EuclideanSpace",
21 startAng = 180.0,
22):
23 """Find the bearing which points towards the minimum maximum distance to
24 some locations
26 This function finds the bearing which points towards the minimum maximum
27 distance (in either Euclidean space or Geodesic space) to some locations
28 around a middle location.
30 Parameters
31 ----------
32 midLon : float
33 the middle longitude (in degrees)
34 midLat : float
35 the middle latitude (in degrees)
36 lons : numpy.ndarray
37 the longitudes (in degrees)
38 lats : numpy.ndarray
39 the latitudes (in degrees)
40 angConv : float, optional
41 the angle change which classifies as converged (in degrees)
42 angHalfRange : float, optional
43 the angle either side of the starting angle to search over (in degrees)
44 debug : bool, optional
45 print debug messages
46 dist : float, optional
47 the distance around the middle location to search over (in degrees or
48 metres)
49 eps : float, optional
50 the tolerance of the Vincenty formula iterations
51 first : bool, optional
52 flag whether this is the first call of an interative/recursive sequence
53 (if the first call then this function just returns the angle with the
54 minimum maximum distance; if not the first call then fit a polynomial
55 degree 2 to the values and return the calculated root)
56 iIter : int, optional
57 the current iteration
58 nAng : int, optional
59 the number of angles around the middle location to search over
60 nIter : int, optional
61 the maximum number of iterations (particularly the Vincenty formula)
62 space : str, optional
63 the geometric space to perform the distance calculation in (either
64 "EuclideanSpace" or "GeodesicSpace")
65 startAng : float, optional
66 the starting angle to search over (in degrees)
68 Returns
69 -------
70 rootAng : float
71 the angle which points towards the minimum maximum distance to some
72 locations (in degrees)
74 Notes
75 -----
76 Copyright 2017 Thomas Guymer [1]_
78 References
79 ----------
80 .. [1] PyGuymer3, https://github.com/Guymer/PyGuymer3
81 """
83 # Import standard modules ...
84 import math
86 # Import special modules ...
87 try:
88 import numpy
89 except:
90 raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None
92 # Import sub-functions ...
93 from .calc_loc_from_loc_and_bearing_and_dist import calc_loc_from_loc_and_bearing_and_dist
94 from .max_dist import max_dist
96 # **************************************************************************
98 # Check arguments ...
99 assert isinstance(lons, numpy.ndarray), "\"lons\" is not a NumPy array"
100 assert isinstance(lats, numpy.ndarray), "\"lats\" is not a NumPy array"
101 if iIter == nIter - 1:
102 raise Exception(f"failed to converge; the middle is currently ({midLon:.6f}°, {midLat:.6f}°); nIter = {nIter:,d}") from None
104 # **************************************************************************
106 if debug:
107 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: The middle is now ({midLon:.6f}°, {midLat:.6f}°) and the minimum maximum distance bearing is now {startAng:.6f}°.")
109 # **************************************************************************
111 # Make fake angular axis ...
112 # NOTE: Taking care not to have the duplicate entries of 0° and 360°.
113 if first:
114 fakeAngs = numpy.linspace(
115 startAng - angHalfRange,
116 startAng + angHalfRange,
117 dtype = numpy.float64,
118 endpoint = False,
119 num = nAng - 1,
120 ) # [°]
121 else:
122 fakeAngs = numpy.linspace(
123 startAng - angHalfRange,
124 startAng + angHalfRange,
125 dtype = numpy.float64,
126 endpoint = True,
127 num = nAng,
128 ) # [°]
130 # **************************************************************************
132 # Initialize location arrays ...
133 angLons = numpy.zeros(
134 fakeAngs.size,
135 dtype = numpy.float64,
136 ) # [°]
137 angLats = numpy.zeros(
138 fakeAngs.size,
139 dtype = numpy.float64,
140 ) # [°]
142 # Loop over angles ...
143 for iAng in range(fakeAngs.size):
144 # Check what space the user wants ...
145 match space:
146 case "EuclideanSpace":
147 # Check arguments ...
148 assert eps is None, "\"eps\" is not None but \"space\" is \"EuclideanSpace\""
150 # Populate location arrays ...
151 # NOTE: Taking care not to stray outside of 0° and 360°.
152 angLons[iAng] = midLon + dist * math.sin(math.radians(fakeAngs[iAng])) # [°]
153 angLats[iAng] = midLat + dist * math.cos(math.radians(fakeAngs[iAng])) # [°]
154 case "GeodesicSpace":
155 # Populate location arrays ...
156 # NOTE: Taking care not to stray outside of 0° and 360°.
157 angLons[iAng], angLats[iAng], _ = calc_loc_from_loc_and_bearing_and_dist(
158 midLon,
159 midLat,
160 (fakeAngs[iAng] + 360.0) % 360.0,
161 dist,
162 eps = eps,
163 nIter = nIter,
164 ) # [°], [°]
165 case _:
166 # Crash ...
167 raise ValueError(f"\"space\" is an unexpected value ({repr(space)})") from None
169 # **************************************************************************
171 # Initialize distance array ...
172 maxDists = numpy.zeros(
173 fakeAngs.size,
174 dtype = numpy.float64,
175 ) # [°] or [m]
177 # Populate distance array ...
178 for iAng in range(fakeAngs.size):
179 maxDists[iAng] = max_dist(
180 lons,
181 lats,
182 angLons[iAng],
183 angLats[iAng],
184 eps = eps,
185 nIter = None if space == "EuclideanSpace" else nIter,
186 space = space,
187 ) # [°] or [m]
189 # **************************************************************************
191 # Check if this is the first time that a solution has been found ...
192 if first:
193 # Find angle with minimum maximum distance ...
194 iAng = maxDists.argmin() # [#]
196 # Return answer ...
197 return find_min_max_dist_bearing(
198 midLon,
199 midLat,
200 lons,
201 lats,
202 angConv = angConv,
203 angHalfRange = 0.5 * angHalfRange,
204 debug = debug,
205 dist = dist,
206 eps = eps,
207 first = False,
208 iIter = iIter + 1,
209 nAng = nAng,
210 nIter = nIter,
211 space = space,
212 startAng = (fakeAngs[iAng] + 360.0) % 360.0,
213 )
215 # Fit a polynomial degree 2 to the values and find the best angle ...
216 eqn = numpy.polynomial.Polynomial.fit(
217 fakeAngs,
218 maxDists,
219 deg = 2,
220 )
221 bestAng = eqn.deriv().roots()[0] # [°]
223 # Check if the answer is converged ...
224 if abs(startAng - bestAng) <= angConv:
225 if debug:
226 print(f"INFO: #{iIter + 1:,d}/{nIter:,d}: The middle is finally ({midLon:.6f}°, {midLat:.6f}°) and the minimum maximum distance bearing is finally {startAng:.6f}°.")
227 return bestAng
229 # Return answer ...
230 return find_min_max_dist_bearing(
231 midLon,
232 midLat,
233 lons,
234 lats,
235 angConv = angConv,
236 angHalfRange = 0.5 * angHalfRange,
237 debug = debug,
238 dist = dist,
239 eps = eps,
240 first = False,
241 iIter = iIter + 1,
242 nAng = nAng,
243 nIter = nIter,
244 space = space,
245 startAng = (bestAng + 360.0) % 360.0,
246 )