Triangular geometría fue una "moda" en su momento.
Se puede hacer si tienes la extensión 3D Analyst, pero si no la tienes y necesitas conocer los principios y puedes usar NumPy y python... aquí vamos.
1 crea un array a partir de la geometría (FeatureClassToNumPyArray es un buen comienzo. He publicado código antes sobre cómo hacer esto
2 triangúlalo (yo uso scipy)
3 obtén el centroide de cada triángulo. (de nuevo, he publicado código para esto.)
4 verifica si los centroides están dentro de la forma original.
Ahora, para complicar esto, los polígonos son convexos o cóncavos, con o sin agujeros y pueden consistir en más de una parte. Así que a veces, necesitas simplificar tu geometría a su forma más simple y recorrerla. Las formas convexas son aburridas, así que las ignoraré.
Aquí hay 3 formas, C, D, A. C es una forma cóncava sin agujeros. D es convexa pero tiene un agujero. A es cóncava y tiene un agujero.
Sus versiones trianguladas.
Y la triangulación restringida. Es decir, la triangulación con los triángulos que no están en el casco exterior original ni en un agujero.
Parece que todo está bien.
Para realizar la triangulación, puedes usar
from scipy.spatial import Delaunay
def triangulate_pnts(pnts):
"""Triangular los puntos y devolver los triángulos.
Parámetros
----------
pnts : array
Puntos para una forma o un grupo de puntos en formato array.
Ya sea geo.shapes o np.ndarray.
out : array
Un array de puntos de triángulos.
.. nota::
Los simples están ordenados en sentido antihorario, esto se invierte en esta
implementación.
"""
pnts = np.unique(pnts, axis=0) # obtener solo los puntos únicos
avg = np.mean(pnts, axis=0)
p = pnts - avg
tri = Delaunay(p)
simps = tri.simplices
# -- contenedor de índices, llenar con índices, repetir el primero y rotar CL
# traducir los puntos de vuelta
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_: # reordenar en sentido horario
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, usado por `areas` y `centroids`.
Las áreas negativas son agujeros. Esto está intencionalmente invertido respecto a
la fórmula `shoelace`.
"""
a = _base_(a)
x0, y1 = (a.T)[:, 1:] # configuración cruzada como sigue
x1, y0 = (a.T)[:, :-1]
e0 = np.einsum('...i,...i->...i', x0, y0) # modificado 2024-03-28
e1 = np.einsum('...i,...i->...i', x1, y1)
return np.sum((e0 - e1) * 0.5)
def _area_centroid_(a):
r"""Calcular área y centroide para un polígono de una sola parte, `a`.
Esto también se usa para calcular área y centroide para las partes de un array Geo.
Notas
-----
Para formas multipartes, usa esta sintaxis:
>>> # rectángulo con agujero
>>> 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])
Aquí está 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]])
Los centroides, áreas y triangulaciones pueden determinarse usando lo anterior.
Probar qué triángulos son parte de la geometría original implica usar una búsqueda de puntos dentro del polígono. Hice esto hace tiempo en mi encarnación anterior.
Punto en Polígono ... Misterios de Geometría - Comunidad Esri
Usa el enfoque del número de vueltas.
Así que si necesitas una triangulación Delaunay restringida para tus objetos geométricos... pruébalo.
Nota:
Si quieres ver más cosas de python, NumPy y ArcPy visita mi sitio github.
Dan Patterson en github
