topic Re: Zonal statistics - percentiles in ArcGIS Spatial Analyst Questions

Zonal statistics - percentiles

toddsams — Mon, 30 Jun 2014 20:35:35 GMT

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?

Re: Zonal statistics - percentiles

curtvprice — Wed, 02 Jul 2014 19:28:36 GMT

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.

The Slice tool is a handy way to do this. If you specify an EQUAL_AREA slice with 20 zones specified, the first class will be in the 5th percentile, the last will be in the 95th.

If your data values are not continuously distributed (that is, many ties, which can make the percentile value less useful) or highly skewed (which can get you into numerical issues) you may not get the results you expect, but if the data are well-distributed this approach should work very well.

Re: Zonal statistics - percentiles

toddsams — Wed, 02 Jul 2014 20:28:00 GMT

Thanks for the Slice suggestion.

I ran the tool using the parameters you suggested. However, the raster attribute table that was returned simply shows 20 rows with "Values" of 1 through 20 and a "Count" of the number of pixels in each of these equal area classes.

I don't see how I can derive the percentile values from this information?

Re: Zonal statistics - percentiles

curtvprice — Thu, 10 Jul 2014 21:22:48 GMT

You then analyze your sliced raster with Zonal Statistics As Table to get the min and max for each class.

Re: Zonal statistics - percentiles

toddsams — Thu, 10 Jul 2014 21:56:11 GMT

It may be that my slice raster was not generated properly (attribute table shown below). It looks like it returned the number of cells in each class, which are all about the same since I used equal area classes.

I don't see how to analyze this sliced raster with zonal statistics to generate a table of percentiles?

Re: Zonal statistics - percentiles

curtvprice — Wed, 23 Jul 2014 19:09:59 GMT

I'm sorry Todd - still working on how to track threads in this new system, hence no reply til now.

If you ask for 20 classes from the Slice tool with the equal area option, the first class is the 5th percentile, ie VALUE times 5 (100 / 20).

If you want the percentile values (breaks), you can run Zonal Statistics As Table using the sliced raster for zones and the input to the slice for values.

Re: Zonal statistics - percentiles

XanderBakker — Fri, 25 Jul 2014 19:45:38 GMT

Hi Todd,

Maybe a little Python can help: Look at the code below:

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


import arcpy

# path to raster (has to be integer!)
ras = r"C:\Path\To\Folder\Or\Geodatabase\With\Raster"

# create dictionary with value vs count
flds = ("VALUE", "COUNT")
dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(ras, flds)}

# calculate number of pixels in raster
cnt_sum = sum(dct.values())

dct_per = {}
cnt_i = 0

# loop through dictionary and create new dictionary with val vs percentile
for val in sorted(dct.keys()):
    cnt_i += dct[val]
    dct_per[val] = cnt_i / cnt_sum

# use dct_per to determine percentiles
perc = 0.05 # (5%)
pixval =  GetPixelValueForPercentile(dct_per, perc)
print "Found for percentile {0}, pixel value {1} with real percentile {2}".format(perc, pixval, dct_per[pixval])

perc = 0.95 # (95%)
pixval =  GetPixelValueForPercentile(dct_per, perc)
print "Found for percentile {0}, pixel value {1} with real percentile {2}".format(perc, pixval, dct_per[pixval])

It will return the following information:

>> Found for percentile 0.05, pixel value 401 with real percentile 0.049861616298

>> Found for percentile 0.95, pixel value 2350 with real percentile 0.949769065993

