area of overlap between polygons in two different feature classes

JohnMorgan · ‎06-26-2012

Hello ArcPy folks,

I am working on a script to get the area of overlap between polygons in two different feature classes. I am able to loop through one polygon by polygon while looping through the other. But, I am stuck at how to get the area of overlap. Any suggestions would be awesome. Here is my code:


import arcpy

arcpy.env.workspace = "C:/path/PADUS/Polys.gdb"

countiesFC = "Counties"
countyRows = arcpy.SearchCursor(countiesFC) 
countyFields = arcpy.ListFields(countiesFC)

PADUSFC = "PADUS"
PADUSRows = arcpy.SearchCursor(PADUSFC) 
PADUSFields = arcpy.ListFields(PADUSFC)

countyGeom = arcpy.Geometry()
countyGeomList = arcpy.CopyFeatures_management("Counties", countyGeom)
PADUSGeom = arcpy.Geometry()
PADUSGeomList = arcpy.CopyFeatures_management("PADUS", PADUSGeom)

countyIndex = 0
#arcpy.AddMessage(str(len(countyGeomList)))
for countyRow in countyRows:
    countyGeom = countyGeomList[countyIndex]
    #arcpy.AddMessage("county: " + countyRow.name + " - " + str(countyGeom.area))
    PADUSFC = "PADUS"
    PADUSRows = arcpy.SearchCursor(PADUSFC) 
    PADUSIndex = 0
    for PADUSRow in PADUSRows:
        PADUSGeom = PADUSGeomList[PADUSIndex]
        if PADUSGeom.overlaps(countyGeom):
            arcpy.AddMessage("PADUS: " + PADUSRow.own_type + "overlaps " + countyRow.NAME)
            #here want to find the area of overlap between the two geometry polygons
        PADUSIndex = PADUSIndex + 1
        #arcpy.Delete_management("C:/temp/test.shp")
    del PADUSFC
    del PADUSRow
    del PADUSRows
    #arcpy.AddMessage(str(countyIndex))
    countyIndex = countyIndex + 1

Thanks,
Derek

WeifengHe · ‎06-26-2012

Don't know much about python, but I think the log is the same as the other language.

When you identify 2 polygons overlap, you can get the intersection of the 2 polygons, and get the area of the intersection.

Here is the code snippet in ArcObject.

ITopologicalOperator to = PADUSGeom;
IArea a = to.Intersect(countyGeom, esriGeometry2Dimension);
double area = a.Area;

JohnMorgan · ‎06-26-2012

Weifeng,

I think you are pointing me in the right direction. However, the syntax must be different for python as you suggested. I have tried the following which appears as if it is going to run. However, it never seems to finish. At least it is been running for about an hour now. Anyone out there know what I may be doing wrong?

a = arcpy.Geometry()
a = arcpy.Intersect_analysis(PADUSGeom, countyGeom, "", "", "")
arcpy.AddMessage("Overlap Area: " + str(a.area))

Cheers,
Derek

NobbirAhmed · ‎06-26-2012

Change your if block as follows:

if PADUSGeom.overlaps(countyGeom):

    diff = PADUSGeom.difference(countyGeom)
    diff_area = diff.area

difference is available in 10.1.

MarcinGasior · ‎06-27-2012

if PADUSGeom.overlaps(countyGeom):

    diff = PADUSGeom.difference(countyGeom)
    diff_area = diff.area

That's great news, the Geometry class has now wonderful methods but... in 10.1

Derek,
If you're still in 10.0, Intersect_analysis tool is reasonable choise, I think.
You need to provide as inputs a list of geometries and an output is a feature class. To make calculations fast, this feature class can be written to 'in_memory' workspace:

if PADUSGeom.overlaps(countyGeom):
    arcpy.Intersect_analysis([countyGeom, PADUSGeom], "in_memory/Intersection")
    intersectionRows = arcpy.SearchCursor("in_memory/Intersection")
    for intersectionRow in intersectionRows:
        intersectionArea = intersectionRow.Shape.area
        print "Overlap Area: " + str(intersectionArea)

    #delete in_memory feature class
    arcpy.Delete_management("in_memory/Intersection")

JohnMorgan · ‎06-27-2012

Hello Marcin,

This is awsome feedback! Now, I really want to get a hold of 10.1. We are currently on 10 but I will inquire today about what it takes to get on 10.1. Do you know if it is easy to upgrade, meaning is it just a download, or do I have to order a DVD?

For now I have adjusted my script to use the in_memory/Intersect_analysis method. However, something still seems a bit off for my script. With the newly adjusted script (below) it runs without error. However, it seemingly just hangs there in the first loop pass. Based on my debug messages it seems as though it might be hanging on Intersect_analysis. Am I doing something wrong with the way that I am calling intersect?

import arcpy

#arcpy.env.workspace = "C:/path/PADUS/Polys.gdb"
arcpy.env.workspace = "C:/path/Polys.gdb"
countiesFC = "Counties"
countyRows = arcpy.SearchCursor(countiesFC) 
countyFields = arcpy.ListFields(countiesFC)

