Comparer deux géométries dans ArcPy?

18

J'essaie de comparer deux classes d'entités distinctes pour identifier les différences entre elles (une sorte de fonction diff). Mon workflow de base:

  1. J'extrait les géométries à l'aide d'un SearchCursor
  2. Enregistrez les géométries des deux classes d'entités en tant que GeoJSON en utilisant une modification __geo_interface__(obtenue de valveLondon return {'type': 'Polygon', 'coordinates': [[((pt.X, pt.Y) if pt else None) for pt in part] for part in self]} ). C'est pour éviter l'objet de géométrie partagée qu'ESRI utilise avec les curseurs et l'impossibilité de faire des copies complètes (certaines discussions ici sur gis.stackexchange en parlent).
  3. Vérifiez les géométries des deux classes d'entités en fonction d'un identifiant unique. Par exemple, comparez la géométrie FC1 OID1 avec la géométrie FC2 OID1. Pour obtenir la géométrie en tant qu'occurrence d'objet ESRI, appelez arcpy.AsShape()(modifié pour lire les polygones avec des trous (voir point 2 ci-dessus) avec return cls(Array([map(lambda p: Point(*p) if p is not None else Point(), part) for part in coordinates])). La comparaison est simplement geom1.equals(geom2)comme indiqué dans la classe de géométrie .

Je m'attends à trouver ~ 140 changements dans les géométries, mais mon script insiste sur le fait qu'il y en a 430. J'ai essayé de vérifier ces représentations GeoJSON et elles sont identiques, mais la classe de géométrie equals () refuse de le dire.

Un exemple est ci-dessous:

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

Le comportement attendu ici doit être True (pas False).

Quelqu'un a-t-il des suggestions avant de tout déplacer vers les géométries ogr? (J'hésite car ogr.CreateGeometryFromGeoJSON () attend une chaîne, et arcpy __geo_interface__retourne un dictionnaire et j'ai l'impression d'ajouter une complexité supplémentaire).

Les ressources suivantes ont été utiles, même si elles ne répondent pas à la question:

  1. question arcpy.Geometry ici sur gis.stackexchange.com qui a été lié ci-dessus dans mon texte.
  2. Erreurs dans la classe Polygon d'arcpy sur les forums arcgis.com (apparemment il y a beaucoup d'erreurs de précision dans ArcGIS 10.0 qui ont théoriquement été corrigées dans 10.1 mais je ne peux pas vérifier que, dans 10.0 SP5, vous obtenez toujours l'erreur).
arcpy geometry Michalis Avraam
la source

Réponses:

12

Le problème est probablement lié à la précision en virgule flottante . Dans votre cas, vous avez déjà extrait les géométries à l'aide d'arcpy et vous les avez associées à votre RUID.

Heureusement, puisque vous avez installé arcpy, vous avez numpy, ce qui facilite la comparaison d'ensembles de tableaux numériques. Dans ce cas, je suggère la fonction numpy.allclose , qui est disponible dans numpy 1.3.0 (installé avec ArcGIS 10).

À partir des échantillons que vous avez donnés ci-dessus

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 atolmot-clé spécifie la valeur de tolérance.

Notez que vous ne devriez pas utiliser arcpy.AsShapedu tout. Déjà. Comme je l'ai noté dans cette question (/ plug sans vergogne), il existe un bogue connu dans ArcGIS qui tronque les géométries lorsqu'elles sont créées sans système de coordonnées (même après avoir défini la env.XYTolerancevariable d'environnement). Il arcpy.AsShapen'y a aucun moyen d'éviter cela. Heureusement, geometry.__geo_interface__il extrait les géométries correctes d'une géométrie existante (bien qu'il ne gère pas les polygones complexes sans le correctif de @JasonScheirer).

om_henners
la source
Je vous remercie. Je n'ai pas pensé à utiliser numpy pour faire ça. Une autre solution semble utiliser le module décimal et résoudre ce problème, mais cela nécessite beaucoup plus de travail.
Michalis Avraam
Je pense qu'il serait important de régler le numpy.allclose() rtolparamètre sur 0. Par défaut, c'est 1e-05 et cela peut conduire à une grande tolérance si les valeurs des tableaux sont grandes voir: stackoverflow.com/a/57063678/1914034
Sous le radar
11

La précision des coordonnées va être une considération importante ici. Les nombres à virgule flottante ne peuvent pas être stockés exactement.

Si vous utilisez l' outil de comparaison des fonctionnalités , obtient-il le résultat attendu en utilisant la tolérance XY par défaut?

blah238
la source
Je n'ai pas vérifié l'outil de comparaison des fonctionnalités car l'outil que je crée compare en fait des fonctionnalités individuelles qui se déplacent entre différentes classes d'entités. C'est-à-dire qu'une entité pourrait passer de CityRoads à CountyRoads, j'ai donc besoin de déterminer si quelque chose a changé dans la géométrie et les attributs autres que la classe d'entités qui la contient. Il existe au total 24 classes d'entités, et les entités peuvent se déplacer entre elles. La comparaison des fonctionnalités ne comparera que 2 classes d'entités, de sorte qu'elle peut me dire si elle n'existe plus dans un FC. Ensuite, je dois encore comparer la fonctionnalité pour m'assurer qu'elle n'a pas changé
Michalis Avraam
J'ai vérifié l'outil de comparaison des fonctionnalités avec la tolérance par défaut (8.983e-009 qui est assez petite mais il s'agit d'un fichier GDB) et il signale quelques changements, mais pas les bons. Plus précisément, il dit qu'il y a 69 changements de géométrie (je suppose que c'est mieux qu'avant) mais il semble supposer que l'OID est le moyen d'identifier des fonctionnalités uniques (recherche l'ancien OID1 et le nouveau OID1), ce qui n'est pas nécessairement vrai (je l'ai configuré pour utiliser mon RUID comme une sorte mais il ne l'aimait pas). Revenons donc à la planche à dessin.
Michalis Avraam
4

à côté de la réponse @ blah328, vous avez le choix de comparer deux tableaux pour signaler les différences et les similitudes avec les valeurs tabulaires et les définitions de champ avec Table Compare .

Exemple:

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
la source
Merci, je vais l'examiner lorsque j'essaierai de comparer les données d'attribut. Pour l'instant, il semble que je ne puisse pas comparer les géométries, ce qui est plus important.
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()fonction ne fonctionne pas comme prévu et / ou si les coordonnées sont légèrement modifiées dans ArcGIS, vous pouvez masser les coordonnées XY, puis comparer l'équivalent chaîne de la géométrie. Remarquez, truncateCoordinates()coupez toutes les valeurs au-delà de la 4ème décimale.

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

geom1 == geom2
klewis
la source
@ klewis- C'est une façon de comparer une géométrie, mais cela ressemble à geometry.equals (geometry) devrait retourner vrai lorsque vous comparez exactement la même géométrie. Tronquer les coordonnées est une sorte de hack dans un sens. ESRI doit peut-être commencer à utiliser le type decimal () au lieu de float s'ils ne peuvent pas gérer correctement les valeurs à virgule flottante en interne mais peuvent les représenter comme des chaînes égales.
Michalis Avraam