topic Re: Python Code issue in Transportation Questions

Python Code issue

Anonymous User — Sat, 11 Dec 2021 14:02:34 GMT

Original User: luke.kaim

Hi Everyone,

I am trying to write python code to detect if there is overlapping polygon from two feature classes. The goal of this code is to get overlapping polygons and there area. I start by trying to get the area of the polygons themselves. This I am trying to do using arcpy.geometry. Next I want to use the overlap geometry operation in 10.1. This is a boolean operation that will speed up the intersection command later. Then I do the intersection and the symmetric difference. The end goal is to get polygon overlap/difference.

I think I have a few bugs in my code, but I am struggling to get the logic and the syntax correct. I tried to model my code off of this thread: http://forums.arcgis.com/threads/60968-area-of-overlap-between-polygons-in-two-different-feature-classes?highlight=arcpy.Geometry%28%29 .

Any help would be greatly appreciated.
Code:

def main():
    try:
        import arcpy, sys, traceback, os, string, locale
        arcpy.env.overwriteOutput = True

        from arcpy import env
        env.workspace =r"C:\Users\Luke Kaim\Documents\University of Maine\Fall_2012\Volunteer Geographic Information\simalar"
        
        VGIFile = r"C:\Users\Luke Kaim\Documents\University of Maine\Fall_2012\Volunteer Geographic Information\simalar\polygon.shp"
        VGIFileShapeName = arcpy.Describe(VGIFile).ShapeFieldName
        VGIRows = arcpy.SearchCursor(VGIFile)
        

        GNIS_FC = r"C:\Users\Luke Kaim\Documents\University of Maine\Fall_2012\Volunteer Geographic Information\simalar\GNIS.shp"
        GNIS_FCShapeName = arcpy.Describe(GNIS_FC).ShapeFieldName
        GNIS_FCRows = arcpy.SearchCursor(GNIS_FC)

        
        #Create an empty Geometry object
        g = arcpy.Geometry()
        area=0
        # get Geometry of every VGI polygon. Then get the area of each polygon
            for  VGIFileShapeName in VGIRows:
                feat =  VGIFileShapeName.getValue(shapeName)
                area += feat.area
print area


        #Create an empty Geometry object
        g = arcpy.Geometry()
        area=0
        # get Geometry of every GNIS polygon. Then get the area of each polygon
             for  GNIS_FCShapeName in GNIS_FCRows:
                feat =  GNIS_FCShapeName.getValue(shapeName)
                area += feat.area
print area





                for VGIRow in VGIRows:
                    VGIGEOM = VGIRow.getValue(VGIFileShapeName)


                    GNIS_Rows = arcpy.SearchCursor(GNIS_FC)
                    for GNIS_Row in GNIS_Rows
                        GNIS_Geom = GNIS_Row.getValue(GNIS_FCShapeName)
                        if GNIS_Geom.overlaps(VGIGEOM):
                            arcpy.Intersect_analysis([GNIS_Geom,VGIGEOM],"in_memory/Intersection")
                            arcpy.SymDiff_analysis([GNIS_Geom,VGIGEOM],"in_memory/SymDiff")
                            
                            intersectionRows = arcpy.SearchCursor("in_memory/Intersection")
                    for intersectionRow in intersectionRows:
                        intersectionArea = intersectionRow.Shape.area
                        print "Overlap Area: " + str(intersectionArea)


                            intersectionRows2 = arcpy.SearchCursor("in_memory/"in_memory/SymDiff"")
                    for intersectionRow1 in intersectionRows2:
                            intersectionArea2 = intersectionRow1.Shape.area
                        print "Overlap Area :in_memory/SymDiff " + str(intersectionArea2)
                        

        del   VGIRows, GNIS_FCRows,GNIS_Rows

Thank you,
Luke Kaim

Thank you
Luke Kaim (SIE)
Lucas.Kaim@maine.edu
(914)263-7866
"Don�??t complain. Just work harder" (Randy Pausch).

Re: Python Code issue

Anonymous User — Sat, 12 Jan 2013 18:39:53 GMT

Original User: luke.kaim

I am still struggling with this code. Does anyone have suggestions?

Re: Python Code issue

Anonymous User — Sat, 11 Dec 2021 14:02:37 GMT

I'm not sure I understand what you are trying to do here. I would just intersect first, then try the SymDiff. But to do what you are trying to do you could try something like this (if you just want to print the areas):

def GetOverlap(infc_a,infc_b):
    try:
        arcpy.overwriteOutput = True
      
        # Geom object A
        Geom_a = arcpy.Geometry()
        Geom_a_list = arcpy.CopyFeatures_management(infc_a,Geom_a)

        # Geom object B
        Geom_b = arcpy.Geometry()
        Geom_b_list = arcpy.CopyFeatures_management(infc_b,Geom_b)

        # Find Overlaps    
        for geom in Geom_a_list:
            for geomb in Geom_b_list:
                if geom.overlaps(geomb):
                    arcpy.Intersect_analysis([geom, geomb], r'in_memory\Intersection')
                    arcpy.SymDiff_analysis(geom,geomb,r'in_memory\SymDiff')
                    Rows = arcpy.SearchCursor(r'in_memory\Intersection')
                    for Row in Rows:
                        intersectionArea = Row.Shape.area
                    Rows2 = arcpy.SearchCursor(r'in_memory\SymDiff')
                    for r in Rows2:
                        symArea = r.Shape.area
                    print "Overlap Area: %s" %str(intersectionArea/43560) # in acres
                    print "SymDiff Area: %s\n"% str(symArea/43560) # in acres
        
        print 'done'
       
        
    except:
        tb = sys.exc_info()[2]
        tbinfo = traceback.format_tb(tb)[0]
        pymsg = "PYTHON ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n     " +  \
                str(sys.exc_type) + ": " + str(sys.exc_value) + "\n"
        msgs = "ARCPY ERRORS:\n" + arcpy.GetMessages(2) + "\n"

        arcpy.AddError(msgs)
        arcpy.AddError(pymsg)

        print msgs
        print pymsg
        
        arcpy.AddMessage(arcpy.GetMessages(1))
        print arcpy.GetMessages(1)



if __name__ == '__main__':  # This will use the variables in this script
                
    import arcpy, sys, traceback
    arcpy.env.overwriteOutput = True

    from arcpy import env
    env.workspace =r"C:\Users\Luke Kaim\Documents\University of Maine\Fall_2012\Volunteer Geographic Information\simalar"
    
    polygon = r"C:\Users\Luke Kaim\Documents\University of Maine\Fall_2012\Volunteer Geographic Information\simalar\polygon.shp"
    GNIS = r"C:\Users\Luke Kaim\Documents\University of Maine\Fall_2012\Volunteer Geographic Information\simalar\GNIS.shp"
    
    try:
        GetOverlap(polygon,GNIS)
        
    except:
        tb = sys.exc_info()[2]
        tbinfo = traceback.format_tb(tb)[0]
        pymsg = "PYTHON ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n     " +        str(sys.exc_type) + ": " + str(sys.exc_value) + "\n"
        msgs = "ARCPY ERRORS:\n" + arcpy.GetMessages(2) + "\n"

        arcpy.AddError(msgs)
        arcpy.AddError(pymsg)

        print msgs
        print pymsg
        
        arcpy.AddMessage(arcpy.GetMessages(1))
        print arcpy.GetMessages(1)

Note, I displayed area in acres because the areas of my test FCs were in feet. Also, look at the indentation. There were a few indentation errors in yours.