Triangulace geometrie byla kdysi "věcí".
Lze ji provést, pokud máte rozšíření 3D Analyst, ale pokud ne a potřebujete znát principy a umíte pracovat s NumPy a pythonem... tady to je.
1 vytvořte pole z geometrie (FeatureClassToNumPyArray je začátek. Už jsem dříve zveřejnil kód, jak to udělat
2 triangulujte ji (používám scipy)
3 získejte centroid každého trojúhelníku. (opět jsem zveřejnil kód pro toto.)
4 zkontrolujte, zda jsou centroidy uvnitř původního tvaru.
Abych to zkomplikoval, polygony jsou buď konvexní nebo konkávní, s dírami nebo bez nich a mohou se skládat z více částí. Takže někdy musíte svou geometrii rozložit na nejjednodušší formu a projít ji cyklicky. Konvexní tvary jsou nudné, takže je ignoruji.
Tady jsou 3 tvary, C, D, A. C je konkávní tvar bez děr. D je konvexní, ale má díru. A je konkávní a má díru.
Jejich triangulované verze.
A omezená triangulace. To znamená triangulace s trojúhelníky, které nejsou v původním vnějším obalu ani v díře.
Vypadá to, že je vše v pořádku.
Pro provedení triangulace můžete použít
from scipy.spatial import Delaunay
def triangulate_pnts(pnts):
"""Triangulujte body a vraťte trojúhelníky.
Parametry
----------
pnts : pole
Body pro tvar nebo skupinu bodů v poli.
Buď geo.shapes nebo np.ndarray.
out : pole
Pole bodů trojúhelníků.
.. poznámka::
Simplices jsou uspořádány proti směru hodinových ručiček, toto je v této
implementaci obráceno.
"""
pnts = np.unique(pnts, axis=0) # získat pouze unikátní body
avg = np.mean(pnts, axis=0)
p = pnts - avg
tri = Delaunay(p)
simps = tri.simplices
# -- držák indexů, naplňte indexy, zopakujte první a posuňte doleva
# přeložte body zpět
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_: # přeuspořádat po směru hodinových ručiček
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, používáno funkcemi `areas` a `centroids`.
Negativní plochy jsou díry. Toto je úmyslně obráceno oproti
vzorci `shoelace`.
"""
a = _base_(a)
x0, y1 = (a.T)[:, 1:] # křížové nastavení následovně
x1, y0 = (a.T)[:, :-1]
e0 = np.einsum('...i,...i->...i', x0, y0) # upraveno 2024-03-28
e1 = np.einsum('...i,...i->...i', x1, y1)
return np.sum((e0 - e1) * 0.5)
def _area_centroid_(a):
r"""Vypočítat plochu a centroid pro jednoprvkový polygon `a`.
Používá se také k výpočtu plochy a centroidu částí Geo pole.
Poznámky
-----
Pro vícedílné tvary použijte tento zápis:
>>> # obdélník s dírou
>>> 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])
Tady je 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]])
Centroidy a plochy i triangulace lze určit pomocí výše uvedeného.
Testování které trojúhelníky patří do původní geometrie vyžaduje použití hledání bodů v polygonu. Udělal jsem to už ve své předchozí inkarnaci.
'Point in Polygon ... Geometry Mysteries - Esri Community' (Bod v polygonu ... Geometrické záhady - Esri komunita)
'Používá přístup založený na winding number.' (Používá přístup založený na čísle obtáčení.)
'Takže pokud potřebujete omezenou Delaunayovu triangulaci svých geometrických objektů... zkuste to.'
'Poznámka:'
'Pokud chcete vidět více pythonu, NumPy a ArcPy věcí podívejte se na můj github.'
'Dan Patterson na githubu'
' '
