Zonal statistics - percentiles

[toddsams]

Hello,

I am looking for an analogous method to Zonal Statistics that is able to generate user-defined percentiles (such as the 5th and 95th percentile) of the input raster data.

I looked at the GME tools and those stats metrics are hardwired in the same fashion as ArcGIS Zonal Statistics tools are.

Any other ideas that I can explore?

[toddsams]

Thanks Xander. I will work with this and attempt to modify for iterating through a list of rasters.

[toddsams]

Hi Xander,

I was not previously able to successfully execute this code, but am now making a second attempt...

I have replaced lines 21 and 22 with my raster and polygons and am receiving "RuntimeError: cannot open 'ras1'" while creating the dictionary at line 79.

Do you have any thoughts on why I am encountering this error?

[XanderBakker]

Can you post a part of your data? It will be difficult to trace the error without it. Is the DEM an integer or is it floating? If it is floating then it will not have an attribute table and the code will fail.

The other things is, today I tried an da updatecursor on an integer raster and it didn't work, while the old update cursor did work. Not sure if that should be changed.

[toddsams]

The raster was floating. I switched to integer and was able to make it further, but was hung up with:

UnboundLocalError: local variable 'pix_val' referenced before assignment

File "D:\Scripts\Python\Raster_Percentiles.py", line 100, in <module>

  pixval =  GetPixelValueForPercentile(dct_per, perc_dec)

File "D:\Scripts\Python\Raster_Percentiles.py", line 10, in GetPixelValueForPercentile

  return pix_val

[XanderBakker]

This would happen when the percentile you're looking for is higher than the first percentile stored in the dictionary. I guess this could be solved by replacing the def GetPixelValueForPercentile(dctper, percentile) for a slightly changed versión (not tested):

def GetPixelValueForPercentile(dctper, percentile):
    """will return last pixel value
      where percentile LE searched percentile."""
    pix_val = sorted(dctper.keys())[0] # initially assign lowest pixel value
    for k in sorted(dctper.keys()):
        perc = dctper
        if perc <= percentile:
            pix_val = k
        else:
            break
    return pix_val

[toddsams]

That works! Many thanks for your help with this.

[XanderBakker]

You're welcome. Glad it worked!

[toddsams]

I ultimately need to iterate through a list of rasters (or a composite raster) to produce multiple sets of percentile results.

My first attempt was to create a list of the rasters and use this in a for loop around lines 21 through 110. However, there were two problems with this:

1) each iteration overwrites (rather than appends to) the previous set of percentile results

2) to keep track of the results, a new field is needed in the output that contains the name of the raster for which the percentiles were derived from

Are you able to provide a solution for this?

I was thinking that output to a non-spatial table with a PolygonID may be an alternate method.

[XanderBakker]

I suppose one could create a related table using the raster value as key so that there can be n raster related to each polygon. It really depends on how you want to use the results.

[toddsams]

Whether the results are in an attribute table associated with a polygon or a non-spatial table doesn't really matter. I would export either table to a spreadsheet to summarize.