Vergleichen von zwei Geometrien in ArcPy?

18

Ich versuche, zwei separate Feature-Classes zu vergleichen, um Unterschiede zwischen ihnen zu identifizieren (eine Art Diff-Funktion). Mein grundlegender Workflow:

  1. Ich extrahiere die Geometrien mit einem SearchCursor
  2. Speichern Sie die Geometrien der beiden Feature-Classes als GeoJSON mit einer geänderten __geo_interface__(von valveLondon return {'type': 'Polygon', 'coordinates': [[((pt.X, pt.Y) if pt else None) for pt in part] for part in self]} ). Dies soll das gemeinsame Geometrieobjekt vermeiden, das ESRI mit Cursorn verwendet, und die Unfähigkeit, tiefe Kopien zu erstellen (einige Diskussionen hier auf gis.stackexchange sprechen darüber).
  3. Überprüfen Sie die Geometrien der beiden Feature-Classes anhand eines eindeutigen Bezeichners. Vergleichen Sie beispielsweise die FC1-OID1-Geometrie mit der FC2-OID1-Geometrie. Um die Geometrie als ESRI-Objektinstanz abzurufen, rufen Sie arcpy.AsShape()(geändert, um Polygone mit Löchern zu lesen (siehe Punkt 2 oben) mit return cls(Array([map(lambda p: Point(*p) if p is not None else Point(), part) for part in coordinates])). Der Vergleich erfolgt einfach geom1.equals(geom2)wie in der Geometrieklasse angegeben .

Ich erwarte etwa 140 Änderungen in den Geometrien, aber mein Skript besteht darauf, dass es 430 gibt. Ich habe versucht, diese GeoJSON-Darstellungen zu überprüfen, und sie sind identisch, aber die Geometrieklasse equals () weigert sich, dies zu sagen.

Ein Beispiel ist unten:

>>> geom1geoJSON 
{'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
>>> geom2geoJSON 
{'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
>>> geom1 = arcpy.AsShape(geom1geoJSON)
>>> geom2 = arcpy.AsShape(geom2geoJSON)
>>> geom1.equals(geom2)
False
>>> geom2.equals(geom1)
False

Das erwartete Verhalten sollte hier True (nicht False) sein.

Hat jemand irgendwelche Vorschläge, bevor ich alles auf OGR-Geometrien verschiebe? (Ich zögere, da ogr.CreateGeometryFromGeoJSON () einen String erwartet und arcpy's __geo_interface__ein Wörterbuch zurückgibt und ich das Gefühl habe, zusätzliche Komplexität hinzuzufügen.)

Fanden die folgenden Ressourcen hilfreich, obwohl sie die Frage nicht beantworten:

  1. arcpy.Geometry Frage hier auf gis.stackexchange.com, die oben in meinem Text verlinkt wurde.
  2. Fehler in der Polygon-Klasse von arcpy aus den arcgis.com-Foren (anscheinend gibt es viele Präzisionsfehler in ArcGIS 10.0, die theoretisch in 10.1 behoben wurden, aber ich kann nicht überprüfen, dass in 10.0 SP5 der Fehler weiterhin auftritt).
arcpy geometry Michalis Avraam
quelle

Antworten:

12

Das Problem ist höchstwahrscheinlich die Gleitkommapräzision . In Ihrem Fall haben Sie die Geometrien bereits mit arcpy extrahiert und mit Ihrer RUID abgeglichen.

Seit Sie arcpy installiert haben, haben Sie glücklicherweise numpy, was das Vergleichen von Sätzen numerischer Arrays vereinfacht. In diesem Fall würde ich die Funktion numpy.allclose vorschlagen , die in numpy 1.3.0 (installiert mit ArcGIS 10) verfügbar ist.

Aus den Proben, die Sie oben angegeben haben

geom1geoJSON = {'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
geom2geoJSON = {'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}

import numpy as np

close = np.allclose(np.array(geom1geoJSON["coordinates"]), np.array(geom2geoJSON["coordinates"]), atol=1e-7)
#Returns True

Das atolSchlüsselwort gibt den Toleranzwert an.

Beachten Sie, dass Sie überhaupt nicht verwenden sollten arcpy.AsShape. Je. Wie in dieser Frage (/ shameless plug) erwähnt, gibt es in ArcGIS einen bekannten Fehler, der Geometrien abschneidet, wenn sie ohne Koordinatensystem erstellt werden (auch nach dem Festlegen der env.XYToleranceUmgebungsvariablen). In arcpy.AsShapegibt es keine Möglichkeit , dies zu vermeiden. Zum Glück geometry.__geo_interface__werden die korrekten Geometrien aus vorhandenen Geometrien extrahiert (obwohl ohne die Korrektur von @JasonScheirer keine komplexen Polygone verarbeitet werden).

om_henners
quelle
Vielen Dank. Ich habe nicht daran gedacht, dies mit Numpy zu tun. Eine andere Lösung scheint darin zu bestehen, das Dezimalmodul zu verwenden und es durchzuarbeiten, aber es erfordert viel mehr Arbeit.
Michalis Avraam
Ich denke, es wäre wichtig, den numpy.allclose() rtolParameter auf 0 zu setzen. Standardmäßig ist es 1e-05 und es kann zu einer großen Toleranz führen, wenn die Array-Werte groß sind. Siehe: stackoverflow.com/a/57063678/1914034
Below the Radar
11

Die Koordinatengenauigkeit wird hier eine wichtige Rolle spielen. Gleitkommazahlen können nicht exakt gespeichert werden.

Wenn Sie das Feature Compare- Tool verwenden, wird das erwartete Ergebnis mit der Standard-XY-Toleranz erzielt?

blah238
quelle
Ich habe das Feature Compare-Tool nicht überprüft, da das von mir erstellte Tool tatsächlich einzelne Features vergleicht, die zwischen verschiedenen Feature-Classes verschoben werden. Das heißt, ein Feature könnte von CityRoads nach CountyRoads verschoben werden. Ich muss also herausfinden, ob sich etwas an der Geometrie und den Attributen geändert hat, außer an der Feature-Class, in der es enthalten ist. Es gibt insgesamt 24 Feature-Classes, zwischen denen Features verschoben werden können. Feature Compare vergleicht nur 2 Feature-Classes, sodass ich feststellen kann, ob es in einer FC nicht mehr vorhanden ist. Dann muss ich das Feature noch vergleichen, um sicherzustellen, dass es sich nicht ändert
Michalis Avraam
Ich habe das Feature Compare-Tool mit der Standardtoleranz (8.983e-009, die recht klein ist, aber eine Datei-GDB ist) überprüft und es werden einige Änderungen gemeldet, jedoch nicht die richtigen. Konkret heißt es, dass es 69 Geometrieänderungen gibt (ich denke besser als zuvor), aber es scheint anzunehmen, dass die OID die Art und Weise ist, eindeutige Merkmale zu identifizieren (sucht nach alter OID1 und neuer OID1), was nicht unbedingt der Wahrheit entspricht (ich habe sie so eingestellt, dass sie verwendet wird) meine RUID als eine Art, aber es hat nicht gefallen). Also zurück zum Zeichenbrett.
Michalis Avraam
4

Neben der Antwort @ blah328 haben Sie die Wahl, zwei Tabellen zu vergleichen, um Unterschiede und Ähnlichkeiten mit Tabellenwerten und Felddefinitionen mit Table Compare zu melden .

Beispiel:

import arcpy
from arcpy import env
arcpy.TableCompare_management(r'c:\Workspace\wells.dbf', r'c:\Workspace
\wells_new.dbf', 'WELL_ID', 'ALL', 'IGNORE_EXTENSION_PROPERTIES', 'WELL_DEPTH 0.001',
'#','CONTINUE_COMPARE', r'C:\Workspace\well_compare.txt'
Aragon
quelle
Vielen Dank, ich werde es untersuchen, wenn ich versuche, Attributdaten zu vergleichen. Im Moment kann ich anscheinend keine Geometrien vergleichen, was wichtiger ist.
Michalis Avraam
3 
def truncateCoordinates(myGeometry)
    trucated_coords = []
    partnum = 0

    for part in (myGeometry):
        for pnt in myGeometry.getPart(partnum):
            if pnt:
                trucated_coords.append("{:10.4f}".format(pnt.X))
                trucated_coords.append("{:10.4f}".format(pnt.Y))
             else:
                continue
        partnum += 1     
    return truncated_coords

Wenn die .equals()Funktion nicht wie erwartet funktioniert und / oder die Koordinaten in ArcGIS geringfügig geändert werden, können Sie die XY-Koordinaten massieren und dann das Zeichenfolgenäquivalent der Geometrie vergleichen. Beachten Sie, dass truncateCoordinates()alle Werte nach der 4. Dezimalstelle abgeschnitten werden.

geom1 = truncateCoordinates(feature1.Shape)
geom2 = truncateCoordinates(feature2.Shape)

geom1 == geom2
klewis
quelle
@klewis- Das ist eine Möglichkeit, eine Geometrie zu vergleichen, aber es fühlt sich an wie Geometrie. Gleichungen (Geometrie) sollten true zurückgeben, wenn Sie genau dieselbe Geometrie vergleichen. Das Abschneiden der Koordinaten ist in gewissem Sinne eine Art Hack. Möglicherweise muss ESRI den Typ decimal () anstelle von float verwenden, wenn sie Gleitkommawerte intern nicht korrekt verarbeiten können, sie aber als gleiche Zeichenfolgen darstellen können.
Michalis Avraam
1

Sie können das Werkzeug Layer nach Speicherort auswählen (Datenverwaltung) mit dem Parameter "ARE_IDENTICAL_TO" verwenden, die Auswahl ändern, die Zeilenzahl überprüfen und dann die Zeilen durchlaufen, um die Objekt-IDs oder andere relevante Informationen zu sammeln.

Brent Edwards
quelle