# Using Zonal Statistics As Table for overlapping polygons in ArcPy?

I'm using ArcGIS 10.2 for Desktop and would like to calculate the population within a specific distance from bus stops.

My input layers are:

1. census block (polygon) with population
2. undissolved bus buffers that overlap each other (polygon, in 1 layer)

I calculated population density and converted my block data into raster, then used Zonal Statistics as Table tool to calculate the sum of population in each buffer. If I use this tool directly, it will ignore 2/3 of my buffers because they overlap with others. Therefore I need code to iterate this through each buffer (about 600 in total).

My situation is the same to the one in this topic Intersecting Overlapping Polygons Individually in ArcGIS

I used the code by @phloem but it returns a table with only one row (which is one of the buffers I have). Could someone look through and give me some advice? I'm new to Python.

``````store_buffs = r"H:\DMLA\NTAS\Map\Shapefiles\2013_A30min_weekday.shp"
table_list = []
with arcpy.da.SearchCursor(store_buffs, ["FacilityID"]) as cursor:
for row in cursor:
exp = '"FacilityID" = ' + str(row[0])
temp_table10 = r"in_memory\temp_table10" + str(row[0])
temp_shp10 = r'in_memory\temp_shp10'
arcpy.Select_analysis(store_buffs, temp_shp10, exp)
arcpy.sa.ZonalStatisticsAsTable(temp_shp10, 'FacilityID', density_ras, temp_table10, 'DATA', 'SUM')
table_list.append(temp_table10)
del row
arcpy.Merge_management(table_list, final_table)
``````

I cannot post comment on phloem's answer because I don't have enough reputation. That's why I have to create another topic.

I found an ESRI supplemental tool called Zonal Statistics As Table 2, which allows me to calculate overlapping polygons. However, this tool takes forever to run (6 hours for me). I would like to have a more efficient way to process data.

• Do you specifically want to stay in arcpy and using rasters? Even if you want to script/model it, in my opinion there are much easier ways to do this by staying in vector and just using your two original input files. Otherwise @phloem might be able to help you troubleshoot the code. Commented Apr 7, 2015 at 23:08
• Agreed, I don't see how using rasters would help here. You could simply use the intersect tool. resources.arcgis.com/en/help/main/10.2/index.html#//… and the summary statistics tool resources.arcgis.com/en/help/main/10.2/index.html#//…
– Fezter
Commented Apr 7, 2015 at 23:18
• Hi Chris W and Fezter, I'm fine with any solution, no need to stick with raster. The thing is my buffer layer contains 600 overlapping polygons (many of them has exactly the same shape and location), so I couldn't use Intersect or the regular statistics table. @phloem could you please help me with the code? Thanks!
– GGT
Commented Apr 8, 2015 at 1:21
• Just another comment on this, when doing this type of work, be aware of the Modifiable Areal Unit Problem (MAUP).
– Fezter
Commented Apr 9, 2015 at 0:44
• Tabulate Statistics might do what you want in one go. Or, make feature layer out of census blocks and enable ratio policy on the pop field. Use Union on the two layers instead of Intersect. This should preserve ID of overlapping buffers, albeit in separate columns, and with the ratio policy the pop should be apportioned by area to the new polygons. A little cleanup and you can use Summary Statistics with BufferID as a case field and the new pop field as a sum statistic to get total apportioned pop per buffer. Commented Apr 11, 2015 at 0:30

You CAN do it using vectors, 100% agree with others. Anyway, I slightly modified your code, replacing very long names and it works as expected:

``````import arcpy
store_buffs = r"D:\Scratch\A30min_weekday.shp"
density_ras = r"D:\Scratch\Raster1.tif"
table_list = []
with arcpy.da.SearchCursor(store_buffs, ["FacilityID"]) as cursor:
for row in cursor:
exp = '"FacilityID" = ' + str(row[0])
temp_table10 = r"in_memory\temp_" + str(row[0])
temp_shp10 = r'in_memory\temp_shp10'
arcpy.Select_analysis(store_buffs, temp_shp10, exp)
arcpy.sa.ZonalStatisticsAsTable(temp_shp10, 'FacilityID', density_ras, temp_table10, 'DATA', 'SUM')
table_list.append(temp_table10)
del row
final_table = r"D:\scratch\all_rows.dbf"
arcpy.Merge_management(table_list, final_table)
``````

