A common task is to calculate the percentage overlap between a polygon feature class overlaying another polygon feature class. One general workflow is to use the (Pairwise) Intersect tool with the Join Attributes parameter set to Only feature IDs, add a field to store the percentage overlap attributes, join the underlay polygons to the output based on Feature IDs, and then use Field Calculator and divide the Shape_Area for the Intersect output by the Shape_Area from the original underlay polygon. Simple stuff really but manually time consuming if you need to perform routinely. Let’s enhance using ArcPy and the Pairwise Intersect tool for ArcGIS Pro to achieve the output from the workflow described above. You can find more about the Pairwise Intersect tool in the Esri documentation here. ArcGIS Pro 3.1.0 was used when implementing the following.
Interested in learning ArcPy? check out our course.
Let’s get the show on the road by importing the ArcPy module.
import arcpy
We’re going to need two required user input parameters; in_features (Data Type: Value Table), the Pairwise Intersect tool only allows to Feature Layers as input and we will restrict these to polygons in our tool’s GUI, and out_feature_class (Data Type: Feature Class) which will be the resulting polygon feature class from our workflow.
## only two input features allowed
in_features = arcpy.GetParameterAsText(0)
## output feature class
out_feature_class = arcpy.GetParameterAsText(1)
Next, we’ll set some required objects to aid with the workflow. The code is commented, read through and get familiar with each statement.
## we only want the FIDs in the intersect attribute table returned
join_attributes = "ONLY_FID"
## output type is the same as input - Polygon, we will limit this is the
## tool GUI
output_type="INPUT"
## the field name with the FIDs of the underlay polygons
## when you run the intersct tool this is the name of the field created with
## the original FIDs of the 1st feature class entered into the in_features
## eg FID_Name_Of_Feature_Class
fid_fld_name_1 = "FID_{0}".format(arcpy.Describe(in_features.split(";")[0].split("\\")[-1]).baseName)
## the field name with the FIDs of the underlay polygons
## when you run the intersct tool this is the name of the field created with
## the original FIDs of the 2nd feature class entered into the in_features
## eg FID_Name_Of_Feature_Class
fid_fld_name_2 = "FID_{0}".format(arcpy.Describe(in_features.split(";")[1].split("\\")[-1]).baseName)
## the name of the 1st table input
in_tbl = in_features.split(";")[0]
## the name of the 2nd table input
join_tbl = in_features.split(";")[1]
## the OID field name for the join_tbl
join_fld = [fld.name for fld in arcpy.ListFields(join_tbl) if fld.type == "OID"][0]
## these lists will be used in the search and update cursors
search_flds = [join_fld, "SHAPE@AREA"]
update_flds = [fid_fld_name_1, fid_fld_name_2, "overlap_percent"]
Now we call the big dawg! The Pairwise Intersect tool and we place the output feature class into the memory workspace.
## use the pairwise intersect tool, this will return a polygon feature class
## that contains the interected portions only.
memory_fc = arcpy.analysis.PairwiseIntersect(in_features, "memory\\pwise_fc", join_attributes, output_type=output_type)
Time to work some magic and add a field to store the percentage overlap values and then populate the attributes. Instead of using ArcPy to join tables together, we will use a Python dictionary like a look up table.
## add a field to contain our percentage of overlap
arcpy.management.AddField(memory_fc, "overlap_percent", "DOUBLE")
## key (fid_fld_name_1, fid_fld_name_2), value Shape_Area
fid_dict = {(row[0],row[1]):row[2] for row in arcpy.da.SearchCursor(memory_fc, [fid_fld_name_1, fid_fld_name_2, "SHAPE@AREA"])}
## SQL expression for where clause when getting necessary features from the
## first of the in_features
sql_exp = "{0} IN ({1})".format(join_fld, ",".join(str(x[1]) for x in fid_dict.keys()))
## iterate through each feature from the sql_exp for the 2nd in_features
with arcpy.da.SearchCursor(join_tbl, search_flds, sql_exp) as cursor:
## for each row (polygon) in the feature class
for row in cursor:
## iterate over the memory_fc records where the fid_fld_name_2 is equal
## to the FID in the current search cursor.
## in the memory fc
with arcpy.da.UpdateCursor(memory_fc, update_flds, "{0} = {1}".format(fid_fld_name_2, row[0])) as u_cursor:
## for each record
for u_row in u_cursor:
## set the percentage_overlap to be the overlap area / the original
## area for the overlapped polygoon from the 2nd fc from in_features
u_row[2] = fid_dict[(u_row[0], u_row[1])]/row[1]
u_cursor.updateRow(u_row)
Our percentage overlaps have been calculated so now we can write the output to disk. The output will simply contain the FIDs of the overlapping polygons and the overlap_percent field that we added. Before writing to disk you could join the tables based on the FIDs for one or both of the original polygon feature classes and also have these in the output, but that’s for another day’s work.
## write output to disk
arcpy.conversion.ExportFeatures(memory_fc, out_feature_class)
And as always, we clean up the memory workspace.
## clean-up memory workspace by deleting the memory_fc
arcpy.management.Delete(memory_fc)
That didn’t take much code now did it? There’s probably more comments in there than actual lines of working code. Save your script and let’s head over to ArcGIS Pro to create our custom tool for use.
Right-click on your toolbox/toolset of choice and select New > Script. The New Script window will appear. In the General tab set Name to polygonOverlapPercentage, Label to Polygon-Polygon Overlap Percentage, and the Description to Polygon-Polygon Overlap Percentage with Pairwise Intersect and ArcPy.
In the Parameters tab set as per below. Set the Input Features with a Data Type of Feature Layer and select Multiple Values, use the Filter to restrict Feature Type to Polygon. For Output Feature Class, set the Data Type to Feature Class, the Direction to Output and the Filter to Feature Type Polygon.
In the Execution tab, click the folder icon in the top-right corner and ad your saved Python script.
In the Validation tab update the updateMessage function as per below to restrict the input to two feature classes only. Click OK and take the tool for a spin.
def updateMessages(self):
# Customize messages for the parameters.
# This gets called after standard validation.
in_count = len(self.params[0].valueAsText.split(";"))
if in_count > 2:
self.params[0].setErrorMessage("Can only accept two input polygon feature classes")
else:
self.params[0].clearMessage()
return
You can download the tool and other custom tools over on this page.
Here’s all the code above for convenience.
import arcpy
################################################################################
## Esri Documentation
## https://pro.arcgis.com/en/pro-app/latest/tool-reference/analysis/pairwise-intersect.htm
## https://pro.arcgis.com/en/pro-app/latest/arcpy/functions/getparameterastext.htm
## https://pro.arcgis.com/en/pro-app/latest/arcpy/functions/listfields.htm
## https://pro.arcgis.com/en/pro-app/latest/tool-reference/data-management/add-field.htm
## https://pro.arcgis.com/en/pro-app/latest/arcpy/data-access/searchcursor-class.htm
## https://pro.arcgis.com/en/pro-app/latest/arcpy/data-access/updatecursor-class.htm
## https://pro.arcgis.com/en/pro-app/latest/tool-reference/conversion/export-features.htm
## https://pro.arcgis.com/en/pro-app/latest/tool-reference/data-management/delete.htm
## https://pro.arcgis.com/en/pro-app/latest/help/analysis/geoprocessing/basics/the-in-memory-workspace.htm
##
## ArcGIS Pro Version: 3.1.0
##
################################################################################
################################################################################
## USER INPUTS
## only two input features allowed
in_features = arcpy.GetParameterAsText(0)
## output feature class
out_feature_class = arcpy.GetParameterAsText(1)
################################################################################
## REQUIRED OBJECTS
## we only want the FIDs in the intersect attribute table returned
join_attributes = "ONLY_FID"
## output type is the same as input - Polygon, we will limit this is the
## tool GUI
output_type="INPUT"
## the field name with the FIDs of the underlay polygons
## when you run the intersct tool this is the name of the field created with
## the original FIDs of the 1st feature class entered into the in_features
## eg FID_Name_Of_Feature_Class
fid_fld_name_1 = "FID_{0}".format(arcpy.Describe(in_features.split(";")[0].split("\\")[-1]).baseName)
## the field name with the FIDs of the underlay polygons
## when you run the intersct tool this is the name of the field created with
## the original FIDs of the 2nd feature class entered into the in_features
## eg FID_Name_Of_Feature_Class
fid_fld_name_2 = "FID_{0}".format(arcpy.Describe(in_features.split(";")[1].split("\\")[-1]).baseName)
## the name of the 1st table input
in_tbl = in_features.split(";")[0]
## the name of the 2nd table input
join_tbl = in_features.split(";")[1]
## the OID field name for the join_tbl
join_fld = [fld.name for fld in arcpy.ListFields(join_tbl) if fld.type == "OID"][0]
## these lists will be used in the search and update cursors
search_flds = [join_fld, "SHAPE@AREA"]
update_flds = [fid_fld_name_1, fid_fld_name_2, "overlap_percent"]
################################################################################
## PAIRWISE INTERSECT TO GET OVERLAPPING POLYGONS
## use the pairwise intersect tool, this will return a polygon feature class
## that contains the interected portions only.
memory_fc = arcpy.analysis.PairwiseIntersect(in_features, "memory\\pwise_fc", join_attributes, output_type=output_type)
################################################################################
## CALCULATE THE PERCENTAGE OF OVERLAP
## add a field to contain our percentage of overlap
arcpy.management.AddField(memory_fc, "overlap_percent", "DOUBLE")
## key (fid_fld_name_1, fid_fld_name_2), value Shape_Area
fid_dict = {(row[0],row[1]):row[2] for row in arcpy.da.SearchCursor(memory_fc, [fid_fld_name_1, fid_fld_name_2, "SHAPE@AREA"])}
## SQL expression for where clause when getting necessary features from the
## first of the in_features
sql_exp = "{0} IN ({1})".format(join_fld, ",".join(str(x[1]) for x in fid_dict.keys()))
## iterate through each feature from the sql_exp for the 2nd in_features
with arcpy.da.SearchCursor(join_tbl, search_flds, sql_exp) as cursor:
## for each row (polygon) in the feature class
for row in cursor:
## iterate over the memory_fc records where the fid_fld_name_2 is equal
## to the FID in the current search cursor.
with arcpy.da.UpdateCursor(memory_fc, update_flds, "{0} = {1}".format(fid_fld_name_2, row[0])) as u_cursor:
## for each record
for u_row in u_cursor:
## set the percentage_overlap to be the overlap area / the original
## area for the overlapped polygoon from the 2nd fc from in_features
u_row[2] = fid_dict[(u_row[0], u_row[1])]/row[1]
u_cursor.updateRow(u_row)
################################################################################
## CALCULATE THE PERCENTAGE OF OVERLAP
## write output to disk
arcpy.conversion.ExportFeatures(memory_fc, out_feature_class)
################################################################################
## CLEAN-UP MEMORY WORKSPACE
## clean-up memory workspace by deleteing the memory_fc
arcpy.management.Delete(memory_fc)
################################################################################