AnsweredAssumed Answered

arcpy - find overlap between polygons and remove the less preferable one

Question asked by Cheng_ANU on Nov 26, 2018

Hi guys, 

 

I have a list of polygons, each with x&y coordinates, and an attribute called WR_ratio. Some of the polygons are overlapped while some are not. What I need to do is to find every pair of overlapped polygons and remove the one with lower WR_ratio. I'm currently doing this:

 

    rmvl = [] # A list containing the indices of the polygons that need to be removed
    for i in range(len(resfc)-1): # "resfc" is the list of polygons
        if i in rmvl:
            continue    
        for j in range(i+1, len(resfc)):
            if j in rmvl:
                continue

            # Get the latitudes and longitudes of i, j
            lati = float(InfoArray[i][1]) # "InfoArray" contains all the attributes of the polygons
            loni = float(InfoArray[i][2])
            latj = float(InfoArray[j][1])
            lonj = float(InfoArray[j][2])

            # Accelerate the Combinations C(n, r)
            if abs(lati - latj)>0.05 or abs(loni - lonj)>0.05:
                continue
           
            # Calculate the intersection
            arcpy.Intersect_analysis(in_features=[resfc[i], resfc[j]], out_feature_class="ijintersect.shp", join_attributes="ALL", cluster_tolerance="Unknown", output_type="INPUT")
            arcpy.AddField_management("ijintersect.shp","area","Double")
            expression1 = "{0}".format("!SHAPE.area@SQUAREKILOMETERS!")       
            arcpy.CalculateField_management("ijintersect.shp", "area", expression1, "PYTHON")
            cursor = arcpy.SearchCursor("ijintersect.shp")
            area = [row.getValue("area") for row in cursor]

            # Compare the water-rock ratio of two polygons
            if area==[]:
                continue
            else:
                ratioi = float(InfoArray[i][3])
                ratioj = float(InfoArray[j][3])

                # Append the "worse" one in rmvl
                if ratioi>=ratioj:
                    rmvl.append(j)
                    print (str(resfc[j]) + " removed")
                else:
                    rmvl.append(i)
                    print (str(resfc[i]) + " removed")
                    break
               
    print (cell,": Remove (removal index): " + str(rmvl))

    # Extract the "better" polygon
    resl = [x for x in resfc if resfc.index(x) not in rmvl]

 

However, if there are 100 polygons I will probably have to calculate the intersection 100*100 times, which is extremely inefficient. Ignoring the pair of polygons that are far away from each other helps a little bit, but it's still not ideal. I know there are many genius here. Can anyone think of a smarter way of doing this? 

 

Any help would be much appreciated! 

Outcomes