Mesuring Percentage of Feature Polygons that Overlap Each Other in Python

DanielInterrante1 · ‎02-13-2015

Hello,

I'm trying to gather preliminary census data within a 15 minute walkshed and a 15 min bikeshed boundary that we have drawn around a development site. As expected, neither of the "shed" layers fit perfectly into the Census Block Groups (CBGs) we have. Now I know this is far from scientific, but we were thinking for the CBGs that only partially overlap the shed layers, calculate the percentage of the overlap, place that percentage in a field, and multiply statistics by that percentage. I've been working on this script for over a week now and I simply cannot get it to work right. Basically what I want it to do is:

-Create an Update Cursor in the CBGs layer

-Select each CBG polygon as it scrolls through

-Clip the "shed" layer to the selected polygon

-Measure the shape area for the clipped layer

-Divide the area of the clipped layer by the CBG

-Set the result as the value in the percentage field in CBG

-Delete the clipped polygon

-Repeat

In its current form the script is completing without error, but with whole numbers that clearly are not the desired percentage. Every time I try to create an AddMessage so that I can see what values are being divided I then get an error and am still unable to see where the values are coming. This is the first time I've ever built a python script without step-by-step instructions so I'm very lost. If anyone sees the errors in my code, any pointers would be greatly appreciated. I have the script written below and have attached a zipped folder with the .py, the bikeshed, a few selected CBGs, and the toolbox I'm executing it with.

import arcpy

censuslayer = arcpy.GetParameter(0)#census unit layer
percentField = arcpy.GetParameter(1)#field where percentage will be input
studylayer = arcpy.GetParameter(2)#walk or bikeshed layer



#currently unused#arcpy.AddField_management (censuslayer, "perc_cov", "float", "3", "3") #adds blank percentage field to census units
censussize = arcpy.AddField_management (censuslayer, "shapearea", "FLOAT")#adds shape area to census
arcpy.CalculateField_management (censuslayer, "shapearea", "!shape.area@squarefeet!", "PYTHON") #calculates the geometry


#for polygon in studylayer:
arcpy.AddField_management (studylayer, "shapearea", "FLOAT") #adds a shape area to walk or bikeshed

  
#search = arcpy.SearchCursor(censuslayer)
update = arcpy.UpdateCursor(censuslayer)#creates update cursor to move through census layers

for row in update:
    ceblocksize=float(row.getValue("shapearea"))#places shape area value in variable
    currentCen=int(row.getValue("OBJECTID"))#places object ID field name in variable
    holder=("in_memory\holder")#output
    whereClause=str("'OBJECTID' = " + "'"+ str(currentCen)+"'")#formula to indicate record to select
    selectBlock=arcpy.Select_analysis(censuslayer, holder, whereClause)#select record where cursor should be
    
    temppolygon = arcpy.Clip_analysis (studylayer, selectBlock, "templayer") #clip shape from shed layer to selected feature
    
    studysize = arcpy.CalculateField_management (temppolygon, "shapearea", "!shape.area@squarefeet!", "PYTHON") #calculate geometry  of clipped feature
    studysearch = arcpy.SearchCursor(temppolygon)# new cursor in clipped feature
    for only in studysearch:
        shedsize = only.getValue("shapearea")#store the shape area in a variable
    
        
    #arcpy.AddMessage(shedsize"/"ceblocksize)
    percent = float(shedsize)/ceblocksize #get the value to place in new field
            
        
    row.setValue(percentField, percent)#sets the value in field
    update.updateRow(row)

    arcpy.Delete_management(temppolygon)#deletes the clipped polygon
    row = update.next

Edit: Formatted Code Block

BlakeTerhune · ‎02-13-2015

I'll just leave this here while someone in the community works on posting an answer to your question

Posting Code blocks in the new GeoNet

DanielInterrante1 · ‎02-13-2015

Thank you!

JoshuaBixby · ‎02-18-2015

