-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgeom.py
359 lines (264 loc) · 10.6 KB
/
geom.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
"""
Helper functions for working with geometries
"""
import math as _math
import hashlib as _hashlib
import copy as _copy
import arcpy
import arcpy as _arcpy
import arcproapi.errors as _errors # noqa
#################
# Classes first #
#################
class _Shape:
def as_list(self):
"""Return shape polygon as list of lists
Examples:
>>> C = Circle([0,0], 1)
>>> C.as_list
[[1.212,1.111], ...]
"""
return list_from_polygon(self.Polygon)
def as_geometry(self, **kwargs) -> _arcpy.Geometry:
"""
Get instance as a Geometry instance.
This is consumed by the InsertCursor
Args:
kwargs: passed to arcpy.Geometry init
Returns:
arcpy.Geometry: Geometry instance of square points
"""
return _arcpy.Geometry('polygon', _arcpy.Array([self.pt_x1y1, self.pt_x1y2, self.pt_x2y2, self.pt_x2y1, self.pt_x1y1]), **kwargs)
def as_array(self) -> _arcpy.Array:
"""get as an arcpy Array instance
Returns:
arcpy.Array: Array instance of square points
"""
return _arcpy.Array([self.pt_x1y1, self.pt_x1y2, self.pt_x2y2, self.pt_x2y1, self.pt_x1y1])
class Circle(_Shape):
""" Stuff around circles ... to expand
Credit:
QFSW at https://stackoverflow.com/a/42879185/5585800 for the core point derivation
"""
def __init__(self, origin: (list, _arcpy.Point) = (0, 0), radius: float = 1, step_size: float = 0.1):
if isinstance(origin, (list, tuple)):
self.pt_x1y1 = _arcpy.Point(*origin)
else:
self.pt_x1y1 = _copy.copy(origin) # deepcopy fails (ESRI object doesnt support necessary interfaces), this just copies the pointer. We dont manipulate origin anyway
self._radius = radius
self._step_size = step_size
self._radius = radius
positions = []
t = 0
while t < 2 * _math.pi:
positions.append((self._radius * _math.cos(t) + self.pt_x1y1.X, self._radius * _math.sin(t) + self.pt_x1y1.Y))
t += self._step_size
self.Polygon = polygon_from_list(positions)
class Square(_Shape):
"""
Create an arcpy Polygon instance representation of a square
with origin at origin, with sides of length side_length.
Origin is typlically in standard cartesian convention (i.e. x1,y1)
Args:
origin (list, tuple, arcpy.Point): A point like object
side_length (float, int): Side length
Members:
Polygon: The arcpy Polygon instance
pt_x?y?: arcpy point objects, defining the square
Examples:
>>> Sq = Square([0,0], 1)
>>> print(Sq.Polygon.centroid)
<Point (0.50006103515625, 0.50006103515625, #, #)>
"""
def __init__(self, origin: (list, _arcpy.Point), side_length):
if isinstance(origin, (list, tuple)):
self.pt_x1y1 = _arcpy.Point(*origin)
else:
self.pt_x1y1 = _copy.copy(origin) # deepcopy fails (ESRI object doesnt support necessary interfaces), this just copies the pointer. We dont manipulate origin anyway
self.pt_x1y2 = _arcpy.Point(self.pt_x1y1.X, self.pt_x1y1.Y + side_length)
self.pt_x2y2 = _arcpy.Point(self.pt_x1y1.X + side_length, self.pt_x1y1.Y + side_length)
self.pt_x2y1 = _arcpy.Point(self.pt_x1y1.X + side_length, self.pt_x1y1.Y)
self.points: list[_arcpy.Point] = [self.pt_x1y1, self.pt_x1y2, self.pt_x2y2, self.pt_x2y1, self.pt_x1y1]
self.Polygon = _arcpy.Polygon(_arcpy.Array(self.points))
###################
# General methods #
###################
def shape_hash(shape: _arcpy.Polygon) -> str:
"""
Given a polygon instance (as read from SearchCursor 'Shape@'), get the hash value
of the wkt representation of the shape.
This is useful if we want to get a unique text value of the shape for aggregate
operations, like Dissolve, which do not support Shape.
Args:
shape: Polygon instance
Returns:
The sha256 has of shape
Examples:
>>> row = arcpy.da.SearchCursor('C:/myshape.shp', config.GeoDatabaseLayersAndTables.ERAMMPCommon.sq, ['Shape@']).next() # noqa
>>> shape_hash(row[0])
'773bd78d5d2c0bf243d0d773272be8e36c25b52fc282f79ebd96afd372e68e58'
"""
return _hashlib.sha256(shape.WKT.encode('ASCII')).hexdigest()
def polygon_from_list(lst: (list, tuple)) -> _arcpy.Polygon:
"""
Args:
lst (any): A list of points, accepts tuples as well
Returns:
arcpy.Polygon: An instance of arcpy.Polygon
Notes:
Include the closing point, so a square is defined by 5 points
Examples:
>>> ply = polygon_from_list([(0,0), (0, 1), (1, 1), (1, 0), (0, 0)]) # noqa
>>> ply.centroid
<Point (0.50006103515625, 0.50006103515625, #, #)>
With Z and M
>>> ply = polygon_from_list([(0, 0, 2, 2), (0, 1, 2, 2), (1, 1, 2, 2), (1, 0, 2, 2), (0, 0, 2, 2)]) # noqa
>>> ply.centroid
<Point (0.50006103515625, 0.50006103515625, #, #)>
"""
pts = _arcpy.Array([_arcpy.Point(*pt) for pt in lst])
poly = _arcpy.Polygon(pts)
return poly
def polygon_from_extent(ext: _arcpy.Extent):
"""Make an arcpy polygon object from an input extent object.
Largely superflous because Extent instances expose the Polygon obj as a property and vica-versa
Args:
ext (_arcpy.Extent): Instance of arcpy.Extent
Returns:
_arcpy.Polygon: A Polygon instance representation of the extent
Examples:
Circular example for illustration
>>> poly = polygon_from_extent(arcpy.da.SearchCursor('C:/my.gdb/countries', ['SHAPE@']).next()[0].extent) # noqa
>>> arcpy.CopyFeatures_management(poly, r'C:/Temp/Project_boundary.shp') # noqa
"""
array = _arcpy.Array()
array.add(ext.lowerLeft)
array.add(ext.lowerRight)
array.add(ext.upperRight)
array.add(ext.upperLeft)
array.add(ext.lowerLeft)
return _arcpy.Polygon(array, ext.spatialReference)
def polyline_from_list(lst: (list[list[(int, float)]])) -> _arcpy.Polyline:
"""
Get instance of an arcpy.Polyline from a list/tuple of points
Args:
lst (any): A list of points, accepts tuples as well
Returns:
arcpy.Polyline: An instance of arcpy.Polyline
Notes:
Order matters.
Should accept any combination of lists or tuples containing ints and/or floats
Examples:
>>> poly = polyline_from_list([(0,0), (0, 1), (1, 1), (1, 0)]) # noqa
>>> poly.centroid
<Point (0.50006103515625, 0.50006103515625, #, #)>
With Z and M
>>> ply = polyline_from_list([(0, 0, 1, 1), (0, 1, 1, 1), (1, 1, 1, 1), (1, 0, 1, 1)]) # noqa
>>> ply.centroid
<Point (0.50006103515625, 0.50006103515625, 0, 0)>
"""
pts = _arcpy.Array([_arcpy.Point(*pt) for pt in lst])
poly = _arcpy.Polyline(pts)
return poly
def polyline_to_polygon(polyline: arcpy.Polyline) -> (arcpy.Polygon, None, tuple[arcpy.Polygon]):
"""
Given a polyline, create a polygon.
Takes the last point of polyline and adds a final point at the start point of polyline
Supports multipart features, but returns them as a tuple of polygons
Args:
polyline: a polyline
Returns:
None: If polyline had no points or no polygon could be generated from the line(s)
arcpy.Polygon: Polygon instance, if polyline was a single polyline feature
tuple[arcpy.Polygon]: If the polyline was a multipart feature, i.e. breaks multipart polylines into individual polygons
"""
line: _arcpy.Array
out = []
for line in polyline:
if points_equal(line[0], line[-1]):
out += [arcpy.Polygon(line)]
else:
line.add(line[0])
out += [arcpy.Polygon(line)]
if not out:
return None
if len(out) == 1:
return out[0]
return out
def points_equal(pt1: _arcpy.Point, pt2: _arcpy.Point) -> bool:
"""
Are two points equal.
pt1 == pt2 doesnt work.
Args:
pt1: arcpy.Point instance
pt2: arcpy.Point instance
Returns:
bool
"""
return pt1.X == pt2.X and pt1.Y == pt2.Y and pt1.Z == pt2.Z
def point_from_list(lst) -> _arcpy.Point:
"""
Get instance of an arcpy.Point from a list/tuple of points
Args:
lst (tuple, list): A list of points, accepts tuples as well
Returns:
arcpy.Point: An instance of arcpy.Point
Notes:
Added for completeness - creating an arcpy point is trivial
Examples:
>>> pt = point_from_list((0, 0)) # noqa
<Point (0.50006103515625, 0.50006103515625, #, #)>
With Z and M
>>> pt = point_from_list((0, 0, 0, 0)) # noqa
<Point (0.50006103515625, 0.50006103515625, 0, 0)>
"""
return _arcpy.Point(*lst)
def array_from_list(lst: (list, tuple)) -> _arcpy.Array:
"""
Get an arcpy.Array instance from a list or tuple
Args:
lst (any): A list of points, accepts tuples as well
Returns:
arcpy.Polygon: An instance of arcpy.Polygon
Notes:
Include the closing point, so a square is defined by 5 points
Examples:
>>> array_from_list([(0,0), (0, 1), (1, 1), (1, 0), (0, 0)])
"""
return _arcpy.Array([_arcpy.Point(x, y) for x, y in lst])
def list_from_polygon(poly: _arcpy.Polygon):
"""
Return arcpy polygon instance as a list of points
Args:
poly: An arcpy polygon object
Returns:
list: 2d-list, e.g. [[0,0], [0, 1], [1, 1], [1, 0], [0, 0]]
Examples:
>>> Poly = polygon_from_list([(0,0), (0, 1), (1, 1), (1, 0), (0, 0)])
>>> list_from_polygon(Poly)
[[0,0], [0, 1], [1, 1], [1, 0], [0, 0]]
"""
return [[pt.X, pt.Y] for pt in poly[0]]
def extent_zoom(extent: _arcpy.Extent, factor_perc: float = None, factor_abs: float = None) -> _arcpy.Extent:
"""
Expand or contract an extent by a percent or absolute factor
Args:
extent:
factor_perc:
factor_abs:
Returns:
"arcpy.Extent": An arcpy extent object
"""
if not (factor_abs or factor_perc):
return extent
if factor_abs and factor_perc:
raise ValueError('factor_perc and factor_abs were passed. Use one or the other.')
if factor_abs:
buff_dist = factor_abs
else:
buff_dist = ((int(abs(extent.lowerLeft.X - extent.lowerRight.X))) * (factor_perc/100))
return extent.polygon.buffer(buff_dist).extent
if __name__ == '__main__':
# quick tests here
ply = polygon_from_list([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)])
extent_zoom(ply.extent, factor_abs=0.4)