Hi Xander Bakker
I have adjusted the paths fro using your code of creating zonal percentiles on 2 rasters or more:
I am getting the following error I suspect it is a matter of wrong path somewhere/
can you please advise?
regards
andries
error:
Runtime error Traceback (most recent call last): File "<string>", line 1, in <module> File "F:/data/UAV/Lodging_model/Geodatabase/Raster_percentile.py", line 171, in <module> main() File "F:/data/UAV/Lodging_model/Geodatabase/Raster_percentile.py", line 59, in main ras_pol = arcpy.sa.ExtractByPolygon(ras, lst_parts, "INSIDE") File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Functions.py", line 1482, in ExtractByPolygon extraction_area) File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Utils.py", line 53, in swapper result = wrapper(*args, **kwargs) File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Functions.py", line 1477, in Wrapper extraction_area) File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\geoprocessing\_base.py", line 510, in <lambda> return lambda *args: val(*gp_fixargs(args, True)) ExecuteError: ERROR 999999: Error executing function. The table was not found. [VAT_after_raster] ERROR 010060: Could not open and write to c:\PROGRA~2\arcgis\DESKTO~1.5\bin\pol83. Failed to execute (ExtractByPolygon).
code:
import arcpy
import os
def main():
# settings
ras_ws = r"F:/data/UAV/Lodging_model/Percentile.gdb/" # input raster workspace
fc = r"F:/data/UAV/Lodging_model/Percentile.gdb/Hermitage_Lodging_Grid" # input polygons
tbl = r"F:/data/UAV/Lodging_model/Percentile.gdb/Zonal_Perc" # output related table
lst_perc = [2, 98] # list of percentiles to be calculated
fld_prefix = "Perc_"
fld_poloid = "POL_OID"
fld_rasname = "RasterName"
# create table
tbl_ws, tbl_name = os.path.split(tbl)
arcpy.CreateTable_management(tbl_ws, tbl_name)
# add fields and fill fields list
flds_tbl = [fld_poloid, fld_rasname]
arcpy.AddField_management(tbl, fld_poloid, "LONG")
arcpy.AddField_management(tbl, fld_rasname, "TEXT", 50)
for perc in lst_perc:
fld_perc = "{0}{1}".format(fld_prefix, perc)
arcpy.AddField_management(tbl, fld_perc, "LONG")
flds_tbl.append(fld_perc)
# Enable Spatial analyst
arcpy.CheckOutExtension("Spatial")
# get list of rasters
lst_ras = getListOfRasterFromWS(ras_ws)
# environments
arcpy.env.workspace = "F:/data/UAV/Lodging_model/Percentile.gdb/"
arcpy.env.overwriteOutput = True
# create dictionary with polygons OID vs lst_parts
dct_pols = createDictPolygons(fc, fld_poloid)
# start insert cursor for table
with arcpy.da.InsertCursor(tbl, flds_tbl) as curs_tbl:
# start loop through rasters
i = 0
for ras_name in lst_ras:
ras = os.path.join(ras_ws, ras_name)
i += 0
print "Processing Raster '{0}'".format(ras_name)
# loop through polygons
for oid, lst_parts in dct_pols.items():
lst_row = [oid, ras_name]
print " - Processing polygon: {0}".format(oid)
# Execute ExtractByPolygon (you can't send the polygon object)
print " - ExtractByPolygon..."
ras_pol = arcpy.sa.ExtractByPolygon(ras, lst_parts, "INSIDE")
outname = "ras{0}pol{1}".format(i, oid)
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 percentiles"
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 in list"
# fld_perc = "{0}{1}".format(fld_prefix, perc)
# row[flds.index(fld_perc)] = pixval
lst_row.append(pixval)
# update row
print " - insert row"
row = tuple(lst_row)
curs_tbl.insertRow(row)
# return SA license
arcpy.CheckInExtension("Spatial")
print "Ready..."
def createDictPolygons(fc, fld_poloid):
flds = ("OID@", "SHAPE@", fld_poloid)
# add pol oid field for link
if not FieldExist(fc, fld_poloid):
arcpy.AddField_management(fc, fld_poloid, "LONG")
dct_pols = {}
with arcpy.da.UpdateCursor(fc, flds) as curs:
for row in curs:
oid = row[0]
polygon = row[1]
row[2] = oid
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)
# add to dict and store pol oid
dct_pols[oid] = lst_parts
curs.updateRow(row)
return dct_pols
def getListOfRasterFromWS(ras_ws):
arcpy.env.workspace = ras_ws
return arcpy.ListRasters()
def GetPixelValueForPercentile(dctper, percentile):
"""will return last pixel value
where percentile LE searched percentile."""
try:
srt_keys = sorted(dctper.keys())
pix_val = srt_keys[0]
for k in srt_keys:
perc = dctper
if perc <= percentile:
pix_val = k
else:
break
return pix_val
except Exception as e:
print " - GetPixelValueForPercentile error: {0}".format(e)
print "dctper:\n{0}".format(dctper)
return -9999
def FieldExist(featureclass, fieldname):
"""Check if field exists"""
import arcpy
fieldList = arcpy.ListFields(featureclass, fieldname)
return len(fieldList) == 1
if __name__ == '__main__':
main()
