nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Zonal Statistics as Table: Overlapping Polygons: Solution Required

PeterWilson

I'm currently using ArcGIS 10 SP2 on Windows 7 64bit. I'm working on a very large railway alignment design project and need to generate "Average Catchment Slopes" for each watershed. The problem that I have is that the "Zonal Statistics as Table" tool does not process overlapping polygons. It's mentions in the help that if the polygons overlap, that should be processed interactively to overcome this. I haven't been able to find any supporting documentation to achieve this. Just a note, I have over 600 watersheds that I need to process and a manual process will not be possible if I'm to meet my deadline. Any advice will be appreciated.

Regards

Vind meer berichten getagd met

Accepted answers

All comments

JamieDrisdelle

Take a look at the attached tool. It will allow you to do zonal statistics as table on overlapping polygons. you can add this tool directly to ArcToolbox.

Jamie

StatisticsForOverlappingZones.zip

curtvprice

We recently published a toolbox (first prize at the ESRI UC 2011 AppFair!) that includes tools that do overlap processing for statistics and area tabulations, among other things.

USGS Open-File Report 2010�??1268
National Water-Quality Assessment (NAWQA) Area-Characterization Toolbox
By Curtis V. Price, Naomi Nakagaki, and Kerie J. Hitt

http://pubs.usgs.gov/of/2010/1268/

The published version has only limited support for 10.x, but I have an updated version that will be posted there soon. If you are interested, let me know and I can send it your way.

curtvprice

Hey Jamie - what's this (from your script tool):

import gc
# ...
#enable garbage collection
gc.enable()

Does this help with looping geoprocessing tools that seem to eat up memory and crash your session?

PeterWilson

Hi Jamie

Thanks so much, your python script has solved my problem. I've made some minor adjustments to your code to tie in with ArcHydro. Thanks again, its greatly appreciated.

Regards

PeterWilson

Hi Curtis

I would really be interested in looking at the toolbox that you have developed. It would be greatly appreciated if you would be willing to email me a copy for ArcGIS 10. Thanks fo offering your assistance.

Regards

JamieDrisdelle

Hi Curtis,

It is the garbage collection module. We were running into some memory issues after several iterations of the tool. The garbage collection module clean up unreferenced objects in memory. Take a look t this link.

http://docs.python.org/library/gc.html

curtvprice

It would be greatly appreciated if you would be willing to email me a copy for ArcGIS 10.

Peter - please email me (contact info in in my sig below)

PeterWilson

Peter - please email me (contact info in in my sig below)

Hi Curtis

My contact details are : Peter.Wilson@aurecongroup.com, Thanks

Regards

DanielSheehan

Jamie, will this run a lot longer than the regular Zonal Stats (it is right now for me, not sure if its working)?

Also, I got an error at first that I had no FID field, my zonal feature had ObjectID instead of FID fields, so I AddField then CalculatedField FID to be the ObjectID values.

Thanks,
Danny

curtvprice

I got an error at first that I had no FID field, my zonal feature had ObjectID instead of FID fields,...

In Jamie's defense, this was kind of a once-off deal working with shapefiles. The more fancy way to do this (to support tables with different object ID field names) is:

instead of hardcoding:

strFID = '"FID" ='
...
    arcpy.SelectLayerByAttribute_management(zone, newSelection, strFID + str(inRow.FID))

get the tables object id field name:

strFID = arcpy.Describe(inputFeatureZone).OIDFieldName
...
    strQry = '"%s" = %s' % (strFID, inRow.GetValue(strFID))
    arcpy.SelectLayerByAttribute_management(zone, newSelection, strQry)

DanielSheehan

Okay, Thanks!!! I'm batch processing a ton of radial and network buffers. Should I be running this in python rather than in Modelbuilder?

curtvprice

@Daniel -

This is just a lot of iterations and will take quite a while. I doubt if wrapping the Python script inside modelbuilder will buy you any speed, as long as you've checked the box to run the python script "in process". One approach to consider is grouping your input polys into some non-overlapping groups so you can them in parallel with a single run of zonal statistics.

I have a script tool that makes an attempt but it's not good enough for me to feel okay putting it here. If interested in helping me with that project by testing it, contact me off line.

DanielSheehan

