¿Comparando dos geometrías en ArcPy?

18

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

  1. Extraigo las geometrías usando un cursor de búsqueda
  2. Guarde las geometrías de las dos clases de entidad como GeoJSON usando un modificado __geo_interface__(lo obtuvo de valveLondon 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 usa con los cursores y la imposibilidad de hacer copias profundas (algunas discusiones aquí en gis.stackexchange hablan de ello).
  3. Verifique las geometrías de las dos clases de entidad en función de 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 arcpy.AsShape()(modificado para leer polígonos con agujeros (consulte 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 la Clase de Geometría .

Espero encontrar ~ 140 cambios en las geometrías, pero mi script insiste en que hay 430. Intenté verificar esas representaciones de GeoJSON y son idénticas, pero la Clase de Geometría igual () se niega a decirlo.

Un ejemplo está abajo:

>>> 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í debería ser Verdadero (no Falso).

¿Alguien tiene alguna sugerencia antes de mover todo a geometrías ogr? (Dudo porque ogr.CreateGeometryFromGeoJSON () espera una cadena, y arcpy's __geo_interface__devuelve un diccionario y siento que estoy agregando complejidad adicional).

Encontraron útiles los siguientes recursos, a pesar de que no responden la pregunta:

  1. arcpy. Pregunta de geometría aquí en gis.stackexchange.com que estaba vinculada anteriormente 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 en teoría se arreglaron en 10.1, pero no puedo verificar que, en 10.0 SP5 todavía se obtenga el error).
arcpy geometry Michalis Avraam
fuente

Respuestas:

12

El problema es muy probablemente uno de precisión de punto flotante . En su caso, ya ha extraído las geometrías usando arcpy, y las ha emparejado con su RUID.

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

De las muestras que diste 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

La atolpalabra clave especifica el valor de tolerancia.

Tenga en cuenta que no debe usar arcpy.AsShapeen absoluto. Nunca. Como señalé en esta pregunta (/ conector descarado), hay un error conocido en ArcGIS que trunca las geometrías cuando se crean sin un sistema de coordenadas (incluso después de establecer la env.XYTolerancevariable de entorno). En arcpy.AsShapeque no hay forma de evitar esto. Afortunadamente geometry.__geo_interface__, extrae las geometrías correctas de las geometrías existentes (aunque no maneja polígonos complejos sin la corrección de @JasonScheirer).

om_henners
fuente
Gracias. No pensé en usar numpy para hacer esto. Parece que otra solución es usar el módulo decimal y resolverlo, pero requiere mucho más trabajo.
Michalis Avraam
Creo que sería importante establecer el numpy.allclose() rtolparámetro en 0. Por defecto es 1e-05 y puede conducir a una gran tolerancia si los valores de los arrays son grandes, consulte: stackoverflow.com/a/57063678/1914034
Debajo del radar
11

La precisión de las coordenadas será una consideración importante aquí. Los números de coma flotante no se pueden almacenar exactamente.

Si utiliza la herramienta de comparación de características , ¿se obtiene el resultado esperado con la tolerancia XY predeterminada?

blah238
fuente
No verifiqué la herramienta de comparación de características ya que la herramienta que estoy construyendo en realidad compara características individuales que se mueven entre diferentes clases de características. Es decir, una entidad podría moverse de CityRoads a CountyRoads, por lo que necesito averiguar si algo cambió en la geometría y los atributos que no sean la clase de entidad que la contiene. Hay un total de 24 clases de entidad, y las entidades pueden moverse entre ellas. La comparación de características solo comparará 2 clases de entidad, por lo que me puede decir si ya no existe en un FC. Entonces todavía necesito comparar la función para asegurarme de que no haya cambiado
Michalis Avraam, el
Verifiqué la herramienta de comparación de características con la tolerancia predeterminada (8.983e-009, que es bastante pequeña pero se trata de un archivo GDB) e informa algunos cambios, pero no los correctos. Específicamente, dice que hay 69 cambios de geometría (supongo que mejor que antes), pero parece suponer que OID es la forma de identificar características únicas (busca OID1 antiguo y nuevo OID1) que no es necesariamente cierto (lo he configurado para usar mi RUID como una especie, pero no me gustó). Así que de vuelta al tablero de dibujo.
Michalis Avraam
4

junto a la respuesta @ blah328, tiene la opción de comparar dos tablas para informar diferencias y similitudes con valores tabulares y definiciones de campo con la Comparación de tablas .

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'
Aragón
fuente
Gracias. Lo investigaré cuando trate de comparar datos de atributos. Por ahora, parece que no puedo comparar geometrías, lo cual es más importante.
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

Si la .equals()función 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 de cadena de la geometría. Observe que truncateCoordinates()corta todos los valores más allá del cuarto decimal.

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

geom1 == geom2
klewis
fuente
@ klewis: esa es una forma de comparar una geometría, pero parece que geometry.equals (geometry) debería ser verdadera cuando se compara exactamente la misma geometría. Truncar las coordenadas es una especie de truco en cierto sentido. Quizás ESRI necesita comenzar a usar el tipo decimal () en lugar de flotante si no pueden manejar los valores de coma flotante internamente pero pueden representarlos como cadenas iguales.
Michalis Avraam