Calculating the area of intersect between features in different layers?

[MarkCooper1]

Hi

Does anyone know if there is a python script (or anything similar) that would enable me to select a layer and calculate the area of intersection between its features and those in another layer. For example I have a layer containing potential development sites and I want to find out the area of each site that intersects a constraint (eg flood plain) and append this to the sites layer. The sites layer has about 300 sites and there are about 80 constraint layers.

At the moment I am using an old script (http://arcscripts.esri.com/details.asp?dbid=10579) but I don't think this will work when we start using v10.

Thanks
Mark

[JakeSkinner]

You could accomplish this by intersecting the sites and floodplain layer, perform a join and then calculate the percent coverage using the 'Shape_Area' fields.  Here is an example:

import arcpy
from arcpy import env
env.workspace = r"C:\temp\python\test.gdb"
env.overwriteOutput = True

fc1 = "DevelopmentSites"
fc2 = "FloodPlain"

# Add field to DevelopmentSites
arcpy.AddField_management(fc1, "Percentage", "Double")

# Perform an intersect between DevelopmentSites & FloodPlain feature classes
arcpy.Intersect_analysis([[fc1], [fc2]], "Intersect")

# Join intersecting feature class with DevelopmentSites
arcpy.MakeFeatureLayer_management(fc1, "development_layer")
arcpy.AddJoin_management("development_layer", "OBJECTID", "Intersect", "FID_DevelopmentSites")

# Calculate Percentage of flood plain coverage
arcpy.CalculateField_management("development_layer", "DevelopmentSites.Percentage", "!Intersect.Shape_Area! / !DevelopmentSites.Shape_Area! * 100", "PYTHON")

# Delete Intersect feature class
arcpy.Delete_management("Intersect")

The result will be a new field called 'Percentage' within the DevelopmentSites feature class with the percent coverage of the flood plain feature class.

[SimchaLevental]

Would there be a way for me to calculate the percentage of categories of flood plain type rather that just the percent intersect?

[T__WayneWhitley]

I have done something similar recently, based on clip output (see the 3rd line), but this will work as well on intersect output - here is the snippet showing the part you're interested in:

import arcpy, time
from operator import itemgetter

arcpy.Clip_analysis('theData', 'buffer500', 'clip500')
 
# define dictionary
summary = dict()
 
rows = arcpy.SearchCursor('clip500')

# type or category is in my feature class called DESCRIPTION.
for row in rows:
        if row.DESCRIPTION in summary:
            summary[row.DESCRIPTION] = summary[row.DESCRIPTION] + row.shape_area
        else:
            summary[row.DESCRIPTION] = row.shape_area

# dictionary now loaded, from here on out are simply acreage conversions, sorting.
sumTotalAc = sum(summary.values())/43560.0

arcpy.AddMessage('\nTotal area: ' + str(sumTotalAc))

sortSummaryAc = sorted(([a, b/43560.0] for a, b in summary.iteritems()), key = itemgetter(1), reverse = True)

arcpy.AddMessage(sortSummaryAc)

EDIT:
ahh, you said you are after percentages, which if you think about it is kind of incidental once you have total acreage.  I have the sum(summary.values()) converted to acres, sumTotalAc, above.  I guess I could have computed percentages in the dictionary, but I just did it on-the-fly since I needed individual DESCRIPTION acreages as well as percentage of the total in some text elements displayed on the map ---- so if this makes sense to you, here's that part (the whole script is really too long to make my point clear):

for each in sortSummaryAc:
        testgrab = elmhndl.pop()
        if round(each[1], 1) > 0.0:
            # here's the line adding both the acreage and percentage
            testgrab.text = each[0] + ' - ' + str(round(each[1], 1)) + '  acres   (' + str(round((each[1]/sumTotalAc*100), 1)) + ' %)'    
            index = assign(testgrab, [initPosX, initPosY])
            initPosY = initPosY - float(0.2122)
        else:
            arcpy.AddMessage('\n...' + each[0] + ' is too small, removing element to remaining set...')
            testgrab.text = '-'
            index = assign(testgrab, [initPosX1, initPosY1])

Hope that's clear.

Enjoy,
Wayne