Curtis, I think I'm using Jamie's script which I benchmarked it in modelbuilder and would take many days to run (if my computer could even chug along for that long). I was just wondering if it would be significantly faster to run things outside of modelbuilder as a .py file vs. running in modelbuilder?

The first workaround I came up with was to run zonal stats on each Feature split into its own Feature Class (so whereas one buffer distance Feature Class had 1541 records, I now have 1541 Feature Classes). Then I set up Iterator in modelbuilder to create Rasters for each FC (n=1541) (did this b/c I guess Zonal Stats converts vector to raster temporarilty to run) and then run through all the Rasters Zonal Stats individually using iterator. This worked until I ran into a few rasters that didn't produce any stats (why I don't know) and caused the program to crash - I'm curious if iterators in modelbuilder have a theoretical limit.

So now I'm thinking of just converting my Kernel Density surface (the file I want to do Zonal Stats for) to a point file and just using Summary Statistics after an intersect of polygons to points. I'm just exhausted using Zonal Stats for 2 reasons, it keeps crashing using iterators AND 13-character limits for file names (I had to reprocess work to create unique IDs that took a whole day b/c my unique IDs were 15 characters long).

Do you see any reason not to do this vector version of Zonal Stats? And how stable is your tool that you don't feel comfortable putting up on the forum? I may not be of help testing b/c I'm not much of a python programmer (I'm self taught modelbuilder and know very little python and even less VBA).

Danny

Danny

SebastianRoberts

Just found this post. Jamie's script worked great. I am using this for zonal stats in a parcel model so only have a small subset of overlapping tax parcels (condos only). In my case, I do an intersect with just one input (the parcel layer) to extract only the overlapping condo tax parcels. I added this to the beginning of Jamie's script, then process only these parcels with the script. In a second model, I do a selection to subtract these condo tax parcels form the county wide parcel feature layer. I can then do a regular zonal statistics on the feature layer for all the remaining non overlapping polygons. Just thought I would share this as a way to speed up the process if your dataset has a small subset of overlapping polygons.

#create overlapping polygons with intersect tool (intersect with self to get overlapping)
arcpy.Intersect_analysis([inputParcelLayer],FeatureZone,"NO_FID","","")

jillianthonell

Hi Curtis,

It is the garbage collection module. We were running into some memory issues after several iterations of the tool. The garbage collection module clean up unreferenced objects in memory. Take a look t this link.

http://docs.python.org/library/gc.html

Hi Jamie

I was interested to see this thread as I have written a script to handle overlapping polygons for zonalstatistics and have been experiencing similar memory issues with iterations and subsequent crashes. I enabled the garbage collection as in your script, but am still have issues. Can you suggest any further solutions?

Thanks
Jill

SebastianRoberts

Jill,
Did you find any solution? As we have added more overlapping parcels to our feature class (we now have 860), I too seem to be running into memory issues. My script was crashing, so I added another loop that uses a smaller cursor and repeatedly deletes the searchCursor. It no longer crashes, but interestingly enough, a "zone" is processed in about 20 seconds at the outset, and then towards the end of the 800 features, it slows down to about 5 minutes to process one zone. Below is the python code that I am using.

#create overlapping polygons with intersect tool (intersect with self, only creates partial parcels where overlapping)
arcpy.Intersect_analysis([inputParcelLayer],FeatureIntersect,"NO_FID","","")

#make feature layer, join with self intersected parcels, use result as overlapping parcels for zonal stats
arcpy.MakeFeatureLayer_management(inputParcelLayer,"parcelsoverlap")
arcpy.AddJoin_management("parcelsoverlap","APN",FeatureIntersect,"APN","KEEP_COMMON")
arcpy.CopyFeatures_management("parcelsoverlap", FeatureZone)


tempDIR = arcpy.GetSystemEnvironment("TEMP")
#create scratch workspace name
tempWorkspace = "tempWorkspaceForZonalStats"

if not arcpy.Exists(tempDIR + os.sep + tempWorkspace):
    arcpy.CreateFolder_management(tempDIR, tempWorkspace)
arcpy.env.workspace = tempDIR + os.sep + tempWorkspace

#input variables
strFIDname = arcpy.Describe(FeatureZone).OIDFieldName
zone = "zone"

arcpy.env.cellSize = 2 #needed to reduce cell size because zonal stats crashes if the polygon is smaller than raster resolution
gc.enable()  #enable garbage collection
arcpy.CheckOutExtension("Spatial")
cellAssign = "CELL_CENTER"
priority = "NONE"
newSelection = "NEW_SELECTION"

# Make a layer from the feature class
arcpy.MakeFeatureLayer_management(FeatureZone, zone) 
#Loop through features
numFeats = int(arcpy.GetCount_management(FeatureZone).getOutput(0))
groups = []
bracket = 30
while(bracket < numFeats):
 groups.append(bracket)
 bracket = bracket + 30
groups.append(numFeats+1)
for x in groups:
    whileClause = '"OBJECTID" >= ' + str(x-30) + ' and "OBJECTID" < ' + str(x)
    print whileClause
    inRows = arcpy.SearchCursor(FeatureZone,whileClause)
    #inRows = arcpy.SearchCursor(FeatureZone)
    inRow = inRows.next()
    while inRow:
        #selct the fid to process
        strID = inRow.getValue(strFIDname)
        print ' processing objID: ' + str(strID) + 'of ' + str(numFeats) + ' total features.'
        strQry = '"%s" = %s' % (strFIDname, strID)
        arcpy.SelectLayerByAttribute_management(zone, newSelection, strQry)
        #create a unique name for the zone
        uniqueZone = arcpy.CreateUniqueName("zone.shp", arcpy.env.workspace)
        uniqueTable = arcpy.CreateUniqueName("ZSasT", arcpy.env.workspace)
        #create a temporary feature to use for zonal stats as table
        arcpy.CopyFeatures_management(zone, uniqueZone)
        outZSaT = ZonalStatisticsAsTable(uniqueZone, joinField, valueRaster, uniqueTable, "NODATA", "ALL")

        #move to next record.
        inRow = inRows.next()
    del inRow
    del inRows

curtvprice

Interestingly enough, a "zone" is processed in about 20 seconds at the outset, and then towards the end of the 800 features, it slows down to about 5 minutes to process one zone.

My experience is that the gp/arcpy scratch name methods can be very slow when the workspace in which you are generating scratchnames gets full of many files. This can slow down all tools!

You may get better performance if you just create an empty folder/workspace and generate names using sequence numbers. I have gotten in the habit with my general purpose script tools (like the [post=124423] NAWQA toolbox[/post]) of creating a temporary folder and workspace so I don't need to worry about file name collisions and can avoid waiting for the scratchname method!

SebastianRoberts

Hi Curtis,
Thanks for your reply. That improved performance greatly. I just moved my temp directory creation inside the loop so that a new directory is created for every loop iteration (every 30 feature clasess), and now it is on target to finish in about 3 hours rather than 30.

Sebastian

for x in groups:
    #create the tempworkspace
    tempWS = tempWorkspacePrefix + str(tempWorkspaceSuffix)
    if not arcpy.Exists(tempDIR + os.sep + tempWS):
        arcpy.CreateFolder_management(tempDIR, tempWS)
    arcpy.env.workspace = tempDIR + os.sep + tempWS
    while......................(rest of the loop as in my previous post)
tempWorkspaceSuffix = tempWorkspaceSuffix + 1;

JoannaWhittier

I am running into memory issues as described in this thread. I see where Curtis and Sebastian recommend and provide coding to create a temporary folder/workspace but am not sure how to work that into the python code I am using. To use ListRasters, I had the impression that you needed to set the env.workspace to the folder containing the raster files. My python code is below....any suggestions on how to modify this code?

I should add that sometimes this script runs and sometimes not. I think the issue is memory but just saw something about a search curser. Not sure what that is.

Thank you,
Jodi

import arcpy, os, tempfile, shutil, sys, time, traceback
from arcpy import env
from arcpy.sa import *

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Set the current workspace, this is folder that contains the .tif files for climate variables
arcpy.env.workspace = "F:\\Data\\Reg09\\tifs"

# Set Geoprocessing environments
arcpy.env.overwriteOutput = True

#Set Local variables
region =('F:\\Data\\NHD\\grids\\reg09_lambi')  # this is the folder with the NHD catchment grid (raster)
    
out = "F:\\Data\\Reg09\\update\\"
    
rasterList = arcpy.ListRasters("*", "TIF")

    
for raster in rasterList:
   name = os.path.basename(raster).strip(".tif")  # removes .tif from the file name
   arcpy.env.cellSize = "MINOF"
   arcpy.gp.ZonalStatisticsAsTable_sa(region, "VALUE", raster, out+name+".dbf","DATA", "ALL")
   print" Raster " + str(raster) + " processing completed at " + time.asctime( time.localtime(time.time()) )

SebastianRoberts

Jodi,
I found that it was not a memory leak issue, but as Curtis said:

gp/arcpy scratch name methods can be very slow when the workspace in which you are generating scratchnames gets full of many files

Creating a new workspace for OUTPUT every 30 or so feature classes greatly increased the speed of my script. You are using listrasters to create the list before you enter the loop. In your loop, you could just change the name of the output workspace after a certain number of iterations. Below is some sample code. I didn't test it, so may have some errors, but should get you started.
Sebastian

#create scratch workspace name
tempDIR = "F:\\Data\\Reg09"
tempWorkspacePrefix = "update"
tempWorkspaceSuffix = 0

count = 0;
for raster in rasterList:
   if count%30 == 0:     # this is modulus operator, so should be true the first iteration and every 30 iterations thereafter.
      tempWorkspaceSuffix = tempWorkspaceSuffix  + 1
      tempWS = tempWorkspacePrefix + str(tempWorkspaceSuffix)
      if not arcpy.Exists(tempDIR + os.sep + tempWS):
          arcpy.CreateFolder_management(tempDIR, tempWS)
      out = tempDIR + os.sep + tempWS + os.sep
   name = os.path.basename(raster).strip(".tif")  # removes .tif from the file name
   arcpy.env.cellSize = "MINOF"
   arcpy.gp.ZonalStatisticsAsTable_sa(region, "VALUE", raster, out+name+".dbf","DATA", "ALL")
   print" Raster " + str(raster) + " processing completed at " + time.asctime( time.localtime(time.time()) ).
   count = count + 1

curtvprice

My python code is below....any suggestions on how to modify this code?
I should add that sometimes this script runs and sometimes not.

What does "sometimes not" mean? Is there any kind of error message? If ArcMap is crashing (this would be a memory leak), an approach that may help make your script more stable is to use the Python garbage collector:

import gc
gc.enable()

While we're talking code, your use of "strip" may not do exactly what you want. The best way to remove file extensions is the os.path.splitext() function:

os.path.basename(os.path.splitext(path)[0])

JoannaWhittier

Thanks to both of you for the quick response. I incorporated the code provided by Sebastian and I think it is running correctly but I get the same error message as before (see below). What made me think it could be a memory issue is that the python code would run for a few files and then throw the below error message. If I delete the partially created table, and close the python shell (am running IDLE, but error also happens when just using command prompt), and the python script - I can reopen the script and run it using the very same file that threw the error. However, the number of files that were processed continued to decline and now it will not run at all. It is true that the files are read only but the script was working previously with read only files.

It may help to know that I am calculating zonal stats for NHD catchments by NHD region using the NHD regional catchment grids and a climate raster for the conterminous US - so fairly large spatial summaries.

I'll try the splittext option because you are correct Curtis, the strip option does not always work.

  File "F:\Data\Reg09\attributetifsReg09.py", line 56, in <module>
    arcpy.gp.ZonalStatisticsAsTable_sa(region, "VALUE", raster, out+name+".dbf","DATA", "ALL")

<class 'arcgisscripting.ExecuteError'>: ERROR 999999: Error executing function.
Workspace or data source is read only.
Workspace or data source is read only. [The F:\Data\Reg09\update1\ workspace is read only.]
ERROR 010067: Error in executing grid expression. Zonal statistics program failed
Failed to execute (ZonalStatisticsAsTable).

ERROR 999999: Error executing function.
Workspace or data source is read only.
Workspace or data source is read only. [The F:\Data\Reg09\update1\ workspace is read only.]
ERROR 010067: Error in executing grid expression. Zonal statistics program failed
Failed to execute (ZonalStatisticsAsTable).

curtvprice

[The F:\Data\Reg09\update1\ workspace is read only.]
ERROR 010067: Error in executing grid expression.
Failed to execute (ZonalStatisticsAsTable).

It looks to me like you are running into file-locking issues. If you are viewing the workspace in another application (say, you're running the tool in IDLE and viewing the workspace in ArcCatalog or ArcMap), and have a table open, if you try to overwrite it using your script, Windows will not give you security to do so. ArcMap for example has set a lock, and IDLE cannot overwrite the table.

If you have Windows Explorer open in the folder, you can often see the lock files that Windows writes to keep track of those locks.

These KB articles have some useful information about file locks.

Why am I getting schema locks on shapefiles?
FAQ: Why do file geodatabase .lock files remain after a process is finished

JoannaWhittier

I am still not quite understanding why there appears to be a lock that did not used to happen but have found a clumsy work around. For whatever reason, my script will run two files and then throw the read-only error. When I changed Sebastian's code from:

if count%30 == 0: # this is modulus operator, so should be true the first iteration and every 30 iterations thereafter.

to if count%2 == 0:

python will continue to run and place 2 zone stats tables in sequentially numbered new folders. Not a great solution but at least the bleeping thing is running.

Thank you again! I'll check back in case you think of something else that could lead to the lock issue I have encountered.

Follow up: my 'solution' only worked for about 10 files and the read only error returned.

HeatherWhitcomb

Hello, I'm a fairly new ArcMap 10 user. I tried using the StatisticsforOverlappingZones tool that was posted above and I got the following error. Any thoughts on what I might be doing wrong? The script ran for more than 4 hours before it returned the error message. I was hoping the fact that it was running for so long meant that it was actually executing. Does this look like a memory issue? I have 448 polygons in my shape file.

Also, how would I create an index for the join field (I'm assuming that is the FID field) in the join table?

Here was my error message:
Executing: StatisticsForOverlappingZones C:\HeatherData\CHP_inMA.shp FID C:\HeatherData\ndvi_20050512.img C:\HeatherData\Greenness_CHP_inMA.shp
Start Time: Tue May 15 12:10:10 2012
Running script StatisticsForOverlappingZones...
<class 'arcgisscripting.ExecuteError'>: Failed to execute. Parameters are not valid.
WARNING 000970: The join field FID in the join table CHP_inMA is not indexed. To improve performance, we recommend that an index be created for the join field in the join table.
ERROR 000728: Field FID does not exist within table
Failed to execute (AddJoin).

Failed to execute (StatisticsForOverlappingZones).
Failed at Tue May 15 16:34:22 2012 (Elapsed Time: 4 hours 24 minutes 12 seconds)

HeatherWhitcomb

We recently published a toolbox (first prize at the ESRI UC 2011 AppFair!) that includes tools that do overlap processing for statistics and area tabulations, among other things.

USGS Open-File Report 2010�??1268
National Water-Quality Assessment (NAWQA) Area-Characterization Toolbox
By Curtis V. Price, Naomi Nakagaki, and Kerie J. Hitt

http://pubs.usgs.gov/of/2010/1268/

The published version has only limited support for 10.x, but I have an updated version that will be posted there soon. If you are interested, let me know and I can send it your way.

Hi Curtis,

I am interested in checking out the tool for use with 10. Is it available on the web?

Heather

LuciaBarbato

Hi Curtis,

We've been following this thread. We too, are having issues with zonal statistics as table with overlapping polygons. We have a very large dataset of over 80,000 land use features with about 2,500 overlaps. While the process runs, it does not output results for the overlapping features. We are using 10.0 SP3. Is there a link to an updated tool that we could download?

Thanks,

Lucia

curtvprice

Is there a link to an updated tool that we could download?

I've been posting updated versions here for review, please use at your own risk and send comments:

ftp://ftpext.usgs.gov/pub/nonvisible/cr/nact

LuciaBarbato

Thanks for the quick reply! We will work on it today and report back.

LuciaBarbato

We just set up and started running the feature statistics to table tool. We checked the box to process overlapping features. After 14 minutes, 31 features of 79,615 were processed. We figure that's about 600 hours of processing. We read in previous posts about setting up scratch workspaces by modifying the code, but have reached the end of our expertise. We just wanted to let you know how things worked out for us. We will continue our work and deal with the overlapping polygons at a later time. Thanks for letting us try your code.

curtvprice

> I am interested in checking out the tool for use with 10. Is it available on the web?

A version of the NAWQA toolbox that runs under 10.x is available in ArcGIS Online here:

http://www.arcgis.com/home/item.html?id=cbb59504f59f4e18b23817fb0ef40e56

Doe mee aan het gesprek

Aangemelde leden kunnen berichten plaatsen, updates volgen en meer. Nieuw hier? Registreer een gratis account.

Snelle Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.