Coverage for pyguymer3/geo/_points2polys.py: 43%
142 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 _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"
80 # **************************************************************************
82 # Create short-hands ...
83 nAng = points.shape[0]
84 mang = (nAng - 1) // 2
86 # Determine if the ring crosses the North Pole ...
87 if (wrapLongitude(point[0]) - tol) < wrapLongitude(points[0, 0]) < (wrapLongitude(point[0]) + tol):
88 crossNorthPole = False
89 else:
90 crossNorthPole = True
92 # Determine if the ring crosses the South Pole ...
93 if (wrapLongitude(point[0]) - tol) < wrapLongitude(points[mang, 0]) < (wrapLongitude(point[0]) + tol):
94 crossSouthPole = False
95 else:
96 crossSouthPole = True
98 # Determine if the ring crosses both poles ...
99 crossBothPoles = crossNorthPole and crossSouthPole
101 # **************************************************************************
103 # Check if the user promises that the ring is huge ...
104 if huge and not crossBothPoles:
105 # Find the keys that index the points from West to East (taking care not
106 # to add a duplicate point) ...
107 keys = points[:-1, 0].argsort()
109 # Find the latitude of the [anti-]meridian crossing ...
110 latCross = interpolate(
111 points[keys[0], 0],
112 points[keys[-1], 0] - 360.0,
113 points[keys[0], 1],
114 points[keys[-1], 1],
115 -180.0,
116 ) # [°]
118 # Initialize list which will hold the ring ...
119 ring = []
121 # Add point before the crossing ...
122 ring.append((-180.0, latCross))
124 # Loop over points ...
125 for iang in range(nAng - 1):
126 # Add point ...
127 ring.append((points[keys[iang], 0], points[keys[iang], 1]))
129 # Add point after the crossing ...
130 ring.append((+180.0, latCross))
132 # Clean up ...
133 del keys
135 # Determine if the ring is around a point in the Northern hemisphere ...
136 if point[1] > 0.0:
137 # Add points around the North Pole ...
138 ring.append((+180.0, +90.0))
139 ring.append((-180.0, +90.0))
140 else:
141 # Add points around the South Pole ...
142 ring.append((+180.0, -90.0))
143 ring.append((-180.0, -90.0))
145 # Make a LinearRing out of the ring ...
146 line = shapely.geometry.polygon.LinearRing(ring)
147 if debug:
148 check(line, prefix = prefix)
150 # Clean up ...
151 del ring
153 # Make a correctly oriented Polygon out of this LinearRing ...
154 poly = shapely.geometry.polygon.orient(shapely.geometry.polygon.Polygon(line))
155 if debug:
156 check(poly, prefix = prefix)
158 # Clean up ...
159 del line
161 # NOTE: Do not call "fillin()" on the Polygon. If the user is calling
162 # this function themselves, then they can also call "fillin()"
163 # themselves. If this function is being called by
164 # "buffer_CoordinateSequence()" then the result of
165 # "shapely.ops.unary_union().simplify()" can be filled in instead
166 # in that function.
168 # Return answer ...
169 return [poly]
171 # **************************************************************************
173 # Calculate the Euclidean distances between each point on the ring ...
174 dists = numpy.hypot(numpy.diff(points[:, 0]), numpy.diff(points[:, 1])) # [°]
176 # Initialize list which will hold the two rings and set the index to point
177 # to the first ring ...
178 rings = [[], []]
179 iring = 0
181 # Populate the start of the ring (taking care not to add duplicate points) ...
182 if crossNorthPole and not crossSouthPole:
183 point = (points[0, 0], +90.0)
184 if point not in rings[iring]:
185 rings[iring].append(point)
186 if crossSouthPole and not crossNorthPole:
187 point = (points[0, 0], -90.0)
188 if point not in rings[iring]:
189 rings[iring].append(point)
190 point = (points[0, 0], points[0, 1])
191 if point not in rings[iring]:
192 rings[iring].append(point)
194 # Loop over angles ...
195 for iang in range(nAng - 1):
196 # Check if the ring has crossed the [anti-]meridian ...
197 if dists[iang] >= 180.0:
198 # Check if it was the anti-meridian or the meridian ...
199 if points[iang + 1, 0] < points[iang, 0]:
200 # Set the longitude and find the latitude of the crossing ...
201 lonCross = +180.0 # [°]
202 latCross = interpolate(
203 points[iang, 0],
204 points[iang + 1, 0] + 360.0,
205 points[iang, 1],
206 points[iang + 1, 1],
207 lonCross,
208 ) # [°]
209 else:
210 # Set the longitude and find the latitude of the crossing ...
211 lonCross = -180.0 # [°]
212 latCross = interpolate(
213 points[iang, 0],
214 points[iang + 1, 0] - 360.0,
215 points[iang, 1],
216 points[iang + 1, 1],
217 lonCross,
218 ) # [°]
220 # Add points before the crossing (taking care not to add duplicate
221 # points) ...
222 point = (lonCross, latCross)
223 if point not in rings[iring]:
224 rings[iring].append(point)
225 if crossNorthPole and not crossSouthPole:
226 point = (lonCross, +90.0)
227 if point not in rings[iring]:
228 rings[iring].append(point)
229 if crossSouthPole and not crossNorthPole:
230 point = (lonCross, -90.0)
231 if point not in rings[iring]:
232 rings[iring].append(point)
234 # Switch to the other ring ...
235 iring = 1 - iring
237 # Add points after the crossing (taking care not to add duplicate
238 # points) ...
239 if crossNorthPole and not crossSouthPole:
240 point = (-lonCross, +90.0)
241 if point not in rings[iring]:
242 rings[iring].append(point)
243 if crossSouthPole and not crossNorthPole:
244 point = (-lonCross, -90.0)
245 if point not in rings[iring]:
246 rings[iring].append(point)
247 point = (-lonCross, latCross)
248 if point not in rings[iring]:
249 rings[iring].append(point)
251 # Add point (taking care not to add duplicate points) ...
252 point = (points[iang + 1, 0], points[iang + 1, 1])
253 if point not in rings[iring]:
254 rings[iring].append(point)
256 # Populate the end of the ring (taking care not to add duplicate points) ...
257 if crossNorthPole and not crossSouthPole:
258 point = (points[-1, 0], +90.0)
259 if point not in rings[iring]:
260 rings[iring].append(point)
261 if crossSouthPole and not crossNorthPole:
262 point = (points[-1, 0], -90.0)
263 if point not in rings[iring]:
264 rings[iring].append(point)
266 # Clean up ...
267 del dists
269 # **************************************************************************
271 # Initialize list which will hold the two Polygons ...
272 polys = []
274 # Loop over rings ...
275 for ring in rings:
276 # Skip this ring if it does not have any points ...
277 if len(ring) == 0:
278 continue
280 # Make a NumPy array of this ring ...
281 tmpRing = numpy.array(ring) # [°]
283 # Skip this ring if it has zero width or zero height ...
284 if abs(tmpRing[:, 0].max() - tmpRing[:, 0].min()) < tol:
285 continue
286 if abs(tmpRing[:, 1].max() - tmpRing[:, 1].min()) < tol:
287 continue
289 # Find the Euclidean distance between each point on the ring ...
290 tmpDist = numpy.hypot(
291 numpy.diff(tmpRing[:, 0]),
292 numpy.diff(tmpRing[:, 1]),
293 ) # [°]
295 # Clean up ...
296 del tmpRing
298 # Make a cleaned copy of the ring (the above if-tests for duplicate
299 # points may fail due to the floating-point precision being better than
300 # the floating-point accuracy in Vincenty's formulae) ...
301 cleanedRing = []
302 cleanedRing.append(ring[0])
303 for i in range(tmpDist.size):
304 if tmpDist[i] < tol:
305 continue
306 cleanedRing.append(ring[i + 1])
308 # Clean up ...
309 del tmpDist
311 # Do one final check for if the start and end points are very close
312 # (because if they are not numerically identical then Shapely will add
313 # the first point to the end of the list to close it, which may tangle
314 # the subsequent LinearRing) ...
315 tmpDist = numpy.hypot(
316 cleanedRing[0][0] - cleanedRing[-1][0],
317 cleanedRing[0][1] - cleanedRing[-1][1],
318 ) # [°]
319 if tmpDist < tol:
320 cleanedRing.pop()
322 # Make a LinearRing out of this (cleaned) ring ...
323 line = shapely.geometry.polygon.LinearRing(cleanedRing)
324 if debug:
325 check(line, prefix = prefix)
327 # Clean up ...
328 del cleanedRing
330 # Make a correctly oriented Polygon out of this LinearRing ...
331 poly = shapely.geometry.polygon.orient(shapely.geometry.polygon.Polygon(line))
332 if debug:
333 check(poly, prefix = prefix)
335 # Clean up ...
336 del line
338 # NOTE: Do not call "fillin()" on the Polygon. If the user is calling
339 # this function themselves, then they can also call "fillin()"
340 # themselves. If this function is being called by
341 # "buffer_CoordinateSequence()" then the result of
342 # "shapely.ops.unary_union().simplify()" can be filled in instead
343 # in that function.
345 # Append Polygon to the list ...
346 polys.append(poly)
348 # **************************************************************************
350 # Check if the two Polygons are holes in a larger Polygon ...
351 if crossBothPoles:
352 # Make a correctly oriented Polygon of the planet ...
353 earth = shapely.geometry.polygon.orient(
354 shapely.geometry.polygon.Polygon(
355 shapely.geometry.polygon.LinearRing(
356 [
357 (-180.0, 90.0),
358 (+180.0, 90.0),
359 (+180.0, -90.0),
360 (-180.0, -90.0),
361 (-180.0, 90.0),
362 ]
363 )
364 )
365 )
366 if debug:
367 check(earth, prefix = prefix)
369 # Loop over Polygons ...
370 for poly in polys:
371 # Subtract this Polygon from the planet ...
372 earth = earth.difference(poly)
374 # Clean up ...
375 del polys
377 # Return answer ...
378 return [earth]
380 # Return answer ...
381 return polys