Scénář polygonové geometrie
Hromada trojúhelníků, některé s dírami. Zelená čísla označují ID hodnoty vnějšího prstence polygonu. Červená čísla jsou pořadová ID bodů vnějších a vnitřních prstenců.
Mělo by být na první pohled zřejmé, že existují místa s ID hodnotami:
- které se vyskytují pouze jednou (napojení segmentů),
- některé mají dvě čísla (začátek a konec polygonu),
- a další mají 3, 4, nebo dokonce 5 přiřazených ID (tam, kde se polygony setkávají v různých konfiguracích).
Od arcpy k numpy
Když pomocí FeatureClassToNumPyArray převedete geometrii do formátu pole NumPy, máte několik možností. Preferuji vypsat ID objektu a hodnoty souřadnic X a Y jako samostatný výstup. Pokud potřebuji pracovat s atributy (textovými nebo jinými číselnými hodnotami), vytvořím samostatný výstup, který lze snadno připojit zpět k libovolnému rekonstruovanému výstupu.
Základní syntaxe je:
import numpy as np
from numpy.lib.recfunctions import structured_to_unstructured as stu # potřeba
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'] # -- pole, která se mají použít
out = FeatureClassToNumPyArray(in_fc, flds, explode_to_points=True)
Samozřejmě pracujete s lokálně uloženými daty, aby byl postup přehlednější, takže si výše uvedený řádek 3 upravte podle svého případu.
Řádek 4 určuje tokeny polí, které se použijí k načtení pole OBJECTID a hodnot X a Y vrcholů polygonu.
Řídím se těmito pravidly:
- Vstupní featureclass ověřuji pomocí validace geometrie OGC, abych se ujistil, že polygony s dírami jsou seřazeny tak, že díry (vnitřní prstence) následují po svém nadřazeném vnějším prstenci. Nedokážu si představit případ, kdy by geometrie díry předcházela jejímu vnějšímu prstenci. Odpadá tak nutnost přiřazovat prstence děr k vnějším prstencům pomocí testů bodu v polygonu.
- Projekce dat a zobrazení, které je alespoň konformní (UTM, MTM), případně v některých případech ekvivalentní ploše. Rozhodně nepoužívejte ty webové projekce, které jsou k určování geometrických atributů k ničemu.
- Vnější prstence jsou orientovány po směru hodinových ručiček, díry (vnitřní prstence) proti směru hodinových ručiček.
- Počáteční a koncový bod polygonu mají stejné souřadnice X, Y, ale odlišné ID.
- a vše ostatní, co specifikuje ESRI.
Přečtěte si Kontrola geometrie a Oprava geometrie , kde najdete důležité informace.
Základní výstup pole
A teď zpět k vytváření pole.
# ---- pokračování výše uvedeného kódu
ids = out['OID@'].astype('int') # -- získání jednoduchého numpy ndarray hodnot
xys = stu(out[['SHAPE@X', 'SHAPE@Y']]) # -- převod strukturovaného pole na 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]
# ... zkráceno
[ 40.0, 0.0]
[ 35.0, 10.0]
[ 40.0, 20.0]
[ 40.0, 0.0]]
Řádek 2 odděluje hodnoty objectid od odvozeného strukturovaného pole jako jednoduchá celá čísla.
Řádek 3 provádí podobnou operaci, kdy načítá souřadnice X, Y jako dvojité numpy ndarray (stuje import z numpy a jeho přejmenování v prvním bloku kódu)
Řádek 4 posouvá souřadnice do přehlednější roviny. Od všech souřadnic odečítám minimální hodnoty X a Y. V mém příkladu je souřadnice (0.0, 0.0) ve skutečnosti (300,000.00, 5,000,000.00). Rovina XY se posunula, takže se v ní lépe orientuje a výpočty jsou přehlednější. Po dokončení výpočtů posuňte souřadnice zpět, takže tím nic nezkazíte.
Řádky 7 a 8 obsahují hodnoty ID prvního útvaru (1), který má 9 bodů tvořících vnější prstenec a jeho díru, a 4 body související s posledním útvarem (5). Řádky 10–24 obsahují souřadnice bodů. Řádky 11 a 14 obsahují počáteční a koncové souřadnice vnějšího prstence (trojúhelníku) a řádky 15–19 obsahují body čtvercové díry.
Co lze nyní dělat
V tuto chvíli máte hromadu bodů představujících polygonovou geometrii. Můžete
- je posunout (jednoduše přičíst nebo odečíst hodnoty X, Y),
- zvýšit počet bodů podél prstenců nebo je zjednodušit (existují k tomu funkce numpy,
- a tak dále
Pokud ale potřebujete pracovat s jednotlivými polygony, musíte je rozdělit. Zde máte několik možností, jak oddělit samotné polygony a/nebo jejich vnitřní a vnější prstence.
def array_to_rings(arr, from_to_ids=False):
"""Rozdělí pole na vnější a vnitřní prstence.
Parametry
----------
arr : array_like
Pole/seznam/geografické pole souřadnic x, y o rozměrech Nx2, které může obsahovat jeden
nebo více vnějších a vnitřních prstenců (děr).
Vrací
-------
Seznam polí, případně vnořených, pokud obsahují vnitřní prstence/díry.
"""
# --
def _area_(a):
"""Upravená verze `_bit_area_`. Záporné plochy představují díry. Vzorec je
obrácený oproti vzorci `shoelace`.
"""
x0, y1 = (a.T)[:, 1:] # součin se sestaví 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 _is_ccw_(a):
"""Kontrola orientace proti směru hodinových ručiček. Vzorec je obrácený oproti vzorci shoelace."""
return 0 if _area_(a) > 0. else 1
# --
#
msg = "\nJe vyžadováno np.ndarray o rozměrech Nx2."
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 # -- počítadlo vnějších prstenců
outer_ = []
inner_ = []
frto_ = [] # -- ID od-do pro rozdělení vstupního pole `arr`
while n < N:
ids = np.nonzero((arr[fst] == arr[fst:]).all(-1))[0]
if len(ids) >= 2:
pair = ids[:2]
else:
print("\nNenalezen odpovídající bod, vstupní pole je vadné.")
return None
fst, lst = pair + n
sub = arr[fst:lst + 1]
n = n + len(sub) # velikost podpole
res = _is_ccw_(sub)
if res != 1:
outer_.append(sub)
out.append([sub])
c_o += 1
else:
inner_.append(sub)
if len(out) > 0: # přidání prstenců děr k existujícímu vnějšímu prstenci
out[-1].append(sub)
else: # pro případ, že se na díru narazí před vnějším prstencem
out.append([sub])
frto_.append([c_o, fst, lst + 1])
fst = lst + 1
# -- vrácení seznamu polí nebo seznamu seznamů polí
out = [i[0] if len(i) == 1 else i for i in out]
if from_to_ids:
return out, frto_
return out
A výstup:
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]])],
[... zkráceno],
[... zkráceno],
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]])]
Výše uvedený výstup ukazuje, že první útvar je seznamem dvou polí: vnějšího prstence a díry. Poslední dvě pole patří útvarům 4 a 5, což jsou jednoduché trojúhelníky.
V tuto chvíli bychom mohli vytvořit výstup bez děr.
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]])]
Pole pak můžete převést na polygony arcpy.
def array_poly(arrs, p_type=None, sr=None):
"""Sestaví polygonové prvky z Geo nebo ndarray.
Používá `geometry_fc`, případně lze použít samostatně.
Parametry
----------
arr : pole
Pole bodů.
p_type : text
POLYGON nebo POLYLINE
sr : prostorová reference
Objekt, název nebo ID prostorové reference.
Poznámky
-----
Z dat pole lze vytvořit liniové nebo polygonové prvky. Tyto prvky
mohou být vícedílné, s vnitřními prstenci nebo bez nich.
Vnější prstence jsou orientovány po směru hodinových ručiček, vnitřní prstence (díry)
proti směru hodinových ručiček.
"""
def _arr_poly_(arr, SR, as_type):
"""Rozdělí pole v místech výskytu hodnot nan."""
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
A potom pomocí výše uvedeného vytvořte featureclass.
def geometry_fc(a, p_type=None, gdb=None, fname=None, sr=None):
"""Vytvoří polygonové prvky ze seznamu polí vytvořeného funkcí `fc_geometry`.
Parametry
----------
a : pole nebo seznam polí
Některá mohou být objektová pole, obvykle vytvořená pomocí ``pnts_arr``
p_type : řetězec
Typ geometrie velkými písmeny, např. POLYGON.
gdb : text
Cesta a název geodatabáze.
fname : text
Název třídy prvků.
sr : prostorová reference
název nebo objekt
Vrací
-------
Jednodílné a/nebo vícedílné třídy prvků.
Poznámky
-----
Práci provádí ``array_poly``.
"""
if p_type is None:
p_type = "POLYGON"
out = array_poly(a, p_type.upper(), sr=sr, IFT=IFT) # volání array_poly
name = gdb + "/" + fname
wkspace = arcpy.env.workspace = r'memory' # starší varianta je 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
Což odpovídá nástroji Eliminate Polygon Part který podivně v Pro vyžaduje licenci Advanced. Proč byste potřebovali licenci Advanced k odstranění/vyplnění díry?
Změnou výstupního typu v uvedeném příkladu na POLYLINE převedete polygony na lomené čáry, stejně jako nástroj Polygon to Line, který také vyžaduje licenci Advanced.
Toto je jen jednoduchý příklad toho, jak můžete pomocí FeatureClassToNumPyArray pracovat s geometrií.
Další nástroje a kód najdete na stránce NumPy Geometry a Tools for ArcGIS Pro.
