Constrained triangulation of polygons
Triangulating geometry was a "thing" at one time.
It can be done if you have the 3D Analyst extension, but if you don't and you need to know the principles and can work NumPy and python... here goes.
1 make an array out of the geometry (FeatureClassToNumPyArray is a starter. I have posted code before on how to do this
2 triangulate it (I use scipy)
3 get the centroid of each triangle. (again, I have posted code for this.)
4 check if the centroids are within the original shape.
Now to complicate this, polygons are either convex or concave, with or without holes and they can consist of more than one part. So sometimes, you need to devolve your geometry to its simplest form and cycle through. Convex shapes are boring, so I will ignore them.
Here are 3 shapes, C, D, A. C is a concave shape with no holes. D is convex but it has a hole. A is concave and has a hole.
Their triangulated versions.
And the constrained triangulation. That is, the triangulation with the triangles not in the original outer hull or in a hole.
Looks like all is good.
To perform the triangulation, you can use
from scipy.spatial import Delaunay
def triangulate_pnts(pnts):
"""Triangulate the points and return the triangles.
Parameters
----------
pnts : array
Points for a shape or a group of points in array format.
Either geo.shapes or np.ndarray.
out : array
An array of triangle points.
.. note::
The simplices are ordered counterclockwise, this is reversed in this
implementation.
"""
pnts = np.unique(pnts, axis=0) # get the unique points only
avg = np.mean(pnts, axis=0)
p = pnts - avg
tri = Delaunay(p)
simps = tri.simplices
# -- indices holder, fill with indices, repeat first and roll CL
# translate the points back
z = np.zeros((len(simps), 4), dtype='int32')
z[:, :3] = simps
z[:, 3] = simps[:, 0]
tmp_ = p[z] + avg
new_pnts= []
for i in tmp_: # reorder clockwise
if _bit_area_(i) < 0.0: # -- 2025_10_27
i = i[::-1]
new_pnts.append(i)
return new_pnts
def _bit_area_(a):
"""Mini e_area, used by `areas` and `centroids`.
Negative areas are holes. This is intentionally reversed from
the `shoelace` formula.
"""
a = _base_(a)
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 _area_centroid_(a):
r"""Calculate area and centroid for a singlepart polygon, `a`.
This is also used to calculate area and centroid for a Geo array's parts.
Notes
-----
For multipart shapes, just use this syntax:
>>> # rectangle with hole
>>> a0 = np.array([[[0., 0.], [0., 10.], [10., 10.], [10., 0.], [0., 0.]],
[[2., 2.], [8., 2.], [8., 8.], [2., 8.], [2., 2.]]])
>>> [_area_centroid_(i) for i in a0]
>>> [(100.0, array([ 5.00, 5.00])), (-36.0, array([ 5.00, 5.00]))]
"""
a = _base_(a)
x0, y1 = (a.T)[:, 1:]
x1, y0 = (a.T)[:, :-1]
e0 = np.einsum('...i,...i->...i', x0, y0)
e1 = np.einsum('...i,...i->...i', x1, y1)
t = e1 - e0
area = np.sum((e0 - e1) * 0.5)
x_c = np.sum((x1 + x0) * t, axis=0) / (area * 6.0)
y_c = np.sum((y1 + y0) * t, axis=0) / (area * 6.0)
return area, np.asarray([-x_c, -y_c])Here is C
C
array([[ 0.00, 0.00],
[ 0.00, 10.00],
[ 10.00, 10.00],
[ 10.00, 8.00],
[ 2.00, 8.00],
[ 2.00, 2.00],
[ 10.00, 2.00],
[ 10.00, 0.00],
[ 0.00, 0.00]])The centroids and areas and triangulations can be determined using the above.
Testing which triangles are part of the original geometry entails using a points-in-polygon search. I did this way back in my previous incarnation.
Point in Polygon ... Geometry Mysteries - Esri Community
It uses the winding number approach.
So if you have a need for a constrained Delaunay triangulation of your geometry objects... give it a try.
Note:
If you want to see more python, NumPy and ArcPy stuff see my github site.
