Coverage for pyguymer3/geo/_points2polys.py: 78%
146 statements
« prev ^ index » next coverage.py v7.10.3, created at 2025-08-16 08:31 +0000
1#!/usr/bin/env python3
3# Define function ...
4def _points2polys(
5 point,
6 points,
7 /,
8 *,
9 debug = __debug__,
10 huge = False,
11 prefix = ".",
12 tol = 1.0e-10,
13):
14 """Convert a buffered point to a list of Polygons
16 This function reads in a coordinate that exists on the surface of the Earth,
17 and an array of coordinates that are the counter-clockwise ring around the
18 coordinate buffered by a constant distance, and returns a list of Polygons
19 which describes the buffer.
21 Parameters
22 ----------
23 point : numpy.ndarray
24 the (2) array of (lon,lat) coordinate (in degrees)
25 points : numpy.ndarray
26 the (nAng, 2) array of (lon,lat) coordinates around the (lon,lat)
27 coordinate (in degrees)
28 debug : bool, optional
29 print debug messages
30 huge : bool, optional
31 if the buffering distance was huge then the points can be turned into a
32 Polygon very easily (as they will definitely cross the [anti-]meridian
33 and a Pole)
34 prefix : str, optional
35 change the name of the output debugging CSVs
36 tol : float, optional
37 the Euclidean distance that defines two points as being the same (in
38 degrees)
40 Returns
41 -------
42 polys : list of shapely.geometry.polygon.Polygon
43 the buffered (lon,lat) coordinate
45 Notes
46 -----
47 According to the `Shapely documentation for the LinearRings
48 <https://shapely.readthedocs.io/en/stable/manual.html#linearrings>`_ , a
49 LinearRing may not cross itself and may not touch itself at a single point.
51 Copyright 2017 Thomas Guymer [1]_
53 References
54 ----------
55 .. [1] PyGuymer3, https://github.com/Guymer/PyGuymer3
56 """
58 # Import special modules ...
59 try:
60 import numpy
61 except:
62 raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None
63 try:
64 import shapely
65 import shapely.geometry
66 except:
67 raise Exception("\"shapely\" is not installed; run \"pip install --user Shapely\"") from None
69 # Import sub-functions ...
70 from .check import check
71 from .wrapLongitude import wrapLongitude
72 from ..interpolate import interpolate
74 # **************************************************************************
76 # Check arguments ...
77 assert isinstance(point, numpy.ndarray), "\"point\" is not a NumPy array"
78 assert isinstance(points, numpy.ndarray), "\"points\" is not a NumPy array"
79 if huge:
80 if debug:
81 print("DEBUG: overruling \"huge\" flag because the feature is currently broken (as of 8/Aug/2025).")
82 huge = False
84 # **************************************************************************
86 # Create short-hands ...
87 nAng = points.shape[0]
88 mang = (nAng - 1) // 2
90 # Determine if the ring crosses the North Pole ...
91 if (wrapLongitude(point[0]) - tol) < wrapLongitude(points[0, 0]) < (wrapLongitude(point[0]) + tol):
92 crossNorthPole = False
93 else:
94 crossNorthPole = True
96 # Determine if the ring crosses the South Pole ...
97 if (wrapLongitude(point[0]) - tol) < wrapLongitude(points[mang, 0]) < (wrapLongitude(point[0]) + tol):
98 crossSouthPole = False
99 else:
100 crossSouthPole = True
102 # Determine if the ring crosses both poles ...
103 crossBothPoles = crossNorthPole and crossSouthPole
105 # **************************************************************************
107 # Check if the user promises that the ring is huge ...
108 if huge and not crossBothPoles:
109 # Find the keys that index the points from West to East (taking care not
110 # to add a duplicate point) ...
111 keys = points[:-1, 0].argsort()
113 # Find the latitude of the [anti-]meridian crossing ...
114 latCross = interpolate(
115 points[keys[0], 0],
116 points[keys[-1], 0] - 360.0,
117 points[keys[0], 1],
118 points[keys[-1], 1],
119 -180.0,
120 ) # [°]
122 # Initialize list which will hold the ring ...
123 ring = []
125 # Add point before the crossing ...
126 ring.append((-180.0, latCross))
128 # Loop over points ...
129 for iang in range(nAng - 1):
130 # Add point ...
131 ring.append((points[keys[iang], 0], points[keys[iang], 1]))
133 # Add point after the crossing ...
134 ring.append((+180.0, latCross))
136 # Clean up ...
137 del keys
139 # Determine if the ring is around a point in the Northern hemisphere ...
140 if point[1] > 0.0:
141 # Add points around the North Pole ...
142 ring.append((+180.0, +90.0))
143 ring.append((-180.0, +90.0))
144 else:
145 # Add points around the South Pole ...
146 ring.append((+180.0, -90.0))
147 ring.append((-180.0, -90.0))
149 # Make a LinearRing out of the ring ...
150 line = shapely.geometry.polygon.LinearRing(ring)
151 if debug:
152 check(line, prefix = prefix)
154 # Clean up ...
155 del ring
157 # Make a correctly oriented Polygon out of this LinearRing ...
158 poly = shapely.geometry.polygon.orient(shapely.geometry.polygon.Polygon(line))
159 if debug:
160 check(poly, prefix = prefix)
162 # Clean up ...
163 del line
165 # NOTE: Do not call "fillin()" on the Polygon. If the user is calling
166 # this function themselves, then they can also call "fillin()"
167 # themselves. If this function is being called by
168 # "buffer_CoordinateSequence()" then the result of
169 # "shapely.ops.unary_union().simplify()" can be filled in instead
170 # in that function.
172 # Return answer ...
173 return [poly]
175 # **************************************************************************
177 # Calculate the Euclidean distances between each point on the ring ...
178 dists = numpy.hypot(numpy.diff(points[:, 0]), numpy.diff(points[:, 1])) # [°]
180 # Initialize list which will hold the two rings and set the index to point
181 # to the first ring ...
182 rings = [[], []]
183 iring = 0
185 # Populate the start of the ring (taking care not to add duplicate points) ...
186 if crossNorthPole and not crossSouthPole:
187 point = (points[0, 0], +90.0)
188 if point not in rings[iring]:
189 rings[iring].append(point)
190 if crossSouthPole and not crossNorthPole:
191 point = (points[0, 0], -90.0)
192 if point not in rings[iring]:
193 rings[iring].append(point)
194 point = (points[0, 0], points[0, 1])
195 if point not in rings[iring]:
196 rings[iring].append(point)
198 # Loop over angles ...
199 for iang in range(nAng - 1):
200 # Check if the ring has crossed the [anti-]meridian ...
201 if dists[iang] >= 180.0:
202 # Check if it was the anti-meridian or the meridian ...
203 if points[iang + 1, 0] < points[iang, 0]:
204 # Set the longitude and find the latitude of the crossing ...
205 lonCross = +180.0 # [°]
206 latCross = interpolate(
207 points[iang, 0],
208 points[iang + 1, 0] + 360.0,
209 points[iang, 1],
210 points[iang + 1, 1],
211 lonCross,
212 ) # [°]
213 else:
214 # Set the longitude and find the latitude of the crossing ...
215 lonCross = -180.0 # [°]
216 latCross = interpolate(
217 points[iang, 0],
218 points[iang + 1, 0] - 360.0,
219 points[iang, 1],
220 points[iang + 1, 1],
221 lonCross,
222 ) # [°]
224 # Add points before the crossing (taking care not to add duplicate
225 # points) ...
226 point = (lonCross, latCross)
227 if point not in rings[iring]:
228 rings[iring].append(point)
229 if crossNorthPole and not crossSouthPole:
230 point = (lonCross, +90.0)
231 if point not in rings[iring]:
232 rings[iring].append(point)
233 if crossSouthPole and not crossNorthPole:
234 point = (lonCross, -90.0)
235 if point not in rings[iring]:
236 rings[iring].append(point)
238 # Switch to the other ring ...
239 iring = 1 - iring
241 # Add points after the crossing (taking care not to add duplicate
242 # points) ...
243 if crossNorthPole and not crossSouthPole:
244 point = (-lonCross, +90.0)
245 if point not in rings[iring]:
246 rings[iring].append(point)
247 if crossSouthPole and not crossNorthPole:
248 point = (-lonCross, -90.0)
249 if point not in rings[iring]:
250 rings[iring].append(point)
251 point = (-lonCross, latCross)
252 if point not in rings[iring]:
253 rings[iring].append(point)
255 # Add point (taking care not to add duplicate points) ...
256 point = (points[iang + 1, 0], points[iang + 1, 1])
257 if point not in rings[iring]:
258 rings[iring].append(point)
260 # Populate the end of the ring (taking care not to add duplicate points) ...
261 if crossNorthPole and not crossSouthPole:
262 point = (points[-1, 0], +90.0)
263 if point not in rings[iring]:
264 rings[iring].append(point)
265 if crossSouthPole and not crossNorthPole:
266 point = (points[-1, 0], -90.0)
267 if point not in rings[iring]:
268 rings[iring].append(point)
270 # Clean up ...
271 del dists
273 # **************************************************************************
275 # Initialize list which will hold the two Polygons ...
276 polys = []
278 # Loop over rings ...
279 for ring in rings:
280 # Skip this ring if it does not have any points ...
281 if len(ring) == 0:
282 continue
284 # Make a NumPy array of this ring ...
285 tmpRing = numpy.array(ring) # [°]
287 # Skip this ring if it has zero width or zero height ...
288 if abs(tmpRing[:, 0].max() - tmpRing[:, 0].min()) < tol:
289 continue
290 if abs(tmpRing[:, 1].max() - tmpRing[:, 1].min()) < tol:
291 continue
293 # Find the Euclidean distance between each point on the ring ...
294 tmpDist = numpy.hypot(
295 numpy.diff(tmpRing[:, 0]),
296 numpy.diff(tmpRing[:, 1]),
297 ) # [°]
299 # Clean up ...
300 del tmpRing
302 # Make a cleaned copy of the ring (the above if-tests for duplicate
303 # points may fail due to the floating-point precision being better than
304 # the floating-point accuracy in Vincenty's formulae) ...
305 cleanedRing = []
306 cleanedRing.append(ring[0])
307 for i in range(tmpDist.size):
308 if tmpDist[i] < tol:
309 continue
310 cleanedRing.append(ring[i + 1])
312 # Clean up ...
313 del tmpDist
315 # Do one final check for if the start and end points are very close
316 # (because if they are not numerically identical then Shapely will add
317 # the first point to the end of the list to close it, which may tangle
318 # the subsequent LinearRing) ...
319 tmpDist = numpy.hypot(
320 cleanedRing[0][0] - cleanedRing[-1][0],
321 cleanedRing[0][1] - cleanedRing[-1][1],
322 ) # [°]
323 if tmpDist < tol:
324 cleanedRing.pop()
326 # Make a LinearRing out of this (cleaned) ring ...
327 line = shapely.geometry.polygon.LinearRing(cleanedRing)
328 if debug:
329 check(line, prefix = prefix)
331 # Clean up ...
332 del cleanedRing
334 # Make a correctly oriented Polygon out of this LinearRing ...
335 poly = shapely.geometry.polygon.orient(shapely.geometry.polygon.Polygon(line))
336 if debug:
337 check(poly, prefix = prefix)
339 # Clean up ...
340 del line
342 # NOTE: Do not call "fillin()" on the Polygon. If the user is calling
343 # this function themselves, then they can also call "fillin()"
344 # themselves. If this function is being called by
345 # "buffer_CoordinateSequence()" then the result of
346 # "shapely.ops.unary_union().simplify()" can be filled in instead
347 # in that function.
349 # Append Polygon to the list ...
350 polys.append(poly)
352 # **************************************************************************
354 # Check if the two Polygons are holes in a larger Polygon ...
355 if crossBothPoles:
356 # Make a correctly oriented Polygon of the planet ...
357 earth = shapely.geometry.polygon.orient(
358 shapely.geometry.polygon.Polygon(
359 shapely.geometry.polygon.LinearRing(
360 [
361 (-180.0, 90.0),
362 (+180.0, 90.0),
363 (+180.0, -90.0),
364 (-180.0, -90.0),
365 (-180.0, 90.0),
366 ]
367 )
368 )
369 )
370 if debug:
371 check(earth, prefix = prefix)
373 # Loop over Polygons ...
374 for poly in polys:
375 # Subtract this Polygon from the planet ...
376 earth = earth.difference(poly)
378 # Clean up ...
379 del polys
381 # Return answer ...
382 return [earth]
384 # Return answer ...
385 return polys