Het trianguleren van geometrie was ooit een "ding".
Het kan gedaan worden als je de 3D Analyst-extensie hebt, maar als je die niet hebt en je moet de principes kennen en kunt werken met NumPy en python... hier gaan we.
1 maak een array van de geometrie (FeatureClassToNumPyArray is een begin. Ik heb eerder code gepost over hoe dit te doen
2 trianguleer het (ik gebruik scipy)
3 bepaal het zwaartepunt van elke driehoek. (ook hiervoor heb ik code gepost.)
4 controleer of de zwaartepunten binnen de oorspronkelijke vorm liggen.
Om dit te compliceren, zijn polygonen ofwel convex of concave, met of zonder gaten en ze kunnen uit meer dan één deel bestaan. Soms moet je dus je geometrie terugbrengen tot de eenvoudigste vorm en doorlopen. Convexe vormen zijn saai, dus die negeer ik.
Hier zijn 3 vormen, C, D, A. C is een concave vorm zonder gaten. D is convex maar heeft een gat. A is concave en heeft een gat.
Hun getrianguleerde versies.
En de beperkte triangulatie. Dat wil zeggen, de triangulatie met de driehoeken die niet in de oorspronkelijke buitenste omtrek of in een gat liggen.
Het lijkt erop dat alles goed is.
Om de triangulatie uit te voeren, kun je gebruiken
from scipy.spatial import Delaunay
def triangulate_pnts(pnts):
"""Trianguleer de punten en retourneer de driehoeken.
Parameters
----------
pnts : array
Punten voor een vorm of een groep punten in array-formaat.
Ofwel geo.shapes of np.ndarray.
out : array
Een array van driehoekspunten.
.. note::
De simplices zijn tegen de klok in geordend, dit wordt omgekeerd in deze
implementatie.
"""
pnts = np.unique(pnts, axis=0) # haal alleen unieke punten op
avg = np.mean(pnts, axis=0)
p = pnts - avg
tri = Delaunay(p)
simps = tri.simplices
# -- houder voor indices, vul met indices, herhaal eerste en rol CL
# vertaal de punten terug
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_: # herschik met de klok mee
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, gebruikt door `areas` en `centroids`.
Negatieve gebieden zijn gaten. Dit is opzettelijk omgekeerd ten opzichte van
de `shoelace` formule.
"""
a = _base_(a)
x0, y1 = (a.T)[:, 1:] # kruis ingesteld als volgt
x1, y0 = (a.T)[:, :-1]
e0 = np.einsum('...i,...i->...i', x0, y0) # 2024-03-28 aangepast
e1 = np.einsum('...i,...i->...i', x1, y1)
return np.sum((e0 - e1) * 0.5)
def _area_centroid_(a):
r"""Bereken gebied en zwaartepunt voor een enkelvoudige polygoon, `a`.
Dit wordt ook gebruikt om gebied en zwaartepunt te berekenen voor delen van een Geo-array.
Notities
-----
Voor multipart-vormen gebruik gewoon deze syntaxis:
>>> # rechthoek met gat
>>> 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, y1)
e1 = np.einsum('...i,...i->...i', x1, y0)
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])
Hier 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]])
De zwaartepunten en gebieden en triangulaties kunnen bepaald worden met bovenstaande methoden.
'Testen welke driehoeken deel uitmaken van de oorspronkelijke geometrie vereist het gebruik van een punten-in-polygon zoekopdracht.' Ik deed dit al lang geleden in mijn vorige incarnatie.
'Punt in Polygoon ... Geometriemysteries - Esri Community'
'Het gebruikt de winding number methode.'
'Dus als je behoefte hebt aan een beperkte Delaunay-triangulatie van je geometrie-objecten... probeer het eens.'
'Opmerking:'
'Als je meer python-, NumPy- en ArcPy-spullen wilt zien, bekijk dan mijn github-site.'
'Dan Patterson op github'