The code requires that your raster is of type integer (to make sure it has an attribute table with pixelvalue and count.

This is what happens in the code:

on line 16 a path to your integer raster is defined (e.g. ras = r'D:\Xander\tmp\ras.gdb\DEMsmallInt')
on line 20 the attribute table is used by a da search cursor to create a dictionary (this helps to speed up the proces)
on line 23 the sum of the count is calculated (the number of pixel with a value in the raster)
on line 29 a dictionary is created that holds the pixel value and corresponding percentile
on line 35 and 39 the function (defined at the top) is used to search the pixel value for a percentile value.

Hope this is of any use.

Kind regards, Xander

Re: Zonal statistics - percentiles

toddsams — Tue, 29 Jul 2014 17:37:03 GMT

Thanks for the Python code for percentiles from a raster Xander.

In my situation, I have multiple polygons that occur in various locations across a raster. I would like to generate the percentile statistics from the various portions of the raster that overlap with each individual polygon. Then, join those statistics with a unique row ID in the polygon layer.

I am essentially in need of "Zonal Statistics as Table" with the ability to include percentile statistics along with Minimum, Mean, and Maximum.

How can this be incorporated via Python?

Re: Zonal statistics - percentiles

XanderBakker — Tue, 29 Jul 2014 22:11:28 GMT

Hi Todd,

One option is that you write a loop through the polygon (Update cursor) to do this for each polygon:

use the arcpy.sa.ExtractByPolygon tool to extract the relevant portion of raster and save that to the IN_MEMORY workspace
use that raster and the previous posted code to determine the percentiles
define a list of percentile values like lst_prec = [2, 5,10, 90, 95, 98] to loop through the values and determine the raster value that corresponds to that percentile. Store the result in previously created fields.

In case there is no overlap between the polygons, you could do the following:

Rasterize your polygon layer
Use the arcpy.sa.Combine tool to create a raster that holds the combination of both values. If the number of unique combinations is less than 65536, this should be no problem.
Adapt the code to read both the polygon layer in an update cursor and the combine raster attribute table in a search cursor to determine the percentiles and store the result

Not sure which option will be the fastest.

Now to throw in some code (which by the way I did not test), it should look something like this:

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

def FieldExist(featureclass, fieldname):
    """Check if field exists"""
    import arcpy
    fieldList = arcpy.ListFields(featureclass, fieldname)
    return len(fieldList) == 1


import arcpy
# settings
ras = r"C:\Path\To\Folder\Or\Geodatabase\With\Raster"
fc = r"C:\Path\To\Folder\Or\Geodatabase\With\PolygonFeatureclass"
lst_perc = [2, 5, 10, 90, 95, 98] # list of percentiles to be calculated
fld_prefix = "Perc_"

arcpy.env.workspace = "IN_MEMORY"

# add fields if they do not exist in fc
# be aware, field that exist will be overwritten!
flds = []
for perc in lst_perc:
    fld_perc = "{0}{1}".format(fld_prefix, perc)
    flds.append(fld_perc)
    if not FieldExist(fc, fld_perc):
        arcpy.AddField_management(fc, fld_perc, "LONG")

# Enable Spatial analyst
arcpy.CheckOutExtension("Spatial")

# loop through polygons
flds.append("SHAPE@")
with arcpy.da.UpdateCursor(fc, flds) as curs:
    for row in curs:
        polygon = row[flds.index("SHAPE@")]
        lst_parts = []
        for part in polygon:
            lst_crds = []
            for crdpair in part:
                x, y = crdpair[0], crdpair[1]
                lst_crds.append(arcpy.Point(x, y))
            lst_parts.append(lst_crds)

        # Execute ExtractByPolygon
        ras_pol = arcpy.sa.ExtractByPolygon(ras, polygon, "INSIDE")
        # not sure if ras_pol should be saved...

        # create dictionary with value vs count
        flds_ras = ("VALUE", "COUNT")
        dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(ras_pol, flds_ras)}

        # calculate number of pixels in raster
        cnt_sum = sum(dct.values())

        # loop through dictionary and create new dictionary with val vs percentile
        dct_per = {}
        cnt_i = 0
        for val in sorted(dct.keys()):
            cnt_i += dct[val]
            dct_per[val] = cnt_i / cnt_sum

        # loop through list of percentiles
        for perc in lst_perc:
            # use dct_per to determine percentiles
            prec_dec = float(perc) / 100
            pixval =  GetPixelValueForPercentile(dct_per, perc)

            # write pixel value to percentile field
            fld_perc = "{0}{1}".format(fld_prefix, perc)
            row[flds.index(fld_perc)] = pixval

        # update row
        curs.updateRow(row)

# return SA license
arcpy.CheckInExtension("Spatial")

If you run into problems, let me know. Maybe you can throw in a piece of your data to test the code.

Kind regards,

Xander

Re: Zonal statistics - percentiles

XanderBakker — Wed, 30 Jul 2014 13:25:13 GMT

OK, so far for posting code without testing it ...

I created a small test dataset, consisting of a DEM (integer raster) and some random polygons. During the test a number of errors in the code were revealed. These are corrected in the code below:

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

def FieldExist(featureclass, fieldname):
    """Check if field exists"""
    import arcpy
    fieldList = arcpy.ListFields(featureclass, fieldname)
    return len(fieldList) == 1


import arcpy
# settings
ras = r"D:\Xander\tmp\ras.gdb\DEMsmpro"
fc = r"D:\Xander\tmp\ras.gdb\pols"
# ras = r"C:\Path\To\Folder\Or\Geodatabase\With\Raster"
# fc = r"C:\Path\To\Folder\Or\Geodatabase\With\PolygonFeatureclass"
lst_perc = [2, 5, 10, 90, 95, 98] # list of percentiles to be calculated
fld_prefix = "Perc_"

