Porównywanie dwóch geometrii w ArcPy?

Próbuję porównać dwie oddzielne klasy funkcji, aby zidentyfikować różnice między nimi (rodzaj funkcji diff). Mój podstawowy obieg pracy:

1. Wyodrębniam geometrie za pomocą SearchCursor
2. Zapisz geometrie dwóch klas elementów jako GeoJSON, używając zmodyfikowanego __geo_interface__( pobranego z ValveLondonreturn {'type': 'Polygon', 'coordinates': [[((pt.X, pt.Y) if pt else None) for pt in part] for part in self]} ). Ma to na celu uniknięcie współdzielonego obiektu geometrii, którego ESRI używa z kursorami, oraz niemożności wykonywania głębokich kopii (niektóre dyskusje na ten temat na gis.stackexchange mówią o tym).
3. Sprawdź geometrie dwóch klas elementów na podstawie unikalnego identyfikatora. Na przykład porównaj geometrię FC1 OID1 z geometrią FC2 OID1. Aby uzyskać geometrię jako instancję obiektu ESRI, wywołaj arcpy.AsShape()(zmodyfikowano tak, aby czytała wielokąty z otworami (patrz punkt 2 powyżej) za pomocą return cls(Array([map(lambda p: Point(*p) if p is not None else Point(), part) for part in coordinates])). Porównanie jest po prostu geom1.equals(geom2)wskazane w klasie geometrii .

Spodziewam się znaleźć ~ 140 zmian w geometrii, ale mój skrypt nalega na 430. Próbowałem sprawdzić te reprezentacje GeoJSON i są one identyczne, ale Klasa Geometrii równa się () nie chce tego powiedzieć.

Przykład jest poniżej:

>>> 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

Oczekiwane zachowanie tutaj powinno być prawdziwe (nie fałszywe).

Czy ktoś ma jakieś sugestie, zanim przeniosę wszystko do geometrii ogr? (Waham się, ponieważ ogr.CreateGeometryFromGeoJSON () oczekuje łańcucha, a arcpy __geo_interface__zwraca słownik i mam wrażenie, że dodam dodatkową złożoność).

Pomocne okazały się następujące zasoby, mimo że nie odpowiadają na pytanie:

1. arcpy.Geometry pytanie tutaj na gis.stackexchange.com, które zostało połączone powyżej w moim tekście.
2. Błędy w klasie Polygon arcpy z forów arcgis.com (podobno w ArcGIS 10.0 istnieje wiele błędów precyzji, które teoretycznie zostały naprawione w 10.1, ale nie mogę tego zweryfikować, w 10.0 SP5 nadal pojawia się błąd).

Problem najprawdopodobniej dotyczy precyzji zmiennoprzecinkowej . W twoim przypadku wyodrębniłeś już geometrie za pomocą arcpy i dopasowałeś je do swojego RUID.

Na szczęście, odkąd masz zainstalowany arcpy, masz numpy, co ułatwia porównywanie zestawów tablic numerycznych. W takim przypadku zasugerowałbym funkcję numpy.allclose , która jest dostępna w numpy 1.3.0 (instalowanej z ArcGIS 10).

Z próbek, które podałeś powyżej

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

Słowo atolkluczowe określa wartość tolerancji.

Pamiętaj, że nie powinieneś w ogóle używać arcpy.AsShape. Zawsze. Jak zauważyłem w tym pytaniu (/ bezwstydna wtyczka), w ArcGIS jest znany błąd, który obcina geometrie, gdy są tworzone bez układu współrzędnych (nawet po ustawieniu env.XYTolerancezmiennej środowiskowej). W arcpy.AsShapenie ma sposobu, aby tego uniknąć. Na szczęście geometry.__geo_interface__wyodrębnia prawidłowe geometrie z istniejących geometrii (choć nie obsługuje złożonych wielokątów bez poprawki z @JasonScheirer).

Ważna będzie tutaj dokładność współrzędnych. Liczb zmiennoprzecinkowych nie można dokładnie zapisać.

Jeśli korzystasz z narzędzia Porównanie funkcji , czy przyniesie on oczekiwany wynik przy użyciu domyślnej tolerancji XY?

oprócz odpowiedzi @ blah328 masz możliwość porównania dwóch tabel pod kątem zgłaszania różnic i podobieństw z wartościami tabelarycznymi i definicjami pól za pomocą Tabeli porównania .

Przykład:

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' 

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

Jeśli .equals()funkcja nie działa zgodnie z oczekiwaniami i / lub współrzędne są nieco zmienione w ArcGIS, możesz masować współrzędne XY, a następnie porównać ciąg znaków w geometrii. Zauważ, że truncateCoordinates()wycina wszystkie wartości poza czwartym miejscem po przecinku.

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

geom1 == geom2

Możesz użyć narzędzia Wybierz warstwę według lokalizacji (zarządzanie danymi) z parametrem „ARE_IDENTICAL_TO” parametru „overlap_type”, przełączać zaznaczenie , sprawdzać liczbę wierszy, a następnie zapętlać rzędy w celu zebrania obiektów objectid lub innych istotnych informacji.