Try the following, see if this is what you are after:

import arcpy

census =        #census unit feature class  
percentField=   #field where percentage will be input  
study =         #walk or bikeshed feature class  

#make feature layer for using with SelectLayerByLocation
arcpy.MakeFeatureLayer_management(census, "census")

#make search cursor over study areas, even if just one.
studyCursor = arcpy.da.SearchCursor(study, ["OID@", "SHAPE@"])
for oid, s_shape in studyCursor:
    #for given study area, select census blocks that intersect
    arcpy.SelectLayerByLocation_management("census", "INTERSECT", s_shape)

    #make update cursor for census blocks selected above
    censusCursor = arcpy.da.UpdateCursor("census", ["SHAPE@", percentField])
    for c_shape, _ in censusCursor:
        #get area of census block
        area = c_shape.area

        #find pct of census block that intersects study area
        pct = c_shape.intersect(s_shape, 4).area / area * 100

        #update record
        censusCursor.updateRow([c_shape, pct])
    del censusCursor
del studyCursor

The code above is designed to be run as a standalone script, not as a Python tool. If the functional code works for you, it can be moved into a Python tool without much effort.

PeterShoemark · ‎05-01-2016

Hi, I just came across this post and have a similar requirement.

I copied the code above and added in a line to add the percentField field, (works fine) however get the following eror message regarding the pct line.

PYTHON ERRORS:

Traceback info:

File "C:\Data\python\Percenatge_overlap.py", line 22, in <module>

pct = c_shape.intersect(s_shape, 4).area / area * 100

Error Info:

ArcPy ERRORS:

fwiw I'm on Win 7 x64, arcgis 10.2.1, Python 2.7.5

Any thoughts?

Cheers, Peter

DanPatterson_Retired · ‎05-01-2016

# actually a check to see if the shape is not null first
area = c_shape.area 
 #find pct of census block that intersects study area 
if not(c_shape.disjoint(s_shape)):

    pct = c_shape.intersect(s_shape, 4).area / area * 100
else:
    area = 0.0
#update record 
 censusCursor.updateRow([c_shape, pct])

whether that is the cause of the error or not.... a check for 0 for the overlap should be made.... think about it... what if the polygons don't overlap? you shouldn't be dividing by 0. So pay attention to c_shape and s_shape to ensure that they exist and are not null (they should be ok), bt the intersect may be the issue

PeterShoemark · ‎05-01-2016

Thanks, I should have thought to add the checker and have now done so.

Unfortunately I still get the same error message.

Note I failed to mention the data sets are in a gdb just in case that makes a difference.

I can confirm both polygon data sets exist and and a 'select where completely within' returns the following result.

Study count = 4230

Census count = 374

(A min bounding rec is nominally the same for both sets of polygons)

Study polygons completely within Census = 2186

I believe that would confirm that there are 2044 that have some degree of overlap of the Census polygons.

Does that sound correct?

JoshuaBixby · ‎05-03-2016

Is there anything after "ArcPy ERRORS?"

Although I don't think Esri's geometry validation tools are that robust, they do occasionally find issues. Where the error is occurring makes me think you may have an invalid geometry in one of the datasets. If you run Check Geometry on the datasets, does it return anything?

Another thought is projections. Do both datasets have projections? Are they the same projection? What if you reproject both datasets into the same projection, does the error go away?

One final thought, any ZM polygons? I have occasionally had issues with ArcPy Geometry classes and Z, M, or ZM polygons.

PeterShoemark · ‎05-04-2016

It looks like bad geometry (self overlaps) as I did manage to run it cleanly on a small subset. I will do a clean up and try again but maybe a day or so yet.

NeilAyres · ‎05-01-2016

When doing something similar I found it necessary to check disjoint & touching...

# if not disjoint and not touching, must intersect
if not geom1.disjoint(geom2) and not geom1.touches(geom2):