Coverage for pyguymer3/geo/bufferSrc/buffer_CoordinateSequence.py: 62%
63 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 buffer_CoordinateSequence(
5 coords,
6 dist,
7 /,
8 *,
9 debug = __debug__,
10 eps = 1.0e-12,
11 fill = 1.0,
12 fillSpace = "EuclideanSpace",
13 nAng = 9,
14 nIter = 100,
15 prefix = ".",
16 ramLimit = 1073741824,
17 simp = 0.1,
18 tol = 1.0e-10,
19):
20 """Buffer a CoordinateSequence
22 This function reads in a CoordinateSequence that exists on the surface of
23 the Earth and returns a [Multi]Polygon of the same CoordinateSequence
24 buffered by a constant distance (in metres).
26 Parameters
27 ----------
28 coords : shapely.coords.CoordinateSequence
29 the CoordinateSequence
30 dist : float
31 the Geodesic distance to buffer each point within the CoordinateSequence
32 by (in metres)
33 debug : bool, optional
34 print debug messages
35 eps : float, optional
36 the tolerance of the Vincenty formula iterations
37 fill : float, optional
38 the Euclidean or Geodesic distance to fill in between each point within
39 the shapes by (in degrees or metres)
40 fillSpace : str, optional
41 the geometric space to perform the filling in (either "EuclideanSpace"
42 or "GeodesicSpace")
43 nAng : int, optional
44 the number of angles around each point within the CoordinateSequence
45 that are calculated when buffering
46 nIter : int, optional
47 the maximum number of iterations (particularly the Vincenty formula)
48 prefix : str, optional
49 change the name of the output debugging CSVs
50 ramLimit : int, optional
51 the maximum RAM usage of each "large" array (in bytes)
52 simp : float, optional
53 how much the final [Multi]Polygons is simplified by; negative values
54 disable simplification (in degrees)
55 tol : float, optional
56 the Euclidean distance that defines two points as being the same (in
57 degrees)
59 Returns
60 -------
61 buffs : shapely.geometry.polygon.Polygon, shapely.geometry.multipolygon.MultiPolygon
62 the buffered CoordinateSequence
64 Notes
65 -----
66 According to the `Shapely documentation for the method object.buffer()
67 <https://shapely.readthedocs.io/en/stable/manual.html#object.buffer>`_ :
69 "Passed a distance of 0, buffer() can sometimes be used to "clean"
70 self-touching or self-crossing polygons such as the classic "bowtie".
71 Users have reported that very small distance values sometimes produce
72 cleaner results than 0. Your mileage may vary when cleaning surfaces."
74 According to the `Shapely documentation for the function
75 shapely.geometry.polygon.orient()
76 <https://shapely.readthedocs.io/en/stable/manual.html#shapely.geometry.polygon.orient>`_ :
78 "A sign of 1.0 means that the coordinates of the product's exterior ring
79 will be oriented counter-clockwise."
81 Copyright 2017 Thomas Guymer [1]_
83 References
84 ----------
85 .. [1] PyGuymer3, https://github.com/Guymer/PyGuymer3
86 """
88 # Import special modules ...
89 try:
90 import numpy
91 except:
92 raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None
93 try:
94 import shapely
95 import shapely.geometry
96 import shapely.ops
97 except:
98 raise Exception("\"shapely\" is not installed; run \"pip install --user Shapely\"") from None
100 # Import sub-functions ...
101 from .._buffer_points_crudely import _buffer_points_crudely
102 from .._points2polys import _points2polys
103 from ..check import check
104 from ..clean import clean
105 from ..fillin import fillin
106 from ...consts import CIRCUMFERENCE_OF_EARTH
107 try:
108 from ...f90 import funcs
109 if debug:
110 print("INFO: Will find the rings using FORTRAN.")
111 fortran = True
112 except ModuleNotFoundError:
113 if debug:
114 print("INFO: Will find the rings using Python (did not find FORTRAN module).")
115 fortran = False
116 except ImportError:
117 if debug:
118 print("INFO: Will find the rings using Python (error when attempting to import found FORTRAN).")
119 fortran = False
121 # **************************************************************************
123 # Check argument ...
124 assert isinstance(coords, shapely.coords.CoordinateSequence), "\"coords\" is not a CoordinateSequence"
125 if debug:
126 check(coords, prefix = prefix)
128 # Check inputs ...
129 assert dist >= 10.0, f"the buffering distance is too small ({dist:,.1f}m < {10.0:,.1f}m)"
130 assert dist <= 0.5 * CIRCUMFERENCE_OF_EARTH, f"the buffering distance is too large ({dist:,.1f}m > {0.5 * CIRCUMFERENCE_OF_EARTH:,.1f}m)"
131 assert nAng >= 9, f"the number of angles is too small ({nAng:,d} < {9:,d})"
132 assert nAng % 2 == 1, f"the number of angles is even ({nAng:,d})"
133 assert (nAng - 1) % 4 == 0, f"the number of angles is not 4n+1 ({nAng:,d})"
135 # **************************************************************************
136 # Step 1: Convert the CoordinateSequence to a NumPy array of the original #
137 # points #
138 # **************************************************************************
140 # Clean the input ...
141 coords = clean(
142 coords,
143 debug = debug,
144 prefix = prefix,
145 tol = tol,
146 ).coords
148 # Check if the user wants to fill in the CoordinateSequence ...
149 if fill > 0.0 and len(coords) > 1:
150 # Convert the filled in CoordinateSequence to a NumPy array ...
151 points1 = numpy.array(
152 fillin(
153 coords,
154 fill,
155 debug = debug,
156 eps = eps,
157 fillSpace = fillSpace,
158 nIter = nIter,
159 prefix = prefix,
160 ramLimit = ramLimit,
161 tol = tol,
162 ).coords
163 ) # [°]
164 else:
165 # Convert the CoordinateSequence to a NumPy array ...
166 points1 = numpy.array(coords) # [°]
168 # Create short-hand ...
169 npoint = int(points1.shape[0]) # [#]
171 # **************************************************************************
172 # Step 2: Buffer the NumPy array of the original points to get a NumPy #
173 # array of the rings around them #
174 # **************************************************************************
176 # Buffer (in Geodesic space) the CoordinateSequence ...
177 if fortran:
178 points2 = funcs.buffer_points_crudely(
179 points1,
180 dist,
181 nAng,
182 ) # [°]
183 else:
184 points2 = _buffer_points_crudely(
185 points1,
186 dist,
187 nAng,
188 eps = eps,
189 nIter = nIter,
190 ramLimit = ramLimit,
191 ) # [°]
193 # **************************************************************************
194 # Step 3: Convert the NumPy array of the rings around the original points #
195 # to a list of Polygons of the buffered original points #
196 # **************************************************************************
198 # Initialize list of Polygons ...
199 buffs = []
201 # Loop over points ...
202 for ipoint in range(npoint):
203 # Add list of Polygons to list of Polygons ...
204 buffs += _points2polys(
205 points1[ipoint, :],
206 points2[ipoint, :, :],
207 debug = debug,
208 huge = dist > 0.25 * CIRCUMFERENCE_OF_EARTH,
209 prefix = prefix,
210 tol = tol,
211 )
213 # Clean up ...
214 del points1, points2
216 # **************************************************************************
217 # Step 4: Create a single [Multi]Polygon that is the union of all of the #
218 # Polygons #
219 # **************************************************************************
221 # Convert list of Polygons to a (unified) [Multi]Polygon ...
222 buffs = shapely.ops.unary_union(buffs).simplify(tol)
223 if debug:
224 check(buffs, prefix = prefix)
226 # Check if the user wants to fill in the [Multi]Polygon ...
227 # NOTE: This is only needed because the "shapely.ops.unary_union()" call
228 # above includes a "simplify()".
229 if simp < 0.0 < fill:
230 # Fill in [Multi]Polygon ...
231 buffs = fillin(
232 buffs,
233 fill,
234 debug = debug,
235 eps = eps,
236 fillSpace = fillSpace,
237 nIter = nIter,
238 prefix = prefix,
239 ramLimit = ramLimit,
240 tol = tol,
241 )
242 if debug:
243 check(buffs, prefix = prefix)
245 # Check if the user wants to simplify the [Multi]Polygon ...
246 # NOTE: This is only needed because the "shapely.ops.unary_union()" call
247 # above might allow more simplification.
248 if simp > 0.0:
249 # Simplify [Multi]Polygon ...
250 buffsSimp = buffs.simplify(simp)
251 if debug:
252 check(buffsSimp, prefix = prefix)
254 # Return simplified answer ...
255 return buffsSimp
257 # Return answer ...
258 return buffs