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?

XanderBakker

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