It takes ages, i.e. almost 4 mins for 100 pgons. Try this one, it takes 25 sec to do the same job. It works from ArcMap, assumes you have layer called 'A30min_weekday' and it has field SUM.

``````import arcpy, traceback, os, sys
from arcpy import env
env.overwriteOutput = True
mxd = arcpy.mapping.MapDocument("CURRENT")
store_buffs=arcpy.mapping.ListLayers(mxd,'A30min_weekday')[0]
density_ras = r"D:\Scratch\Raster1.tif"
parID="FacilityID"
parID2="FacilityID_1"
env.workspace = "in_memory"
dbf="stat"
joinLR="SD"

try:
def showPyMessage():
arcpy.AddMessage(str(time.ctime()) + " - " + message)
def Get_V(aKey):
try:
return smallDict[aKey]
except:
return (-1)
arcpy.SpatialJoin_analysis(store_buffs, store_buffs, joinLR,"JOIN_ONE_TO_MANY")
dictFeatures = {}
with arcpy.da.SearchCursor(store_buffs, parID) as cursor:
for row in cursor:
dictFeatures[row[0]]=()
del row, cursor
with arcpy.da.SearchCursor(joinLR, (parID,parID2)) as cursor:
for row in cursor:
aKey=row[0]
aList=dictFeatures[aKey]
aList+=(row[1],)
dictFeatures[aKey]=aList
del row, cursor
runNo=0
while (True):
toShow,toHide=(),()
nF=len(dictFeatures)
for item in dictFeatures:
if item not in toShow and item not in toHide:
toShow+=(item,)
toHide+=(dictFeatures[item])
m=len(toShow)
quer='"FacilityID" IN '+str(toShow)
if m==1:
quer='"FacilityID" = '+str(toShow[0])
store_buffs.definitionQuery=quer
runNo+=1
arcpy.AddMessage("Run %i, %i polygon(s) found" % (runNo,m))
arcpy.gp.ZonalStatisticsAsTable_sa(store_buffs, parID, density_ras, dbf, "DATA", "SUM")
smallDict={}
with arcpy.da.SearchCursor(dbf, (parID,"SUM")) as cursor:
for row in cursor:
smallDict[row[0]]=row[1]
del row, cursor
with arcpy.da.UpdateCursor(store_buffs, (parID,"SUM")) as cursor:
for row in cursor:
aKey=row[0]
row[1]=Get_V(aKey)
cursor.updateRow(row)
del row, cursor
for item in toShow:
del dictFeatures[item]
m=len(dictFeatures)
if m==0:
break
store_buffs.definitionQuery=""
except:
message = "\n*** PYTHON ERRORS *** "; showPyMessage()
message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()
``````
• Hi FelixIP, thanks so much for the code! I've been running the first one from last night and it has not finished @_@. The second one returns this message: Runtime error Traceback (most recent call last): File "<string>", line 74, in <module> NameError: name 'traceback' is not defined Could you please share the way to do it in vector, too? It would be extremely helpful because I don't want to spend 3-4 hours to run each buffer layer (600 polygons for each layer).
– GGT
Commented Apr 8, 2015 at 14:16
• 2nd script assumes that Facilityid is integer. If this is not a case, create field ABC integer and populate it using field calculator, e.g. FID+1. Inside script replace Faciltyid with ABC. Commented Apr 8, 2015 at 20:14
• you may not have the traceback module,put it in "D:\Python27\ArcGIS10.2\Lib\site-packages"(it's mine),download traceback module and install Commented Nov 18, 2018 at 8:54