arcpy.env.workspace = "IN_MEMORY"
arcpy.env.overwriteOutput = True

# add fields if they do not exist in fc
# be aware, field that exist will be overwritten!
print "filling field list and adding fields"
flds = []
for perc in lst_perc:
    fld_perc = "{0}{1}".format(fld_prefix, perc)
    flds.append(fld_perc)
    if not FieldExist(fc, fld_perc):
        arcpy.AddField_management(fc, fld_perc, "LONG")
print "flds={0}".format(flds)

# Enable Spatial analyst
arcpy.CheckOutExtension("Spatial")

# loop through polygons
print "loop through polygons"
flds.append("SHAPE@")
i = 0
with arcpy.da.UpdateCursor(fc, flds) as curs:
    for row in curs:
        i += 1
        print "Processing polygon: {0}".format(i)
        polygon = row[flds.index("SHAPE@")]
        lst_parts = []
        if polygon.partCount == 1:
            for part in polygon:
                for pnt in part:
                    x, y = pnt.X, pnt.Y
                    lst_parts.append(arcpy.Point(x, y))
        else:
            for part in polygon:
                lst_crds = []
                for pnt in part:
                    x, y = pnt.X, pnt.Y
                    lst_crds.append(arcpy.Point(x, y))
                lst_parts.append(lst_crds)
            print " - lst_parts={0}".format(lst_parts)

        # Execute ExtractByPolygon (you can't send the polygon object)
        print " - ExtractByPolygon..."
        ras_pol = arcpy.sa.ExtractByPolygon(ras, lst_parts, "INSIDE")
        outname = "ras{0}".format(i)
        ras_pol.save(outname)
        print " - saved raster as {0}".format(outname)

        # create dictionary with value vs count
        print " - fill dict with Value x Count"
        flds_ras = ("VALUE", "COUNT")
        dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(outname, flds_ras)}

        # calculate number of pixels in raster
        print " - determine sum"
        cnt_sum = sum(dct.values())
        print " - sum={0}".format(cnt_sum)

        # loop through dictionary and create new dictionary with val vs percentile
        print " - create percentile dict"
        dct_per = {}
        cnt_i = 0
        for val in sorted(dct.keys()):
            cnt_i += dct[val]
            dct_per[val] = cnt_i / cnt_sum

        # loop through list of percentiles
        print " - iterate perceniles"
        for perc in lst_perc:
            # use dct_per to determine percentiles
            perc_dec = float(perc) / 100
            print "   - Perc_dec for is {0}".format(perc_dec)
            pixval =  GetPixelValueForPercentile(dct_per, perc_dec)
            print "   - Perc for {0}% is {1}".format(perc, pixval)

            # write pixel value to percentile field
            print "   - Store value"
            fld_perc = "{0}{1}".format(fld_prefix, perc)
            row[flds.index(fld_perc)] = pixval

        # update row
        print " - update row"
        curs.updateRow(row)

# return SA license
arcpy.CheckInExtension("Spatial")

What you get are a number of fields in your polygon featureclass:

When you add some labels to display the values, you'll see that they actually correspond to the areas:

So Todd, I hope this will work for you.

Kind regards,

Xander

Re: Zonal statistics - percentiles

toddsams — Thu, 04 Sep 2014 17:14:06 GMT

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

Re: Zonal statistics - percentiles

toddsams — Wed, 14 Jan 2015 22:40:45 GMT

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?

Re: Zonal statistics - percentiles

XanderBakker — Thu, 15 Jan 2015 00:35:22 GMT

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.

Re: Zonal statistics - percentiles

toddsams — Thu, 15 Jan 2015 18:23:45 GMT

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

Re: Zonal statistics - percentiles

XanderBakker — Sat, 11 Dec 2021 20:05:24 GMT

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

Re: Zonal statistics - percentiles

toddsams — Thu, 15 Jan 2015 19:05:44 GMT

That works! Many thanks for your help with this.

Re: Zonal statistics - percentiles

XanderBakker — Thu, 15 Jan 2015 19:16:01 GMT

You're welcome. Glad it worked!

Re: Zonal statistics - percentiles

toddsams — Thu, 15 Jan 2015 19:58:53 GMT

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.

Re: Zonal statistics - percentiles

XanderBakker — Thu, 15 Jan 2015 21:02:59 GMT

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.

Re: Zonal statistics - percentiles

toddsams — Thu, 15 Jan 2015 22:02:20 GMT

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.