Coverage for pyguymer3/geo/fillinSrc/fillin_CoordinateSequence.py: 78%
64 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 fillin_CoordinateSequence(
5 coords,
6 fill,
7 /,
8 *,
9 debug = __debug__,
10 eps = 1.0e-12,
11 fillSpace = "EuclideanSpace",
12 nIter = 100,
13 prefix = ".",
14 ramLimit = 1073741824,
15):
16 """Fill in a CoordinateSequence
18 This function reads in a CoordinateSequence that exists on the surface of
19 the Earth and returns a LineString of the same CoordinateSequence filled in
20 by a constant distance: either in degrees in Euclidean space; or in metres
21 in Geodesic space.
23 Parameters
24 ----------
25 coords : shapely.coords.CoordinateSequence
26 the CoordinateSequence
27 fill : float
28 the Euclidean or Geodesic distance to fill in between each point within
29 the shape by (in degrees or metres)
30 debug : bool, optional
31 print debug messages
32 eps : float, optional
33 the tolerance of the Vincenty formula iterations
34 fillSpace : str, optional
35 the geometric space to perform the filling in (either "EuclideanSpace"
36 or "GeodesicSpace")
37 nIter : int, optional
38 the maximum number of iterations (particularly the Vincenty formula)
39 prefix : str, optional
40 change the name of the output debugging CSVs
41 ramLimit : int, optional
42 the maximum RAM usage of each "large" array (in bytes)
44 Returns
45 -------
46 fills : shapely.geometry.linestring.LineString
47 the filled in CoordinateSequence
49 Notes
50 -----
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 ..calc_dist_between_two_locs import calc_dist_between_two_locs
71 from ..check import check
72 from ..great_circle import great_circle
74 # **************************************************************************
76 # Check argument ...
77 assert isinstance(coords, shapely.coords.CoordinateSequence), "\"coords\" is not a CoordinateSequence"
78 if debug:
79 check(coords, prefix = prefix)
81 # Convert the CoordinateSequence to a NumPy array ...
82 points1 = numpy.array(coords) # [°]
84 # Check what space the user wants to fill in ...
85 match fillSpace:
86 case "EuclideanSpace":
87 # Find the Euclidean distance between each original point ...
88 dr = numpy.hypot(numpy.diff(points1[:, 0]), numpy.diff(points1[:, 1])) # [°]
89 case "GeodesicSpace":
90 # Find the Geodesic distance between each original point ...
91 dr = numpy.zeros(points1.shape[0] - 1, dtype = numpy.float64) # [m]
92 for ipoint in range(dr.size):
93 dr[ipoint], _, _ = calc_dist_between_two_locs(
94 points1[ipoint , 0],
95 points1[ipoint , 1],
96 points1[ipoint + 1, 0],
97 points1[ipoint + 1, 1],
98 eps = eps,
99 nIter = nIter,
100 ) # [m]
101 case _:
102 # Crash ...
103 raise ValueError(f"\"fillSpace\" is an unexpected value ({repr(fillSpace)})") from None
105 # Find the number of filled segments required between each original point ...
106 # NOTE: This is the number of fence panels, not the number of fence posts.
107 ns = (dr / fill).astype(numpy.int32) # [#]
108 numpy.place(ns, ns < 1, 1) # [#]
110 # Clean up ...
111 del dr
113 # Find the total number of filled points required ...
114 # NOTE: This is the total number of fence posts, not the total number of
115 # fence panels.
116 nsTot = ns.sum() + 1 # [#]
118 # Check array size ...
119 if nsTot * 2 * 8 > ramLimit:
120 raise Exception(f"\"points2\" is going to be {nsTot * 2 * 8:,d} bytes, which is larger than {ramLimit:,d} bytes") from None
122 # Create empty array ...
123 points2 = numpy.zeros((nsTot, 2), dtype = numpy.float64) # [°]
125 if debug:
126 print(f"INFO: There are x{points2.shape[0] / points1.shape[0]:,.1f} more points due to filling in.")
128 # Initialize index ...
129 ifill = 0 # [#]
131 # Check what space the user wants to fill in ...
132 match fillSpace:
133 case "EuclideanSpace":
134 # Loop over original points ...
135 for ipoint in range(ns.size):
136 # Check if no filling in is actually needed ...
137 if ns[ipoint] == 1:
138 # Fill in point ...
139 points2[ifill, :] = points1[ipoint, :] # [°]
140 else:
141 # Fill in points ...
142 points2[ifill:ifill + ns[ipoint], 0] = numpy.linspace(
143 points1[ipoint , 0],
144 points1[ipoint + 1, 0],
145 endpoint = False,
146 num = ns[ipoint],
147 ) # [°]
148 points2[ifill:ifill + ns[ipoint], 1] = numpy.linspace(
149 points1[ipoint , 1],
150 points1[ipoint + 1, 1],
151 endpoint = False,
152 num = ns[ipoint],
153 ) # [°]
155 # Increment index ...
156 ifill += ns[ipoint] # [#]
157 case "GeodesicSpace":
158 # Loop over original points ...
159 for ipoint in range(ns.size):
160 # Check if no filling in is actually needed ...
161 if ns[ipoint] == 1:
162 # Fill in point ...
163 points2[ifill, :] = points1[ipoint, :] # [°]
164 else:
165 # Find the great circle connecting the two original points
166 # and convert the LineString to a NumPy array ...
167 arc = great_circle(
168 points1[ipoint , 0],
169 points1[ipoint , 1],
170 points1[ipoint + 1, 0],
171 points1[ipoint + 1, 1],
172 debug = debug,
173 eps = eps,
174 maxdist = None,
175 nIter = nIter,
176 npoint = int(ns[ipoint]) + 1,
177 prefix = prefix,
178 ramLimit = ramLimit,
179 )
180 arcCoords = numpy.array(arc.coords) # [°]
182 # Clean up ...
183 del arc
185 # Fill in points ...
186 points2[ifill:ifill + ns[ipoint], :] = arcCoords[:-1, :] # [°]
188 # Clean up ...
189 del arcCoords
191 # Increment index ...
192 ifill += ns[ipoint] # [#]
193 case _:
194 # Crash ...
195 raise ValueError(f"\"fillSpace\" is an unexpected value ({repr(fillSpace)})") from None
197 # Clean up ...
198 del ns
200 # Fill in last point ...
201 points2[-1, :] = points1[-1, :] # [°]
203 # Clean up ...
204 del points1
206 # Convert array of points to a LineString ...
207 fills = shapely.geometry.linestring.LineString(points2)
208 if debug:
209 check(fills, prefix = prefix)
211 # Clean up ...
212 del points2
214 # Return answer ...
215 return fills