Szenario zur Polygongeometrie
Eine Reihe von Dreiecken, einige davon mit Löchern. Die grünen Zahlen bezeichnen die ID-Werte der äußeren Polygonringe. Die roten Zahlen sind die fortlaufenden Punkt-ID-Werte für die äußeren und inneren Ringe.
Es sollte auf den ersten Blick klar sein, dass es Stellen gibt, an denen ID-Werte:
- nur einmal vorkommen (Segmentverbindungen),
- manche zwei Zahlen haben (Polygonanfang und -ende),
- und andere 3, 4 oder sogar 5 zugehörige IDs haben (an denen Polygone in verschiedenen Konfigurationen zusammentreffen).
Von arcpy zu numpy
Wenn Sie mit FeatureClassToNumPyArray die Geometrie im NumPy-Array-Format ausgeben, haben Sie mehrere Möglichkeiten. Am liebsten gebe ich die Objekt-ID sowie die Werte der X- und Y-Koordinaten separat aus. Wenn ich mit Attributen arbeiten muss (Text- oder andere Zahlenwerte), erstelle ich eine separate Ausgabe, die sich problemlos wieder mit jeder rekonstruierten Ausgabe verknüpfen lässt.
Die grundlegende Syntax lautet:
import numpy as np
from numpy.lib.recfunctions import structured_to_unstructured as stu # erforderlich
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'] # -- die zu verwendenden Felder
out = FeatureClassToNumPyArray(in_fc, flds, explode_to_points=True)
Natürlich arbeiten Sie mit lokal gespeicherten Daten, damit alles übersichtlicher bleibt. Passen Sie daher Zeile 3 oben an Ihren Anwendungsfall an.
In Zeile 4 werden die Feld-Token angegeben, mit denen das Feld OBJECTID sowie die X- und Y-Werte der Polygonstützpunkte eingelesen werden.
Meine Regeln:
- Ich validiere die Eingabe-Feature-Class mithilfe der OGC-Geometrievalidierung, um sicherzustellen, dass Polygone mit Löchern so geordnet sind, dass die Löcher (Innenringe) auf ihren übergeordneten Außenring folgen. Mir fällt einfach kein Anwendungsfall ein, bei dem die Geometrie eines Lochs vor dem zugehörigen Außenring kommt. So muss man Löcherringe nicht mithilfe von Punkt-in-Polygon-Prüfungen ihrem Außenring zuordnen.
- Projizierte Daten und eine Projektion, die mindestens konform ist (UTM, MTM) oder in manchen Fällen flächentreu. Auf keinen Fall diese Web-Projektionen, die für die Bestimmung von Geometrieattributen nutzlos sind.
- Außenringe sind im Uhrzeigersinn orientiert, Löcher (Innenringe) gegen den Uhrzeigersinn.
- Der Anfangs- und Endpunkt eines Polygons haben dieselben X- und Y-Werte, aber unterschiedliche ID-Nummern.
- zuzüglich all der anderen Vorgaben von ESRI.
Lesen Sie Geometrie überprüfen und Geometrie reparieren für relevante Informationen.
Grundlegende Array-Ausgabe
Nun zurück zur Erzeugung des Arrays.
# ---- Fortsetzung von oben
ids = out['OID@'].astype('int') # -- ein einfaches NumPy-Array mit den Werten erstellen
xys = stu(out[['SHAPE@X', 'SHAPE@Y']]) # -- strukturiertes Array in ein ndarray umwandeln
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]
# ... gekürzt
[ 40.0, 0.0]
[ 35.0, 10.0]
[ 40.0, 20.0]
[ 40.0, 0.0]]
Zeile 2 löst die Objekt-ID-Werte aus dem abgeleiteten strukturierten Array heraus und wandelt sie in einfache Ganzzahlen um.
Zeile 3 führt einen ähnlichen Vorgang aus und übernimmt die X- und Y-Koordinaten als zweidimensionales NumPy-Array mit Gleitkommazahlen (stuist ein NumPy-Import mit Umbenennung im ersten Codeblock)
Zeile 4 verschiebt die Koordinaten in einen übersichtlicheren planaren Raum. Ich ziehe die minimalen X- und Y-Werte von allen Koordinatenwerten ab. In meinem Beispiel entspricht die Koordinate (0.0, 0.0) tatsächlich (300.000,00, 5.000.000,00). Die XY-Ebene wurde verschoben, was die Darstellung übersichtlicher und die Berechnungen einfacher macht. Wenn Sie fertig sind, verschieben Sie die Koordinaten der Berechnungen wieder zurück. Das ist also unproblematisch.
Die Zeilen 7 und 8 enthalten die ID-Werte für die erste Geometrie (1), die 9 Punkte umfasst, die den Außenring und sein Loch bilden, sowie für die 4 Punkte der letzten Geometrie (5). Die Zeilen 10–24 enthalten die Punktkoordinaten. Zeile 11 und 14 enthalten die Anfangs- und Endkoordinaten des Außenrings (Dreieck), und Zeile 15–19 die Punkte des quadratischen Lochs.
Was sich jetzt machen lässt
An diesem Punkt haben Sie eine Reihe von Punkten, die Polygongeometrien darstellen. Sie könnten
- sie verschieben (einfach Werte zu den X- und Y-Werten addieren oder davon abziehen),
- die Anzahl der Punkte entlang der Ringe erhöhen oder sie vereinfachen (dafür gibt es NumPy-Funktionen,
- und so weiter.
Wenn Sie jedoch mit einzelnen Polygonen arbeiten müssen, müssen Sie diese aufteilen. Dafür gibt es mehrere Möglichkeiten, die Polygone selbst und/oder ihre Innen- und Außenringe zu isolieren.
def array_to_rings(arr, from_to_ids=False):
"""Teilt ein Array in Außen- und Innenringe auf.
Parameter
----------
arr : array_like
Ein Nx2-Array/eine Liste/ein Geo-Array mit x- und y-Koordinaten, das
einen oder mehrere Außen- und Innenringe (Löcher) enthalten kann.
Rückgabewert
------------
Eine Liste von Arrays, gegebenenfalls verschachtelt, wenn Innenringe/Löcher vorhanden sind.
"""
# --
def _area_(a):
"""Angepasstes `_bit_area_`. Negative Flächen sind Löcher. Umgekehrt zur
`shoelace`-Formel.
"""
x0, y1 = (a.T)[:, 1:] # Kreuzprodukt wie folgt angesetzt
x1, y0 = (a.T)[:, :-1]
e0 = np.einsum('...i,...i->...i', x0, y0) # 2024-03-28 geändert
e1 = np.einsum('...i,...i->...i', x1, y1)
return np.sum((e0 - e1) * 0.5)
def _is_ccw_(a):
"""Prüfung auf Orientierung gegen den Uhrzeigersinn. Umgekehrt zur Shoelace-Formel."""
return 0 if _area_(a) > 0. else 1
# --
#
msg = "\nEin Nx2-np.ndarray ist erforderlich."
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 # -- Zähler für Außenringe
outer_ = []
inner_ = []
frto_ = [] # -- Start-End-IDs zum Aufteilen des Eingabearrays `arr`
while n < N:
ids = np.nonzero((arr[fst] == arr[fst:]).all(-1))[0]
if len(ids) >= 2:
pair = ids[:2]
else:
print("\nKein übereinstimmender Punkt gefunden, das Eingabearray ist fehlerhaft.")
return None
fst, lst = pair + n
sub = arr[fst:lst + 1]
n = n + len(sub) # Größe des Teilarrays
res = _is_ccw_(sub)
if res != 1:
outer_.append(sub)
out.append([sub])
c_o += 1
else:
inner_.append(sub)
if len(out) > 0: # Lochringe an einen vorhandenen Außenring anhängen
out[-1].append(sub)
else: # für den Fall, dass ein Loch vor einem Außenring vorkommt
out.append([sub])
frto_.append([c_o, fst, lst + 1])
fst = lst + 1
# -- je nach Fall eine Liste von Arrays oder eine Liste von Listen mit Arrays zurückgeben
out = [i[0] if len(i) == 1 else i for i in out]
if from_to_ids:
return out, frto_
return out
Und die Ausgabe:
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]])],
[... gekürzt],
[... gekürzt],
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]])]
Die obige Ausgabe zeigt, dass die erste Geometrie aus einer Liste mit zwei Arrays besteht: dem Außenring und dem Loch. Die letzten beiden Arrays gehören zu den Geometrien 4 und 5, einfachen Dreiecken.
An dieser Stelle könnten wir eine Ausgabe erzeugen, in der die Löcher entfernt werden.
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]])]
Anschließend können Sie die Arrays in arcpy-Polygone umwandeln.
def array_poly(arrs, p_type=None, sr=None):
"""Erstellt Polygon-Features aus Geo-Arrays oder ndarrays.
Wird von `geometry_fc` verwendet oder kann separat eingesetzt werden.
Parameter
----------
arr : Arrays
Punktarray(s).
p_type : Text
POLYGON oder POLYLINE
sr : Raumbezug
Raumbezugsobjekt, Name oder ID.
Hinweise
-------
Aus den Array-Daten können Linien- oder Polygon-Features erstellt werden. Die
Features können Multipart-Geometrien mit oder ohne Innenringe sein.
Außenringe sind im Uhrzeigersinn angeordnet, Innenringe (Löcher) dagegen
gegen den Uhrzeigersinn.
"""
def _arr_poly_(arr, SR, as_type):
"""Teilt das Array an den NaN-Werten auf."""
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
Anschließend lässt sich daraus mithilfe des obigen Codes eine Feature-Class erstellen.
def geometry_fc(a, p_type=None, gdb=None, fname=None, sr=None):
"""Erstellt Polygon-Features aus der von `fc_geometry` erzeugten Array-Liste.
Parameter
----------
a : Array oder Liste von Arrays
Einige können Objekt-Arrays sein, die normalerweise von ``pnts_arr`` erstellt werden.
p_type : Zeichenfolge
Geometrietyp in Großbuchstaben, z. B. POLYGON.
gdb : Text
Pfad und Name der Geodatabase.
fname : Text
Name der Feature-Class.
sr : Raumbezug
Name oder Objekt
Rückgabewert
-------
Einzelteilige und/oder mehrteilige Feature-Classes.
Hinweise
-----
Die Arbeit wird von ``array_poly`` erledigt.
"""
if p_type is None:
p_type = "POLYGON"
out = array_poly(a, p_type.upper(), sr=sr, IFT=IFT) # array_poly aufrufen
name = gdb + "/" + fname
wkspace = arcpy.env.workspace = r'memory' # veraltet: 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
Das entspricht dem Werkzeug Eliminate Polygon Part , für das man in Pro merkwürdigerweise eine Advanced-Lizenz benötigt. Warum braucht man eine Advanced-Lizenz, um ein Loch zu entfernen oder zu füllen?
Wenn Sie den Ausgabetyp oben in POLYLINE ändern, wandeln Sie Polygone in Polylinien um – genauso wie mit Polygon to Line , für das ebenfalls eine Advanced-Lizenz erforderlich ist.
Dies ist nur ein einfaches Beispiel dafür, wie Sie FeatureClassToNumPyArray für die Arbeit mit Geometrien verwenden können.
Weitere Werkzeuge und Code finden Sie unter NumPy Geometry und Tools for ArcGIS Pro.
