Accepted Solutions
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
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?
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
