features with identical geometry WKT representation not equal

JamieKass

The results of arcpy.FindIdentical() with "Shape" selected as field option show two features in my dataset to be identical in geometry, yet when I run a search cursor through it, and append the WKT of the two geometries into a list, these two strings are not equal (they even have different lengths). I thought the WKT was a representation of the geometry in string form. The JSON representation is also unequal.

JamieKass

Yes, I developed something similar that creates a unique string representation of area, centroid X, and centroid Y. It is imperfect, as different datasets may require tweaking of the decimal tolerances, but it works alright the way it is.

def repGeomString(shp):
    return str(round(shp.area)) + '|' + str(round(shp.centroid.X, 2)) + '|' + str(round(shp.centroid.Y, 2))

I call this in a search cursor loop and add it to a dictionary with the OID as key, much in the same way you did it. I have experienced finding polygons with mysterious null centroids, but if you account for them somehow that's the only trip-up. This wouldn't be necessary if FindIdentical worked on large datasets. I wish the internal tolerances of this tool were exposed so we could improve these imperfect geometry estimators. I've included a sample function I've wrote that returns a dictionary relating each OID to its overlapping OIDs, based on if repGeomString is equal. Buffering often leads to small changes in shapes that should be identical, and therefore this function gets thrown off sometimes, but usually works quite well.

def fasterFindID(fc):
    """Records geometrical attributes of all features, then returns a dictionary:
        overlapDict: key = OID; value = identical overlap OIDs. The parameter fc
        is assumed to be the product of an overlap operation like Union or Intersect."""
    # this creates a representation of a geometry string by concatenating various attributes to create a "unique" identifier
    def repGeomString(shp):
        return str(round(shp.area)) + '|' + str(round(shp.centroid.X, 2)) + '|' + str(round(shp.centroid.Y, 2))

    # two dictionaries to keep track of duplicate geometries and their OIDs
    geomDictByOID = {}  # key = OID; value = represented geometry string
    uniqueGeoms = {}  # key = represented geometry string; value = list of OIDs with same geometry

    oidFieldName = arcpy.Describe(fc).oidFieldName

    nullCentroids = []

    # gather geometries, and if centroid null, append to list
    with arcpy.da.SearchCursor(fc, ("OID@", "SHAPE@")) as cursor:
        for row in cursor:
            try:
                # get full k,v pair of OID and geometries
                geomDictByOID[row[0]] = repGeomString(row[1])
                # initialize empty list as value of unique geometry key
                uniqueGeoms[repGeomString(row[1])] = []
            except:
                nullCentroids.append(row[0])

    # if null centroids exist, convert features to point and reread centroids
    # post-conversion, centroids will exist... for some reason
    if len(nullCentroids) > 0:
        # Format to string tuple for SQL query, and remove trailing comma for len = 1
        if len(nullCentroids) > 1:
            nullCentroidsFormat = str(tuple(nullCentroids))
        else:
            nullCentroidsFormat = str(tuple(nullCentroids)).replace(',', '')
        arcpy.MakeFeatureLayer_management(fc, 'nLyr', '"{0}" IN {1}'.format(oidFieldName, nullCentroidsFormat))
        nullCentroidsPts = arcpy.FeatureToPoint_management('nLyr', arcpy.Describe(fc).name + 'NullCentroids')
        # loop through null centroid features and read new centroids
        with arcpy.da.SearchCursor(nullCentroidsPts, ("OID@", "SHAPE@")) as cursor:
            for row in cursor:
                # get full k,v pair of OID and geometries
                geomDictByOID[row[0]] = repGeomString(row[1])
                # initialize empty list as value of unique geometry key
                uniqueGeoms[repGeomString(row[1])] = []

    # for each geometry key in uniqueGeoms, append all OIDs that match
    for oid, g in geomDictByOID.iteritems():
        uniqueGeoms.append(oid)

    # create dictionary that relates each oid to its overlaps
    overlapDict = {}
    for oid, g in geomDictByOID.iteritems():
        if len(uniqueGeoms) > 1:
            overlapDict[oid] = [k for k in uniqueGeoms if k != oid]

    # Delete potentially large dicts
    del geomDictByOID
    del uniqueGeoms

    return overlapDict