Simplify your poly* geometry ... part 2

495
0
03-11-2023 08:27 PM
Labels (1)
DanPatterson
MVP Esteemed Contributor
3 0 495

This is a short missive that emerged as part of the Fun with Clipping series

Fun with clipping ... part 1 

Let's keep it simple.  You have a poly* feature like the one on the left and you want to get rid of all the extra points that are obviously on the straight line segments and keep the main points needed to maintain the geometry.

simplify.png

What could be easier.  Grab your points into numpy.  Apply a little numpy-magic and you have the image on the right.

 

def simplify(arr, tol=1e-6):
    """Remove redundant points on a poly perimeter."""
    x1, y1 = arr[:-2].T
    x2, y2 = arr[1:-1].T
    x3, y3 = arr[2:].T
    result = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)
    whr = np.nonzero(np.abs(result) >= tol)[0]
    bits = arr[1:-1][whr]
    keep = np.concatenate((arr[0, None], bits, arr[-1, None]), axis=0)
    return keep

 

The x1,y1 points represent the start coordinates of the sliding test with x3,y3 the end and obviously x2,y2 the test point in the middle. 

So as line 6 slides along the line it is doing a redundancy check as to whether the middle point is actually needed... that is are all three points on the same sliding segment.

Line 7 and 8 do the big where-oh-where-is test.  To avoid floating point issues (which there shouldn't be, but to be safe) a small inclusion tolerance.  The absolute differences are compared to the tolerance (or 0.0) (line 7) which results in a boolean array ( True - False or 1 - 0 ).

NumPy's nonzero allows you to get the indices of where the test yielded "True/1".

You then need to slice the original array (arr) skipping the first and last point (arr[1:-1) because of the sliding thing.  Piece the array together using "concatenate" because they are the keepers. 

The bits in the concatenate are for homework, but your clue is that the first and last points will have and array dimension of 1 while the middle points have a dimension of 2.  Ok... some hints.

 

arr[0]  
Geo([  3.50,   3.00])

arr[0].ndim, arr[-1].ndim
(1, 1)

arr[1:-1].ndim
2

# -- Now add a new axis, np.newaxis or None does the trick
arr[0, None]
Geo([[  3.50,   3.00]])  # -- the square brackets now number 2
arr[0, None].ndim
2

 

Arrays have to be "broadcastable" in order to perform many operations on them.  Just don't read up on broadcasting quite yet in your array journey.

Code for converting featureclasses or arcpy or json shapes is in my github site...

NumPy and Geometry 

 

 

About the Author
Retired Geomatics Instructor (also DanPatterson_Retired). Currently working on geometry projects (various) as they relate to GIS and spatial analysis. I use NumPy, python and kin and interface with ArcGIS Pro.
Labels