geometry - Comparing two geometries in ArcPy?

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

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