Zonal statistics - percentiles

28038
65
Jump to solution
06-30-2014 01:35 PM
toddsams
New Contributor III
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?
65 Replies
AxelThomas
New Contributor III

Hi Xander,

perhaps the original raster I have been working with would be helpful (please see the attachment).

Regards,

Axel

XanderBakker
Esri Esteemed Contributor

Indeed that will help... Thanks!

XanderBakker
Esri Esteemed Contributor

Made some changes correction to the code by simply skipping the geometry and after that I encountered a memory error, so I included some delete statements and that seems to work.

import arcpy

def main():
    # settings
    ras = r"C:\GeoNet\ZonalStatsMultiPart\test\test.gdb\test"
    fc = r"C:\GeoNet\ZonalStatsMultiPart\gwbodygeom_gk3_20140318Kopie.shp"
    lst_perc = [2, 5, 10, 25, 50, 75, 90, 95, 98] # list of percentiles to be calculated
    fld_prefix = "Perc_"
    start_at = 0

    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)
            if i >= start_at:
                polygon = row[flds.index("SHAPE@")]
                lst_parts = []
                for part in polygon:
                    for pnt in part:
                        try:
                            x, y = pnt.X, pnt.Y
                            lst_parts.append(arcpy.Point(x, y))
                        except:
                            pass

                # Execute ExtractByPolygon (you can't send the polygon object)
                print " - ExtractByPolygon..."
                ras_pol = arcpy.sa.ExtractByPolygon(ras, lst_parts, "INSIDE")
                del lst_parts, polygon
                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)}
                del ras_pol

                # 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
                del dct

                # 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
                del dct_per

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

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

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

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

if __name__ == '__main__':
    main()

See attached shapefile.

AxelThomas
New Contributor III

Hi Xander,

thanks for this one - worked perfect.

One additional note: I initially received an error:

RuntimeError: ERROR 010240: Raster-Dataset in_memory\ras61 could not be saved in format MEM      (my translation)

I assumed a memory problem on my pc and changed

arcpy.env.workspace = "IN_MEMORY"

to my local default gdb, e.g.

arcpy.env.workspace =r"W:\user\thomasa\Default.gdb"

This worked but then the code appears to run somewhat slower.

Thanks again, Xander, your code saves me quite some time!

Best regards,

Axel

XanderBakker
Esri Esteemed Contributor

That is a very good solution, but maybe you should include some code to clean up the gdb, since all temporal rasters will be stored there.

AxelThomas
New Contributor III

Right - I think of creating a new gdb at the begin of the code which will be deleted at the end of the code.

Thanks again, Xander

XanderBakker
Esri Esteemed Contributor

Either that or use the scratch workspace...

ZengYi-Chao
New Contributor

Hi Xander

I used this code that you issued on 2015/8/27 but i found the results of all  percentiles are the same value. Do you have any ideas what might cause the error?

0 Kudos
XanderBakker
Esri Esteemed Contributor

If you can post the data (or a part of it) then I can take a look at what might be causing this. 

ZengYi-Chao
New Contributor

Xander

Thanks for your help. The attached file includes raster(DEM) and shape file. Wish that can help you to solve this error.

Best regards

0 Kudos