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

Encontre mais publicações marcadas com

Accepted answers

All comments

curtvprice

Are there plans to add this significant enhancement of the tool to ArcMap 10.8? Is it even a possibility?

SarmisthaChatterjee

Thank you for your question and this discussion.
We are happy to announce that, in ArcGIS Pro 2.5 you can now use the Zonal Statistics as Table tool to calculate statistics for overlapping polygon zones. The zonal analysis will be performed for each individual polygon.

ThomasStanley

I gave Zonal Statistics as Table 2 from v 1.3 a try. It was very easy to use, but it still skips a lot of polygons. (Only 3279 out of 3550 buffers made it into the final table.) It's too bad, because I ran the Model Builder method overnight and it didn't complete.

curtvprice

Esri has added a new tool to the Spatial Analyst Supplemental Tools that works with overlapping polygons. Haven't looked at it yet but I'm sure it's worth a try.

ArcGIS Blog: New Spatial Analyst Supplemental tools, v1.2

Spatial Analyst Supplemental Tools

curtvprice

Zonal statistics could not be computed for feature 26 (zone = 26)
ERROR 000601: Cannot delete C:\Research\Projects\xxxwk1\xyz.gdb\workfeat. May be locked by another application.
Failed to execute (CopyFeatures).

Manisha, I don't think the above error has do with my tool. There have been reports of antivirus software locking ArcGIS files while they are being scanned, and making things miserable for everyone.

I honestly don't have any good ideas for you. Can you try on a different machine that may have different security set up?

manisharai

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

Hi Curtis,

I am a PhD student and was using the NAWQA tool for ArcGIS 9.3 for my research. However, I upgraded to version 10.1 and couldn't run the NAWQA tool anymore. So, I downloaded the file that you listed in your post for version 10.x. I am trying to calculate mean elevations around some sites in CA. I use the Feature statistics to table tool from the NACT toolbox and set the buffer distance to say 500m. The tool runs fine for some sites and then I start getting these error messages :

Zonal statistics could not be computed for feature 26 (zone = 26)
ERROR 000601: Cannot delete C:\Research\Projects\xxxwk1\xyz.gdb\workfeat. May be locked by another application.
Failed to execute (CopyFeatures).

I do not have any other application running (arc catalogue, IDLE) that could possibly be causing a lock. I have no experience with coding in ArcGIS so I am not able to figure out the problem on my own by looking at the code.

I hope you could help me with this.

Thanks.

Manisha

SteveLynch

I did, but I get errors

Traceback (most recent call last):
File "Z:\Extraction of Bands\Extraction\Python for extraction\Trial\ZonalStatsOverlappingPolys (2)\makeZones.py", line 170, in <module>
arcpy.CopyFeatures_management(inFeat, inFeatures)
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\management.py", line 2226, in CopyFeatures
raise e
ExecuteError: ERROR 000210: Cannot create output Z:\Extraction of Bands\Extraction\Python for extraction\Trial\Scratch.mdb\50m_buffer_P_selection
ERROR 000361: The name starts with an invalid characterFailed to execute (CopyFeatures).

Failed to execute (ZonesWOverlap).
Failed at Wed Oct 31 10:37:07 2012 (Elapsed Time: 8.00 seconds)

Output in a geodatabase must start with an alphabetic character, yours starts with the number 5

BrendanCollins

Hey can people vote for this idea to enhance zonal tools:

http://ideas.arcgis.com/ideaView?id=0873000000086QNAAY

MarkBoucher

I did finish the tool I was working on. The tool is faster than the individual polygon solution. Originally I produced 3 different solutions for breaking up the polygons and ended up using two of those methods. I have uploaded the tool and documented it.

The tool uses Spatial Join and also breaks up the feature class into smaller pieces.

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

Every one of my polygons overlap at least one other polygon. Will this tool work for me? Basically, every polygon would have to be process separately, correct? The grouping technique would not work.

JamieKass

