Polygon geometry scenario
A bunch of triangles, some with holes. The green numbers denote the polygon outer ring ID values. The red numbers are the sequential point ID values for the outer and inner rings.
It should be visually obvious that there are locations that have ID values :
- only appear once (segment connections),
- some that have two numbers (polygon start and end),
- and others that have 3, 4 or even 5 associated IDs (where polygons meet with various configurations).
From arcpy to numpy
When you use FeatureClassToNumPyArray to output the geometry to NumPy array format, you have several options. My preference is to output the object id and values for the X and Y coordinates as a separate output. If I need to work with attributes (text or other numeric values), I produce a separate output, which can readily be joined back to any reconstructed output.
The basic syntax is:
import numpy as np
from numpy.lib.recfunctions import structured_to_unstructured as stu # needed
import arcpy
from arcpy.da import FeatureClassToNumPyArray
in_fc = r"C:\path_to_gdb\the_geodatabase.gdb\featureclass_name"
flds = ['OID@', 'SHAPE@X', 'SHAPE@Y'] # -- the fields to use
out = FeatureClassToNumPyArray(in_fc, flds, explode_to_points=True)
Of course, you are working with some locally stored data to make things cleaner, so modify line 3 above for your use case.
Line 4 specifies the field tokens to use to bring in the OBJECTID field and values for the polygon vertices X and Y values.
The rules I use:
- I validate the input featureclass using OGC geometry validation to ensure that polygons with holes are ordered such that the holes (inner rings) follow their parent outer ring. I just can't see a use case where the geometry for a hole comes before its outer ring. It saves having to assign hole rings to their outer ring with point-in-polygon checks.
- Projected data and a projection that is at least Conformal (UTM, MTM) or in some cases equal area. Definitely not those webbie projections which are useless for determining geometry attributes.
- Outer rings are clockwise oriented, holes (inner rings) are counter clockwise.
- The start and end point for a polygon has duplicate X, Y values but a different ID number.
- plus all the other stuff that ESRI specifies.
Read Check Geometry and Repair Geometry for pertinent information.
Basic array output
Now back to the array generation.
# ---- carrying on from above
ids = out['OID@'].astype('int') # -- get a simple numpy ndarray of values
xys = stu(out[['SHAPE@X', 'SHAPE@Y']]) # -- convert structured array to ndarray
LL = np.min(xys, axis=0)
arr = xys - LL
#
ids
[1 1 1 1 1 1 1 1 1 ... 5 5 5 5]
arr
[[ 10.0, 0.0]
[ 20.0, 20.0]
[ 30.0, 0.0]
[ 10.0, 0.0]
[ 17.5, 5.0]
[ 22.0, 5.0]
[ 22.0, 10.0]
[ 17.5, 10.0]
[ 17.5, 5.0]
# ... snip
[ 40.0, 0.0]
[ 35.0, 10.0]
[ 40.0, 20.0]
[ 40.0, 0.0]]
Line 2 splits off the objectid values from the derived structured array as simple integers.
Line 3 performs a similar operation by pulling in the X, Y coordinates as a double numpy ndarray (stuis a numpy import and rename in the first code block)
Line 4 shifts the coordinates to a less onerous planar space. I subtract the minimum X and Y values from all the coordinate values. In my example, coordinate (0.0, 0.0) is actually (300,000.00, 5,000,000.00). The XY plane has been shifted making things easier to read and calculations cleaner. You shift any calculations to coordinate back when you are done, so there is no harm.
Lines 7, 8 are the id values for the first shape (1) which has 9 points forming the outer ring and its hole and the 4 points associated with the last shape (5). Lines 10-24 are the coordinates for the points. Line 11 and 14 are the start-end coordinates for the outer ring (triangle) and 15-19 are the points for the square hole.
What can be done now
At this point you have a bunch of points representing polygon geometry. You could
- shift them (simply add/subtract to/from the X, Y values)
- densify the number of points along the rings or simplify them (there are numpy functions for this
- etcetera
But if you have to work with individual polygons, then you have to split them. This is where there are several options to isolate the polygons themselves and/or their inner and outer rings.
def array_to_rings(arr, from_to_ids=False):
"""Partition an array into outer and inner rings.
Parameters
----------
arr : array_like
An Nx2 array/list/geo array of x, y coordinates, which may contain one
or more outer and inner (holes) rings.
Returns
-------
A list of arrays, possibly nested if there are inner rings/holes.
"""
# --
def _area_(a):
"""Modified `_bit_area_`. Negative areas are holes. Reversed from
the `shoelace` formula.
"""
x0, y1 = (a.T)[:, 1:] # cross set up as follows
x1, y0 = (a.T)[:, :-1]
e0 = np.einsum('...i,...i->...i', x0, y0) # 2024-03-28 modified
e1 = np.einsum('...i,...i->...i', x1, y1)
return np.sum((e0 - e1) * 0.5)
def _is_ccw_(a):
"""Counterclockwise check. Reversed from shoelace formula."""
return 0 if _area_(a) > 0. else 1
# --
#
msg = "\nAn Nx2 np.ndarray is required."
if not isinstance(arr, np.ndarray):
print(msg)
return None
if hasattr(arr, "IFT"):
arr = np.copy(arr, subok=False)
if len(arr.shape) != 2:
print(msg)
return None
#
N = len(arr)
fst = 0
out = []
n = 0
c_o = 0 # -- outer ring counter
outer_ = []
inner_ = []
frto_ = [] # -- from-to ids to split the input array `arr`
while n < N:
ids = np.nonzero((arr[fst] == arr[fst:]).all(-1))[0]
if len(ids) >= 2:
pair = ids[:2]
else:
print("\nNo matching point, input array is defective.")
return None
fst, lst = pair + n
sub = arr[fst:lst + 1]
n = n + len(sub) # size of sub array
res = _is_ccw_(sub)
if res != 1:
outer_.append(sub)
out.append([sub])
c_o += 1
else:
inner_.append(sub)
if len(out) > 0: # append hole rings to an existing outer ring
out[-1].append(sub)
else: # just in case a hole is encountered before an outer ring
out.append([sub])
frto_.append([c_o, fst, lst + 1])
fst = lst + 1
# -- return a list of arrays or a list of lists of arrays in some form
out = [i[0] if len(i) == 1 else i for i in out]
if from_to_ids:
return out, frto_
return out
And the output:
array_to_rings(A5)
([[array([[ 10.000, 0.000],
[ 20.000, 20.000],
[ 30.000, 0.000],
[ 10.000, 0.000]]),
array([[ 17.500, 5.000],
[ 22.000, 5.000],
[ 22.000, 10.000],
[ 17.500, 10.000],
[ 17.500, 5.000]])],
[... snip],
[... snip],
array([[ 30.000, 0.000],
[ 35.000, 10.000],
[ 40.000, 0.000],
[ 30.000, 0.000]]),
array([[ 40.000, 0.000],
[ 35.000, 10.000],
[ 40.000, 20.000],
[ 40.000, 0.000]])]
The output above shows the first shape is a list of two arrays, the outer ring and the hole. The last two arrays are for shapes 4 and 5 which are simple triangles.
We could at this point produce an output that removes the holes.
out = []
for i in arr:
if isinstance(i, list):
out.append(i[0])
else:
out.append(i)
out
[array([[ 10.000, 0.000],
[ 20.000, 20.000],
[ 30.000, 0.000],
[ 10.000, 0.000]]),
array([[ 30.000, 0.000],
[ 20.000, 20.000],
[ 40.000, 20.000],
[ 30.000, 0.000]]),
array([[ 10.000, 0.000],
[ 0.000, 20.000],
[ 20.000, 20.000],
[ 10.000, 0.000]]),
array([[ 30.000, 0.000],
[ 35.000, 10.000],
[ 40.000, 0.000],
[ 30.000, 0.000]]),
array([[ 40.000, 0.000],
[ 35.000, 10.000],
[ 40.000, 20.000],
[ 40.000, 0.000]])]
You can then convert the arrays to arcpy polygons.
def array_poly(arrs, p_type=None, sr=None):
"""Assemble poly features from Geo or ndarrays.
Used by `geometry_fc` or it can be used separately.
Parameters
----------
arr : arrays
Points array(s).
p_type : text
POLYGON or POLYLINE
sr : spatial reference
Spatial reference object, name or id.
Notes
-----
Polyline or polygon features can be created from the array data. The
features can be multipart with or without interior rings.
Outer rings are ordered clockwise, inner rings (holes) are ordered
counterclockwise.
"""
def _arr_poly_(arr, SR, as_type):
"""Slice the array where nan values appear, splitting them off."""
aa = [Point(*pairs) for pairs in arr]
if as_type.upper() == 'POLYGON':
poly = Polygon(Array(aa), SR)
elif as_type.upper() == 'POLYLINE':
poly = Polyline(Array(aa), SR)
return poly
# --
polys = []
if isinstance(arrs, np.ndarray):
if len(arrs.shape) == 2:
arrs = [arrs]
elif arrs.dtype == 'O':
arrs = arrs
for a in arrs:
p = _arr_poly_(a, sr, p_type)
polys.append(p)
return polys
And then to a featureclass using the above.
def geometry_fc(a, p_type=None, gdb=None, fname=None, sr=None):
"""Form poly features from the list of arrays created by `fc_geometry`.
Parameters
----------
a : array or list of arrays
Some can be object arrays, normally created by ``pnts_arr``
p_type : string
Uppercase geometry type eg POLYGON.
gdb : text
Geodatabase path and name.
fname : text
Featureclass name.
sr : spatial reference
name or object
Returns
-------
Singlepart and/or multipart featureclasses.
Notes
-----
The work is done by ``array_poly``.
"""
if p_type is None:
p_type = "POLYGON"
out = array_poly(a, p_type.upper(), sr=sr, IFT=IFT) # call array_poly
name = gdb + "/" + fname
wkspace = arcpy.env.workspace = r'memory' # legacy is in_memory
CreateFeatureclass(wkspace, fname, p_type, spatial_reference=sr)
AddField(fname, 'ID_arr', 'LONG')
with InsertCursor(fname, ['SHAPE@', 'ID_arr']) as cur:
for row in out:
cur.insertRow(row)
CopyFeatures(fname, name)
return None
Which is equivalent to the Eliminate Polygon Part tool which strangely requires an Advanced License for Pro. Now why would you need an advanced license to remove/fill a hole.
By changing the output type in the above to POLYLINE, you change polygons to polylines just like Polygon to Line which also requires an Advanced license.
This is just a simple example of how you can use FeatureClassToNumPyArray to work with geometry.
More tools and code available at NumPy Geometry and Tools for ArcGIS Pro.
