Hello, thanks a lot for replying.
Attached the file (named ras1.tif instead of line1.tif). This is the code:
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"""
import arcpy
fieldList = arcpy.ListFields(featureclass, fieldname)
return len(fieldList) == 1
import arcpy, sys
from arcpy import env
from arcpy.sa import *
env.workspace = r"C:\Users\Flood Depth"
f = open('log_file.txt', 'w')
# Enable Spatial analyst
arcpy.CheckOutExtension("Spatial")
# settings
ras = r"C:\Users\Dataset\Flood_Depth\LiD_MA_Sal_cm.tif"
fc = r"C:\Users\Dataset\Flood_Depth\C_Saletto.shp"
# 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_"
f.write("DEM File input: " + ras + "\n")
f.write("Shapefile input: " + fc + "\n" + "\n")
# 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)
# 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)
f.write("Processing polygon: {0}".format(i) + "\n")
polygon = row[flds.index("SHAPE@")]
# Execute ExtractByMask (you can't send the polygon object)
try:
print " - ExtractByMask..."
ras_pol = ExtractByMask(ras, polygon)
outname = "ras{0}".format(i) + ".tif"
ras_pol.save(outname)
print " - saved raster as {0}".format(outname)
except:
print arcpy.GetMessages(2) #"Unexpected error:", sys.exc_info()[0]
# 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] = float(cnt_i) / float(cnt_sum)
sort_dct_per = sorted(dct_per)
#print sort_dct_per
print dct_per
# loop through list of percentiles
print " - iterate perceniles"
f.write(" - iterate perceniles" + "\n")
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)
f.write(" - Perc_dec for is {0}".format(perc_dec) + "\n")
pixval = GetPixelValueForPercentile(dct_per, perc_dec)
print " - Perc for {0}% is {1}".format(perc, pixval)
f.write(" - Perc for {0}% is {1}".format(perc, pixval)+ "\n")
# 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)
del(outname)
# return SA license
arcpy.CheckInExtension("Spatial")
f.close()
