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.
