2

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" 
density_ras = r"D:\Downloads\block_clipped_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_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 
final_table = r"D:\Downloads\2013_A30min_weekday.dbf" 
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.

  • 1
    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. – Chris W Apr 7 '15 at 23:08
  • 1
    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 Apr 7 '15 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 Apr 8 '15 at 1:21
  • Just another comment on this, when doing this type of work, be aware of the Modifiable Areal Unit Problem (MAUP). – Fezter Apr 9 '15 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. – Chris W Apr 11 '15 at 0:30
1

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:
        arcpy.AddMessage(row[0])
        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.AddMessage("Defining neighbours...")
    arcpy.SpatialJoin_analysis(store_buffs, store_buffs, joinLR,"JOIN_ONE_TO_MANY")
    arcpy.AddMessage("Creating empty dictionary")
    dictFeatures = {}
    with arcpy.da.SearchCursor(store_buffs, parID) as cursor:
        for row in cursor:
            dictFeatures[row[0]]=()
        del row, cursor
    arcpy.AddMessage("Assigning neighbours...")
    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    
    arcpy.AddMessage("Defining non-overlapping subsets...")
    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.AddMessage("Running Statistics...")
        arcpy.gp.ZonalStatisticsAsTable_sa(store_buffs, parID, density_ras, dbf, "DATA", "SUM")
        arcpy.AddMessage("Data transfer...")
        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 Apr 8 '15 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. – FelixIP Apr 8 '15 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 – user132058 Nov 18 '18 at 8:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.