¿Comparar dos geometrías en ArcPy?

Question

Estoy tratando de comparar dos clases de características separadas para identificar las diferencias entre ellas (una especie de función diff). Mi flujo de trabajo básico:

1. Extraigo las geometrías utilizando un BuscarCursor
2. Guarde las geometrías de las dos clases de características como GeoJSON utilizando un __geo_interface__ (lo obtuve de válvulaLondres return {'type': 'Polygon', 'coordinates': [[((pt.X, pt.Y) if pt else None) for pt in part] for part in self]} ). Esto es para evitar el objeto de geometría compartida que ESRI utiliza con los cursores y la incapacidad de hacer copias profundas (algunas discusiones aquí en gis.stackexchange hablan de ello).
3. Compruebe las geometrías de las dos clases de características basándose en un identificador único. Por ejemplo, compare la geometría FC1 OID1 con la geometría FC2 OID1. Para obtener la geometría como una instancia de objeto ESRI, llame a arcpy.AsShape() (modificado para leer los polígonos con agujeros (véase el punto 2 anterior) con return cls(Array([map(lambda p: Point(*p) if p is not None else Point(), part) for part in coordinates])) . La comparación es simplemente geom1.equals(geom2) como se indica en el Clase de geometría .

Espero encontrar ~140 cambios en las geometrías, pero mi script insiste en que hay 430. Intenté comprobar esas representaciones GeoJSON y son idénticas, pero la clase Geometry equals() se niega a decirlo.

Un ejemplo es el siguiente:

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

El comportamiento esperado aquí debe ser True (no False).

¿Alguien tiene alguna sugerencia antes de pasar todo a geometrías ogr? (Estoy dudando, ya que ogr.CreateGeometryFromGeoJSON() espera una cadena, y arcpy __geo_interface__ devuelve un diccionario y siento que estoy añadiendo complejidad extra).

Los siguientes recursos son útiles, aunque no responden a la pregunta:

1. Pregunta sobre arcpy.Geometry aquí en gis.stackexchange.com que estaba enlazado arriba en mi texto.
2. Errores en la clase Polygon de arcpy de los foros de arcgis.com (aparentemente hay muchos errores de precisión en ArcGIS 10.0 que teóricamente se solucionaron en 10.1 pero no puedo verificarlo, en 10.0 SP5 todavía se produce el error).

Answer 1

Lo más probable es que se trate de una cuestión de precisión en coma flotante . En tu caso ya has extraído las geometrías usando arcpy, y las has emparejado con tu RUID.

Afortunadamente, ya que tienes arcpy instalado, tienes numpy, que hace que la comparación de conjuntos de matrices numéricas sea fácil. En este caso yo sugeriría el numpy.allclose que está disponible en numpy 1.3.0 (instalado con ArcGIS 10).

A partir de las muestras que has dado arriba

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

Le site atol especifica el valor de la tolerancia.

Tenga en cuenta que no debe utilizar arcpy.AsShape en absoluto. Siempre. Como señalé en esta pregunta (/enchufe sin vergüenza) hay un error conocido en ArcGIS que trunca las geometrías cuando se crean sin un sistema de coordenadas (incluso después de establecer el env.XYTolerance variable de entorno). En arcpy.AsShape no hay manera de evitarlo. Afortunadamente geometry.__geo_interface__ extrae las geometrías correctas de una geometría existente (aunque no maneja polígonos complejos sin el arreglo de @JasonScheirer ).

Answer 2

La precisión de las coordenadas va a ser una consideración importante aquí. Los números en coma flotante no pueden almacenarse con exactitud.

Si utiliza el Comparación de características ¿se obtiene el resultado esperado utilizando la tolerancia XY por defecto?

Answer 3

Además de la respuesta de @blah328, tiene la opción de comparar dos tablas para informar de las diferencias y similitudes con los valores tabulares y las definiciones de campo con Tabla de comparación .

Ejemplo:

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'

Answer 4

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

Si el .equals() no funciona como se esperaba y/o las coordenadas están ligeramente alteradas en ArcGIS, puede masajear las coordenadas XY, y luego comparar el equivalente en cadena de la geometría. Aviso, truncateCoordinates() corta todos los valores más allá del cuarto decimal.

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

geom1 == geom2

Answer 5

Puede utilizar el Seleccionar capa por ubicación (gestión de datos) con el parámetro de tipo de solapamiento "ARE_IDENTICAL_TO", cambiar la selección Compruebe el recuento de filas y, a continuación, realice un bucle a través de las filas para recoger los objectids o cualquier otra información pertinente.