Densify by distance
Densification
Densification of polygon boundaries or polyline segments, is a common task.
In a previous incarnation, I wrote about densifying geometry based on a "factor", that is, double, triple... the number of points along a segment.
Planar densification can be carried out in Editing Tools, with several options. It is not available at the Basic license level for some reason.
Densify (Editing)—ArcGIS Pro | Documentation
This missive concentrates on densification based on a planar distance step along the line/perimeter. It is the most common usage.
It will included a couple of my toolsets (listed below).
The code is illustrative of the power of numpy in geoprocessing.
def _pnts_on_line_(a, spacing=1, is_percent=False): # densify by distance
"""Add points, at a fixed spacing, to an array representing a line.
Parameters
----------
a : array
A sequence of `points`, x,y pairs, representing the bounds of a polygon
or polyline object.
spacing : number
Spacing between the points to be added to the line.
is_percent : boolean
Express the densification as a percent of the total length.
"""
a = _get_base_(a)
N = len(a) - 1 # segments
dxdy = a[1:, :] - a[:-1, :] # coordinate differences
leng = np.sqrt(np.einsum('ij,ij->i', dxdy, dxdy)) # segment lengths
if is_percent: # as percentage
spacing = abs(spacing)
spacing = min(spacing / 100, 1.)
steps = (sum(leng) * spacing) / leng # step distance
else:
steps = leng / spacing # step distance
deltas = dxdy / (steps.reshape(-1, 1)) # coordinate steps
pnts = np.empty((N,), dtype='O') # construct an `O` array
for i in range(N): # cycle through the segments and make
num = np.arange(steps[i]) # the new points
pnts[i] = np.array((num, num)).T * deltas[i] + a[i]
a0 = a[-1].reshape(1, -1) # add the final point and concatenate
return np.concatenate((*pnts, a0), axis=0)
I think my favorite 2 lines are ...
dxdy = a[1:, :] - a[:-1, :] # coordinate differences
leng = np.sqrt(np.einsum('ij,ij->i', dxdy, dxdy)) # segment lengths
Einstein was a smart dude.
His notation syntax is implemented in many languages and causes brain squint until you get used to it (fodder for another blog).
In short, dxdy represents the sequential differences in an array's coordinates (read polygon/polyline coordinates).
None of this reading the geometry and subtracting each x and y coordinate in turn. One line and you are done.
Now a[1:, :] means from the second coordinate pair onward and a[:-1, :] means from the first coordinate upto but not including the last (indexing is zero-based).
Head swimming? You have been using Pandas or shapely too long.
Here is what a simple polygon looks like represented like an array (even an arcpy.Array )
a # ---- an array of x, y coordinates... first and last are the same
array([[ 0.00, 0.00],
[ 2.00, 8.00],
[ 8.00, 10.00],
[ 10.00, 10.00],
[ 10.00, 8.00],
[ 9.00, 1.00],
[ 0.00, 0.00]])
dxdy # ---- sequential, delta x and delta y between the pairs
array([[ 2.00, 8.00],
[ 6.00, 2.00],
[ 2.00, 0.00],
[ 0.00, -2.00],
[ -1.00, -7.00],
[ -9.00, -1.00]])
leng # ---- the inter-point pair distance
array([ 8.25, 6.32, 2.00, 2.00, 7.07, 9.06])
Now that the magic is done, the rest is just stuff to determine whether you want to densify based on an absolute distance or as a percentage of the total distance.
The for loop section can be simplified, however, I kept it as such since it shows how each segment gets densified and the results collected and then assembled into the final result,
So what does it look like? The simple polygon shown above and the densified version are as follows.
Plots done in Spyder using Matplotlib.
Finally, I will add this capability to my Polyline/Polygon tools and Free Tools tools shortly.
Free_Tools in Tools_for_ArcGIS_Pro (github.com)
PolygonLine Tools in Tools_for_ArcGIS_Pro (github.com)
If you are interested in geometry, computational geometry and numpy and its interface with arcpy and arcgis pro, then my Numpy Geometry link may be of interest. It is the foundation for all the toolset in my site.
numpy_geometry: A numpy geometry class and functions that work with arcpy and ESRI featureclasses
That's all for now. Give me a buzz if you need something you would like added... especially if it requires and Advance license 😉.
--------------------------
Previous work
Densification.... sometimes being dense is a good ... - GeoNet, The Esri Community