Even if you used my current (10x-compatible) version (ftp://ftpext.usgs.gov/pub/nonvisible/cr/nact) it would not be any faster than Jamie's code,as i also do one at a time. Have you tried Sue's tool? I haven't tried it yet, but it seems like a really good solution.

My original code does not process one feature at a time, rather it attributes all features with an "explode" field that breaks them up into non-overlapping groups, which Zonal Statistics can then be iterated over. My method involves as few arcpy functions as possible to save time, as processing large features often tie up geoprocessing tools that are managing spatial characteristics as well as table data. Essential functions like Identity are run, then cursors pull table data into Python dictionaries, where most of the processing takes place. These functions were built to handle large datasets, so they should perform well even with many features.
I have a newer version of the code I posted last time that's a bit faster (uses da cursor functions and iterates over the dicts less) and cleaner, which I'll post below. It comes in a couple of functions, which can be called separately if need be. If anyone needs help implementing it, please post to this thread.

def makeOverlapDict(fc, returnID = None):
    """This function creates a dictionary relating each feature to a list of its
    overlapping features. The optional returnID parameter can be used to create
    an overlapDict with keys equal to a different unique ID than OID, which is the
    default."""
    fcName = arcpy.Describe(fc).name
    if not returnID:
        returnID = 'FID_' + fcName
        returnID2 = 'FID_' + fcName + '_1'
    else:
        returnID2 = returnID + '_1'

    identity = arcpy.Identity_analysis(fc, fc, fcName + 'ID')

    overlapDict = {}
    with arcpy.da.SearchCursor(identity, (returnID, returnID2)) as cursor:
        for row in cursor:
            if row[0] <> row[1]:
                if not overlapDict.get(row[0]):
                    overlapDict[row[0]] = [row[1]]
                else:
                    overlapDict[row[0]].append(row[1])

    arcpy.Delete_management(identity)

    return overlapDict

def makeExplodeDict(overlapDict):
    """This function creates a new dictionary relating each feature to
    its explode value based on the overlapDict."""

    # set initial explode value to 1
    expl = 1
    # create dictionary to hold explode values of oids
    explDict = {}
    overlapDict2 = overlapDict.copy()
    # while there are still keys in overlapDict2 not also in explDict
    while len([k for k in overlapDict2.iterkeys() if not explDict.get(k)]) > 0:
        # create list to hold overlapping oids
        ovSet = set()
        # loop over overlapDict2
        for oid, ovlps in overlapDict2.iteritems():
            if not explDict.get(oid) and oid not in ovSet:
                explDict[oid] = expl
                for o in ovlps:
                    ovSet.add(o)
        for oid in explDict.iterkeys():
            if overlapDict2.get(oid):
                del overlapDict2[oid]
        print expl
        expl += 1
    return explDict

def defineOverlaps(fc):
    """This function creates three fields: 'overlaps', which lists the OBJECTIDs
    of all other polygons overlapping each polygon, and 'ovlpCount', which counts
    the number of overlaps per feature, and 'expl', which separates all features
    into unique non-overlapping groups by value."""

    # Build dicts used to relate each feature to all identical, overlapping shapes
    overlapDict = makeOverlapDict(fc)

    # Loop through overlapDict and give each oid an expl value that separates all oids into
    # non-overlapping groups
    explDict = makeExplodeDict(overlapDict)

    # Add remaining non-overlapping oids into explDict and give them an expl value
    # equal to the max expl; the max should have the least members, and this preserves
    # a good distribution of values so that each processing step is somewhat balanced
    maxExpl = max(explDict.itervalues())
    query = str(tuple(explDict.keys()))
    fillVals = []
    with arcpy.da.SearchCursor(fc, ("OID@"), '"OBJECTID" NOT IN ' + query) as cursor:
        for row in cursor:
            fillVals.append(row[0])
    for val in fillVals:
        explDict[val] = maxExpl

    # Add overlaps, ovlpCount, and expl fields
    arcpy.AddField_management(fc,'overlaps',"TEXT",'','',1000) # the 1000 is the field character limit... increase if need be
    arcpy.AddField_management(fc,'ovlpCount',"SHORT")
    arcpy.AddField_management(fc,'expl',"SHORT")

    # For each row with an overlapping poly (i.e. with an entry in idDict)
    # give overlaps value of all overlapping ids, ovlpCount value of how
    # many overlaps occur, and expl value that separates all features into
    # unique non-overlapping groups (explode value)

    with arcpy.da.UpdateCursor(fc,("OID@","overlaps","ovlpCount","expl")) as cursor:
        for row in cursor:
            # If oid in overlapDict, write overlapping oids and their count
            # to the appropriate fields
            if overlapDict.get(row[0]):
                row[1] = str(overlapDict[row[0]]).strip('[]')
                row[2] = len(overlapDict[row[0]])
            # Write expl value to expl field
            row[3] = explDict[row[0]]
            cursor.updateRow(row)

If you're using the Python window, once these 3 function definitions have been entered, simply call the operation like this, where 'fc' is your feature class:

defineOverlaps(fc)

Please note that the next step would be to select by attributes based on the 'expl' value and run Zonal Statistics for each value, which will be at most in the tens as long as the features are not too piled up. This can be performed easily with a simple for loop.

Xin_RanDing

I believe the error (from CopyFeatures, not Sue's script) is telling you to avoid naming datasets starting with a number. I also suggest avoiding spaces in pathnames, it works most of the time but is just asking for trouble -- something will break when you most need it not to.

000361 : The name starts with an invalid character

Hi

I renamed the files and got past that stage but now I'm getting an error of during zonal statistics, apparently it can't read the raster file

Executing: ZonesWOverlap S_buffer OBJECTID Band_1.tif "Z:\Extraction of Bands\Trial Extraction\Trail1\StatisticsForOverlappingZones" F:\Data\xiding\Documents\ArcGIS\Default.gdb\S_buffer_ZonesWOverlap 7
Start Time: Thu Nov 01 13:53:15 2012
Running script ZonesWOverlap...
Selecting the [fCount] = 0...
Feature layer number 1 of 7
Iteration 1, Spatial Join
Iteration 1, Join
Select S_buffer.OBJECTID = j1.OBJECTID_1
Removing Join...
Iteration 1, Zonal Statistics
<class 'arcgisscripting.ExecuteError'>: ERROR 010160: Unable to open raster zLayer. Zonal statistics program failed
ERROR 010067: Error in executing grid expression.
Failed to execute (ZonalStatisticsAsTable).

Failed to execute (ZonesWOverlap).
Failed at Thu Nov 01 13:54:01 2012 (Elapsed Time: 46.00 seconds)

curtvprice

ERROR 000361: The name starts with an invalid character
Failed to execute (CopyFeatures).

I believe the error (from CopyFeatures, not Sue's script) is telling you to avoid naming datasets starting with a number. I also suggest avoiding spaces in pathnames, it works most of the time but is just asking for trouble -- something will break when you most need it not to.

000361 : The name starts with an invalid character

Xin_RanDing

Even if you used my current (10x-compatible) version (ftp://ftpext.usgs.gov/pub/nonvisible/cr/nact) it would not be any faster than Jamie's code,as i also do one at a time. Have you tried Sue's tool? I haven't tried it yet, but it seems like a really good solution.

I did, but I get errors

Traceback (most recent call last):
File "Z:\Extraction of Bands\Extraction\Python for extraction\Trial\ZonalStatsOverlappingPolys (2)\makeZones.py", line 170, in <module>
arcpy.CopyFeatures_management(inFeat, inFeatures)
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\management.py", line 2226, in CopyFeatures
raise e
ExecuteError: ERROR 000210: Cannot create output Z:\Extraction of Bands\Extraction\Python for extraction\Trial\Scratch.mdb\50m_buffer_P_selection
ERROR 000361: The name starts with an invalid character
Failed to execute (CopyFeatures).

Failed to execute (ZonesWOverlap).
Failed at Wed Oct 31 10:37:07 2012 (Elapsed Time: 8.00 seconds)

AudraPayne1_2-1740422807555.png

curtvprice

Hi

I'm trying to extraction band information from raster images using overlapping buffers. I tried using Jamie's code and it takes a long time. I can't ran Curtis' code. Does anymore have a code that works faster?

Thanks

Even if you used my current (10x-compatible) version (ftp://ftpext.usgs.gov/pub/nonvisible/cr/nact) it would not be any faster than Jamie's code,as i also do one at a time. Have you tried Sue's tool? I haven't tried it yet, but it seems like a really good solution.

Xin_RanDing

Hi

I'm trying to extraction band information from raster images using overlapping buffers. I tried using Jamie's code and it takes a long time. I can't ran Curtis' code. Does anymore have a code that works faster?

Thanks

AbbyFlory

I did finish the tool I was working on. The tool is faster than the individual polygon solution. Originally I produced 3 different solutions for breaking up the polygons and ended up using two of those methods. I have uploaded the tool and documented it.

The tool uses Spatial Join and also breaks up the feature class into smaller pieces.

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

Hi Sue,

Thanks for sharing your code! I was wondering though, would you mind explaining the logic behind the following sentence taken from your official tool post (your link above) under Description > Solution: "If a feature class is spatially joined to itself, the polygon id in the target and join will match in the top-most layer of polygons. This set of polygons will not overlap any other polygons." I am specifically confused by the phrase "top-most layer of polygons," and how this process produces non-overlapping features.

I tried manual testing to better understand it, and when I spatially joined a feature class to itself, the output was identical to the original feature class (with the exception of added fields from the join). Further, joining this back to the original feature class on polygon id did nothing but add more fields - as all the polygons matched/joined because they were all the same as the originals.

I have reviewed your code and am still stumped by this process. If you have the time and wouldn't mind going into the logic details, I'd love to understand what's happening here. I have been trying to implement a similar routine (but not using spatial join) that results in grouping points (via a flag field) into non-overlapping subsets based on either their buffer or distance to other points. I'm needing it for the Tabulate Area tool, which has the same overlapping areas problem as the Zonal Stats as Table tool. Any help would be appreciated! Thanks for your time!

Abby

AbbyFlory

Hi Sue,

Thanks for sharing your code! I was wondering though, would you mind explaining the logic behind the following sentence taken from your official tool post under Description > Solution: "If a feature class is spatially joined to itself, the polygon id in the target and join will match in the top-most layer of polygons. This set of polygons will not overlap any other polygons." I am specifically confused by the phrase "top-most layer of polygons," and how this process produces non-overlapping features.

I tried manual testing to better understand it, and when I spatially joined a feature class to itself, the output was identical to the original feature class (with the exception of added fields from the join). Further, joining this back to the original feature class on polygon id did nothing but add more fields - as all the polygons matched/joined because they were all the same as the originals.

I have reviewed your code and am still stumped by this process. If you have the time and wouldn't mind going into the logic details, I'd love to understand what's happening here. I have been trying to implement a similar routine (but not using spatial join) that results in grouping points (via a flag field) into non-overlapping subsets based on either their buffer or distance to other points. I'm needing it for the Tabulate Area tool, which has the same overlapping areas problem as the Zonal Stats as Table tool. Any help would be appreciated! Thanks for your time!

Abby

PeterWilson

As promised

please find attached my Toolbox (ArcGIS 10.0) with attached model that iterates through each feature within a feature class and passes it onto zonal statistics to calculate the average slope for each polygon (Watershed Catchment) and outputs it into a single table.

Please note that I've used the following model to process over 1200 polygons with no problems. It runs pretty quickly.

I'll post the updated model which identifies the overlapping features and processes them separately to minimise the amount of iterates shortly. problems with my internet connection at the office.

Regards

IterateZonalStatistics.zip

PeterWilson

I wouldn't say it's a thing of the past. Setting up the model is not trivial for people new to ModelBuilder, especially how you aggregate the results. Please do post your model, IMHO this is a good "real-world" example of how iteration can be used to great advantage in Arc 10.

The tool Sue just posted is better than a brute-force iteration though - it separates your overlapping polygons into non-overlapping groups - if you have thousands of overlapping polygons, this could mean a few hundred iterations instead of thousands.

Hi Curtis

I agree that using a iterate to deal with the overlapping is a bit of a brute force, but as a suggestion I could alter the tool to only deal with feature that overlap and pass them as single "in_memory" process" to zonal statistics thereby eliminating the amount of iterations required. I'll post my original model tonight when I get home from work as well as an adjusted model to deal with the overlaps only.

Regards

curtvprice

The problem with overlapping polygons for Zonal Statistics is of the past. If you iterate over your features within a feature class and pass each one to zonal statists it works perfectly and writes out the results into a single table. If anyone would like I would gladly post my model.

I wouldn't say it's a thing of the past. Setting up the model is not trivial for people new to ModelBuilder, especially how you aggregate the results. Please do post your model, IMHO this is a good "real-world" example of how iteration can be used to great advantage in Arc 10.

The tool Sue just posted is better than a brute-force iteration though - it separates your overlapping polygons into non-overlapping groups - if you have thousands of overlapping polygons, this could mean a few hundred iterations instead of thousands.

PeterWilson

Thats all I get from this tool.

<type 'exceptions.RuntimeError'>: Row: Field FID does not exist
Failed to execute (StatisticsForOverlappingZones)

Hi everyone

With the release of ArcGIS 10 a list of iterate functions have been added to ModelBuilder. The problem with overlapping polygons for Zonal Statistics is of the past. If you iterate over your features within a feature class and pass each one to zonal statists it works perfectly and writes out the results into a single table. If anyone would like I would gladly post my model.

Regards

bodoMalowerschnig

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

Thats all I get from this tool.

<type 'exceptions.RuntimeError'>: Row: Field FID does not exist
Failed to execute (StatisticsForOverlappingZones)

SueClark

I did finish the tool I was working on. The tool is faster than the individual polygon solution. Originally I produced 3 different solutions for breaking up the polygons and ended up using two of those methods. I have uploaded the tool and documented it.

The tool uses Spatial Join and also breaks up the feature class into smaller pieces.

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

JamieKass

If anyone is interested, I developed some functions in arcpy that break up the overlapping features into distinct non-overlapping groups. After the functions are run, one simply has to iterate over all the unique values in the "expl" field and run Zonal Stats on each. This method reduces the time spent iterating considerably, as the number of non-overlapping groups is a mere fraction of the total number of features. Please see the answer with the green check mark.

http://gis.stackexchange.com/questions/32217/explode-overlapping-polygons-to-new-non-overlapping-features

SueClark

Thanks Curt! I'm still relatively new to arcpy. The code I published is a special case. I've been trying to make it run as a fully implemented tool in my limited spare time. Right now I'm having some problems getting the tool to run with a FGDB raster.

Sue

curtvprice

Neat, Sue! I will definitely look into this - thanks for sharing your code.

I noticed something you may be interested in:

Here's your (nicely documented) code:

    # Note the select string will be different depending on the workspace
    desc = arcpy.Describe(env.workspace)
    if desc.extension == 'mdb':
        selStr = "[" + fCount + "] = " 
    else:
        selStr = "\"" + fCount + "\" = "

Esri has provided us a more

'>robust arcpy method to solve this problem:

    # Note the select string will be different depending on the workspace
    dfCount = arcpy.AddFieldDelimiters(env.workspace,fCount)
    selStr = dfCount + " = "

SueClark

I was working on separating a feature class with overlapping polygons into a minimal set of feature classes with non-overlapping polygons when I was asked to work on a solution for the Zonal Statistics problem described here. I came up with a different solution, probably because I was working on this from a different angle. The routine I wrote is recursive and uses spatial join.

I'm only attaching the python code. Feel free to use it or modify it for your needs.

Sue

The code wasn't as efficient as I would like, so I rewrote it. It's no longer recursive. Although you might have to think about the methodology a bit, I think the code is a little clearer. The code requires one spatial join and uses the joined id field to seperate the features into non-overlapping groups.

makeZones.py.zip

SueClark

I uploading the full tool to arcgis online, http://www.arcgis.com/home/item.html?id=b859b33c616a47d2b99b5e133942db02

Sue

SueClark

I was working on separating a feature class with overlapping polygons into a minimal set of feature classes with non-overlapping polygons when I was asked to work on a solution for the Zonal Statistics problem described here. I came up with a different solution, probably because I was working on this from a different angle. The routine I wrote is recursive and uses spatial join.

I'm only attaching the python code. Feel free to use it or modify it for your needs.

Sue

makeZones.py.zip

curtvprice

How do you get this installed into ArcMap10??

In Arc 10, you can just access the toolbox directly in the catalog window. No need to "add toolbox".

StevenDempsey

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

How do you get this installed into ArcMap10??

Participe da Conversa

Membros conectados podem postar, seguir atualizações e mais. Novo aqui? Registre uma conta gratuita.

Links Rápidos

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