Comparing two geometries in ArcPy?

Question

I am trying to compare two separate feature classes to identify differences between them (sort of a diff function). My basic workflow:

1. I extract the geometries using a SearchCursor
2. Save the geometries of the two feature classes as GeoJSON using a modified __geo_interface__ (got it from valveLondon return {'type': 'Polygon', 'coordinates': [[((pt.X, pt.Y) if pt else None) for pt in part] for part in self]}). This is to avoid the shared geometry object that ESRI uses with cursors and the inability to make deep copies (some discussions here on gis.stackexchange talk about it).
3. Check the geometries of the two feature classes based on a unique identifier. For example, compare the FC1 OID1 geometry with the FC2 OID1 geometry. To get the geometry as an ESRI object instance, call arcpy.AsShape() (modified to read polygons with holes (see point 2 above) with return cls(Array([map(lambda p: Point(*p) if p is not None else Point(), part) for part in coordinates])). The comparison is simply geom1.equals(geom2) as indicated in the Geometry Class.

I expect to find ~140 changes in geometries, but my script insists there are 430. I tried to check those GeoJSON representations and they are identical, yet the Geometry Class equals() refuses to say so.

An example is below:

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

The expected behavior here should be True (not False).

Does anyone have any suggestions before I move everything to ogr geometries? (I am hesitant as ogr.CreateGeometryFromGeoJSON() expects a string, and arcpy's __geo_interface__ returns a dictionary and I feel like I am adding extra complexity).

Found the following resources helpful, even though they do not answer the question:

1. arcpy.Geometry question here on gis.stackexchange.com which was linked above in my text.
2. Errors in arcpy's Polygon class from the arcgis.com forums (apparently there are a lot of precision errors in ArcGIS 10.0 which theoretically got fixed in 10.1 but I cannot verify that, in 10.0 SP5 you still get the error).

Answer 1

The issue is most likely one of floating point precision. In your case you've already extracted the geometries using arcpy, and you've matched them with your RUID.

Happily since you've got arcpy installed you've got numpy, which makes comparing sets of numeric arrays easy. In this case I'd suggest the numpy.allclose function, which is available in numpy 1.3.0 (installed with ArcGIS 10).

From the samples you gave above

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

The atol keyword specifies the tolerance value.

Note that you shouldn't use arcpy.AsShape at all. Ever. As I noted in this question (/shameless plug) there's a known bug in ArcGIS that truncates geometries when they are created without a coordinate system (even after setting the env.XYTolerance environment variable). In arcpy.AsShape there's no way to avoid this. Fortunately geometry.__geo_interface__ does extract the correct geometries from an existing geometries (though it doesn't handle complex polygons without the fix from @JasonScheirer ).

Answer 2

The coordinate precision is going to be an important consideration here. Floating point numbers cannot be stored exactly.

If you use the Feature Compare tool, does it come up with the expected result using the default XY tolerance?

Answer 3

beside @blah328 answer, you hava choice to compare two tables for reporting differences and similarities with tabular values and field definitions with Table Compare.

Example:

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

If the .equals() function is not working as expected and/or the coordinates are slightly altered in ArcGIS, you can massage the XY coordinates, then compare the String equivalent of the Geometry. Notice, truncateCoordinates() chops off all values beyond the 4th decimal place.

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

geom1 == geom2

Answer 5

You could use the Select Layer By Location (Data Management) tool with the "ARE_IDENTICAL_TO" overlap_type parameter, switch the selection, check the row count and then loop through the rows to collect the objectids or any other pertinent information.