I am trying to compare two separate feature classes to identify differences between them (sort of a diff function). My basic workflow:
- I extract the geometries using a SearchCursor
- Save the geometries of the two feature classes as GeoJSON using a modified
__geo_interface__
(got it from valveLondonreturn {'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). - 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) withreturn cls(Array([map(lambda p: Point(*p) if p is not None else Point(), part) for part in coordinates]))
. The comparison is simplygeom1.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:
- arcpy.Geometry question here on gis.stackexchange.com which was linked above in my text.
- 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).