PADUSFC = "PADUS"
PADUSRows = arcpy.SearchCursor(PADUSFC) 
PADUSFields = arcpy.ListFields(PADUSFC)

countyGeom = arcpy.Geometry()
countyGeomList = arcpy.CopyFeatures_management("Counties", countyGeom)
PADUSGeom = arcpy.Geometry()
PADUSGeomList = arcpy.CopyFeatures_management("PADUS", PADUSGeom)

countyIndex = 0
#arcpy.AddMessage(str(len(countyGeomList)))
for countyRow in countyRows:
    countyGeom = countyGeomList[countyIndex]
    arcpy.AddMessage("county: " + countyRow.name + " - " + str(countyGeom.area))
    PADUSFC = "PADUS"
    PADUSRows = arcpy.SearchCursor(PADUSFC) 
    PADUSIndex = 0
    for PADUSRow in PADUSRows:
        if arcpy.Exists("in_memory/Intersection"):
            arcpy.Delete_management("in_memory/Intersection")
            arcpy("del in_mem")
        PADUSGeom = PADUSGeomList[PADUSIndex]
        if PADUSGeom.overlaps(countyGeom):
  arcpy.AddMessage("PADUS: " + PADUSRow.own_type + "overlaps " + countyRow.NAME)
        #here want to find the area of overlap between the two geometry polygons
  arcpy.Intersect_analysis([countyGeom, PADUSGeom], "in_memory/Intersection")
  intersectionRows = arcpy.SearchCursor("in_memory/Intersection")
  for intersectionRow in intersectionRows:
   intersectionArea = intersectionRow.Shape.area
   arcpy.AddMessage("Overlap Area: " + str(intersectionArea))
    PADUSIndex = PADUSIndex + 1
    del PADUSFC
    del PADUSRow
    del PADUSRows
    #arcpy.AddMessage(str(countyIndex))
    countyIndex = countyIndex + 1

MarcinGasior · ‎06-27-2012

I didn't upgraded to 10.1. But the web help system for 10.1 is already online.

For your code: you dont need to create additional list of geometries - you have an access directly through SHAPE fiels in search cursor.
Here's my full testing script (with extended error cathing part) which is working. When arcpy.env.overwriteOutput = True you don't need to care about in_memory feature class - it will be overwritten every time.

def main():
    try:
        import arcpy, sys, traceback, os
        arcpy.env.overwriteOutput = True
        countiesFC = r"C:\tmp\Test.gdb\Counties"
        countyShapeName = arcpy.Describe(countiesFC).ShapeFieldName
        countyRows = arcpy.SearchCursor(countiesFC)

        polyFC = r"C:\tmp\t.gdb\Polygons"
        polyShapeName = arcpy.Describe(polyFC).ShapeFieldName

        for countyRow in countyRows:
            countyGeom = countyRow.getValue(countyShapeName)

            polyRows = arcpy.SearchCursor(polyFC)
            for polyRow in polyRows:
                polyGeom = polyRow.getValue(polyShapeName)

                if countyGeom.contains(polyGeom):
                    arcpy.Intersect_analysis([countyGeom, polyGeom], "in_memory/Intersection")
                    intersectionRows = arcpy.SearchCursor("in_memory/Intersection")
                    for intersectionRow in intersectionRows:
                        intersectionArea = intersectionRow.Shape.area
                        print "Overlap Area: " + str(intersectionArea)

        del countyRows, countyRow, polyRows

    except:
        print arcpy.GetMessages()
        # Get the traceback object
        tb = sys.exc_info()[2]
        tbinfo = traceback.format_tb(tb)[0]

        # Concatenate information together concerning the error into a
        #   message string
        pymsg = tbinfo + "\n" + str(sys.exc_type)+ ": " + str(sys.exc_value)

        # Return python error messages for use with a script tool
        arcpy.AddError(pymsg)

        # Print Python error messages for use in Python/PythonWin
        print pymsg

if __name__ == '__main__':
    main()

JohnMorgan · ‎06-27-2012

Hello Marcin,

You sir are a scholar and a gentleman! Really though I appreciate your help and answers to my question. And I am quite impressed that the comparison of polygon by area in this way can be done with ArcPy.

Cheers,
Derek

DanLee · ‎06-27-2012

While you guys are going through the trouble fixing the code, I wonder if you are aware of the Intersect tool which gives you the desired result:
http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00080000000p000000

You can call Intersect tool in arcpy; once you obtain the overlapping areas in the output of Intersect tool, you can do further analysis. Have a look of the examples in the help page above. Does that simply your workflow and coding?

JohnMorgan · ‎06-27-2012

Hi Dan,

Good point. I had originally looked at the intersect tool. I would have preferred to use that method as it does everything out-of-the-box for me so to speak. However, the polygon data-set I am working with has a huge amount of records (polys) and I need to do a national-scale analysis. As I remember the intersect tool didn't seem to work so well with this huge national-scale amount of data. At least this was the case on my sluggish laptop. So, I thought it would be nice to be able to script it. I realize I also could have clipped out zonal regions to probably make smaller data-sets. Also, a coworker had suggested converting the polygons to raster and doing zonal statistics. All great ideas, but I am happy to learn how to do this with python.

Thanks,
Derek