19

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

5 Answers 5

14
+100

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

2
  • Thank you. I didn't think of using numpy to do this. Another solution seems to be using the decimal module and work through that, but it requires a lot more work. Nov 13, 2012 at 19:58
  • I think it would be important to set the numpy.allclose() rtol parameter to 0. By default it's 1e-05 and it can lead to a large tolerance if the arrays values are large see: stackoverflow.com/a/57063678/1914034 Jul 16, 2019 at 19:06
11

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?

2
  • I did not check the Feature Compare tool as the tool I am building actually compares individual features that move between different feature classes. That is to say, a feature could move from CityRoads to CountyRoads, so I need to figure out if anything changed in the geometry and the attributes other than the feature class holding it. There are a total of 24 feature classes, and features can move between them. Feature Compare will compare only 2 feature classes, so it can tell me if it no longer exists in a FC. Then I still need to compare the feature to ensure it didn't change Nov 7, 2012 at 16:49
  • I checked the Feature Compare tool with the default tolerance (8.983e-009 which is quite small but this is a File GDB) and it reports some changes, but not the right ones. Specifically, it says there are 69 geometry changes (I guess better than before) but it seems to assume that OID is the way to identify unique features (searches old OID1 and new OID1) which is not necessarily true (I have set it to use my RUID as a sort but it didn't like it). So back to the drawing board. Nov 7, 2012 at 17:12
4

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' 
1
  • Thank you I will look into it when I try to compare attribute data. For now, it looks like I cannot compare geometries, which is more important. Nov 7, 2012 at 17:14
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

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
1
  • @klewis- That is one way to compare a geometry, but it feels like geometry.equals(geometry) should return true when you are comparing the exact same geometry. Truncating the coordinates is kind of a hack in a sense. Perhaps ESRI needs to start using the decimal() type instead of float if they can't handle floating point values correctly internally but can represent them as equal strings. Nov 12, 2012 at 17:19
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.