Scenario met polygoongeometrie
Een heleboel driehoeken, sommige met gaten. De groene getallen geven de ID-waarden van de buitenring van de polygonen aan. De rode getallen zijn de opeenvolgende punt-ID-waarden voor de buiten- en binnenringen.
Visueel is duidelijk te zien dat er locaties zijn met ID-waarden die:
- slechts één keer voorkomen (segmentverbindingen),
- twee getallen hebben (begin- en eindpunt van een polygoon),
- en andere met 3, 4 of zelfs 5 bijbehorende ID's (waar polygonen elkaar in verschillende configuraties raken).
Van arcpy naar numpy
Wanneer je FeatureClassToNumPyArray gebruikt om de geometrie naar de indeling van een NumPy-array uit te voeren, heb je verschillende opties. Ik geef er de voorkeur aan de object-ID en de waarden voor de X- en Y-coördinaten als afzonderlijke uitvoer te produceren. Als ik met attributen moet werken (tekst of andere numerieke waarden), maak ik daarvoor een aparte uitvoer, die eenvoudig weer aan elke gereconstrueerde uitvoer kan worden gekoppeld.
De basissyntaxis is:
import numpy as np
from numpy.lib.recfunctions import structured_to_unstructured as stu # nodig
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'] # -- de te gebruiken velden
out = FeatureClassToNumPyArray(in_fc, flds, explode_to_points=True)
Je werkt natuurlijk met lokaal opgeslagen gegevens om het overzichtelijk te houden, dus pas regel 3 hierboven aan voor jouw situatie.
Regel 4 geeft de veldtokens op waarmee je het OBJECTID-veld en de waarden voor de X- en Y-coördinaten van de polygoonvertices ophaalt.
De regels die ik hanteer:
- Ik valideer de invoerfeatureclass met OGC-geometrievalidering om ervoor te zorgen dat polygonen met gaten zo zijn geordend dat de gaten (binnenringen) na hun bijbehorende buitenring komen. Ik kan geen situatie bedenken waarin de geometrie van een gat vóór de buitenring komt. Zo hoef je gaten niet met point-in-polygon-controles aan hun buitenring toe te wijzen.
- Geprojecteerde gegevens en een projectie die ten minste conform is (UTM, MTM) of in sommige gevallen een equivalente projectie. Absoluut geen van die webprojecties, die onbruikbaar zijn voor het bepalen van geometrieattributen.
- Buitenringen zijn met de klok mee georiënteerd; gaten (binnenringen) tegen de klok in.
- Het begin- en eindpunt van een polygoon hebben dezelfde X- en Y-waarden, maar een verschillend ID-nummer.
- plus alle andere zaken die ESRI voorschrijft.
Lees Geometrie controleren en Geometrie herstellen voor relevante informatie.
Basisuitvoer van de array
Nu terug naar het genereren van de array.
# ---- ga verder vanaf hierboven
ids = out['OID@'].astype('int') # -- haal een eenvoudige numpy-ndarray met waarden op
xys = stu(out[['SHAPE@X', 'SHAPE@Y']]) # -- zet de gestructureerde array om naar een 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]
# ... ingekort
[ 40.0, 0.0]
[ 35.0, 10.0]
[ 40.0, 20.0]
[ 40.0, 0.0]]
Regel 2 haalt de object-ID-waarden uit de afgeleide gestructureerde array en zet ze om in gewone gehele getallen.
Regel 3 voert een vergelijkbare bewerking uit door de X- en Y-coördinaten op te halen als een dubbele numpy-ndarray (stuis een numpy-import met een andere naam in het eerste codeblok)
Regel 4 verschuift de coördinaten naar een overzichtelijker vlak. Ik trek de minimale X- en Y-waarden af van alle coördinaatwaarden. In mijn voorbeeld is coördinaat (0.0, 0.0) in werkelijkheid (300,000.00, 5,000,000.00). Het XY-vlak is verschoven, waardoor alles makkelijker af te lezen is en berekeningen eenvoudiger worden. Als je klaar bent, kun je de berekeningen terugschuiven naar de oorspronkelijke coördinaten; dat kan dus geen kwaad.
Regels 7 en 8 bevatten de ID-waarden voor de eerste vorm (1), met 9 punten die de buitenring en het gat vormen, en voor de 4 punten van de laatste vorm (5). Regels 10-24 bevatten de coördinaten van de punten. Regels 11 en 14 bevatten de begin- en eindcoördinaten van de buitenring (driehoek); regels 15-19 bevatten de punten van het vierkante gat.
Wat je nu kunt doen
Op dit punt heb je een heleboel punten die de polygoongeometrie voorstellen. Je kunt ze
- verschuiven (eenvoudigweg waarden optellen bij of aftrekken van de X- en Y-waarden)
- het aantal punten langs de ringen vergroten of ze vereenvoudigen (hiervoor zijn numpy-functies
- enzovoort
Maar als je met afzonderlijke polygonen moet werken, moet je ze splitsen. Er zijn verschillende opties om de polygonen zelf en/of hun binnen- en buitenringen afzonderlijk te maken.
def array_to_rings(arr, from_to_ids=False):
"""Verdeel een array in buiten- en binnenringen.
Parameters
----------
arr : array_like
Een Nx2-array/lijst/geo-array met x-, y-coördinaten, die een of
meer buiten- en binnenringen (gaten) kan bevatten.
Returns
-------
Een lijst met arrays, eventueel genest als er binnenringen/gaten zijn.
"""
# --
def _area_(a):
"""Aangepaste versie van `_bit_area_`. Negatieve oppervlakken zijn gaten. Omgekeerd ten opzichte van
de `shoelace`-formule.
"""
x0, y1 = (a.T)[:, 1:] # kruisproduct als volgt opgesteld
x1, y0 = (a.T)[:, :-1]
e0 = np.einsum('...i,...i->...i', x0, y0) # aangepast op 2024-03-28
e1 = np.einsum('...i,...i->...i', x1, y1)
return np.sum((e0 - e1) * 0.5)
def _is_ccw_(a):
"""Controle op tegen-de-klok-inrichting. Omgekeerd ten opzichte van de shoelace-formule."""
return 0 if _area_(a) > 0. else 1
# --
#
msg = "\nEen Nx2 np.ndarray is vereist."
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 # -- teller voor buitenringen
outer_ = []
inner_ = []
frto_ = [] # -- van-tot-ID's om de invoerarray `arr` te splitsen
while n < N:
ids = np.nonzero((arr[fst] == arr[fst:]).all(-1))[0]
if len(ids) >= 2:
pair = ids[:2]
else:
print("\nGeen overeenkomend punt; de invoerarray is defect.")
return None
fst, lst = pair + n
sub = arr[fst:lst + 1]
n = n + len(sub) # grootte van de subarray
res = _is_ccw_(sub)
if res != 1:
outer_.append(sub)
out.append([sub])
c_o += 1
else:
inner_.append(sub)
if len(out) > 0: # voeg gatenringen toe aan een bestaande buitenring
out[-1].append(sub)
else: # voor het geval er een gat wordt aangetroffen vóór een buitenring
out.append([sub])
frto_.append([c_o, fst, lst + 1])
fst = lst + 1
# -- geef een lijst met arrays of een lijst met lijsten van arrays in een bepaalde vorm terug
out = [i[0] if len(i) == 1 else i for i in out]
if from_to_ids:
return out, frto_
return out
En de uitvoer:
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]])],
[... ingekort],
[... ingekort],
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]])]
De bovenstaande uitvoer laat zien dat de eerste vorm een lijst met twee arrays is: de buitenring en het gat. De laatste twee arrays zijn voor vormen 4 en 5, eenvoudige driehoeken.
We kunnen nu een uitvoer maken waarin de gaten zijn verwijderd.
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]])]
Vervolgens kun je de arrays omzetten naar arcpy-polygonen.
def array_poly(arrs, p_type=None, sr=None):
"""Stel polygoonobjecten samen uit Geo- of ndarrays.
Gebruikt door `geometry_fc` of kan afzonderlijk worden gebruikt.
Parameters
----------
arr : arrays
Puntenarray(s).
p_type : tekst
POLYGON of POLYLINE
sr : ruimtelijke referentie
Ruimtelijk-referentieobject, naam of ID.
Opmerkingen
---------
Polyline- of polygoonobjecten kunnen worden gemaakt op basis van de arraygegevens. De
objecten kunnen multipart zijn, met of zonder binnenringen.
Buitenringen zijn met de klok mee geordend; binnenringen (gaten) zijn
tegen de klok in geordend.
"""
def _arr_poly_(arr, SR, as_type):
"""Knip de array op de plaatsen waar NaN-waarden voorkomen en splits de delen af."""
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
Vervolgens zet je deze met behulp van het bovenstaande om naar een featureclass.
def geometry_fc(a, p_type=None, gdb=None, fname=None, sr=None):
"""Maak polygoonobjecten van de lijst met arrays die door `fc_geometry` is gemaakt.
Parameters
----------
a : array of lijst met arrays
Sommige kunnen objectarrays zijn, die doorgaans door ``pnts_arr`` worden gemaakt.
p_type : string
Geometrie-type in hoofdletters, bijvoorbeeld POLYGON.
gdb : tekst
Pad en naam van de geodatabase.
fname : tekst
Naam van de featureclass.
sr : ruimtelijke referentie
naam of object
Retourneert
-------
Enkelvoudige en/of meerdelige featureclasses.
Opmerkingen
-----
Het werk wordt uitgevoerd door ``array_poly``.
"""
if p_type is None:
p_type = "POLYGON"
out = array_poly(a, p_type.upper(), sr=sr, IFT=IFT) # array_poly aanroepen
name = gdb + "/" + fname
wkspace = arcpy.env.workspace = r'memory' # legacy is 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
Wat gelijkwaardig is aan de tool Polygon Part verwijderen waarvoor vreemd genoeg een Advanced-licentie voor Pro vereist is. Waarom zou je een Advanced-licentie nodig hebben om een gat te verwijderen of op te vullen?
Door het uitvoertype hierboven te wijzigen in POLYLINE, zet je polygonen om in polylijnen, net als met Polygon naar lijn waarvoor ook een Advanced-licentie vereist is.
Dit is slechts een eenvoudig voorbeeld van hoe je FeatureClassToNumPyArray kunt gebruiken om met geometrie te werken.
Meer tools en code zijn beschikbaar op NumPy Geometry en Tools voor ArcGIS Pro.
