Latest Contributions by AndriesPotgieter

Dear Xander Bakker‌ I am trying to run the python script outside ArcGIS. It takes forever to run in ArcGIS so wanted to speed things up too. I am getting the following error when executing the first line from command line terminal in windows: A:\Andries_MainData\data\UAV\boxsite_2016_17\Data>C:\Python27\ArcGIS10.5\python.exe raster_percentile.py Traceback (most recent call last): File "raster_percentile.py", line 182, in <module> main() File "raster_percentile.py", line 23, in main arcpy.CreateTable_management(tbl_ws, tbl_name) File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\management.py", line 17411, in CreateTable raise e arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid. ERROR 000732: Table Location: Dataset A:/data/UAV/boxsite_2016_17/Data/Geodatabase.gdb does not exist or is not supported Failed to execute (CreateTable). Can you please advise what I am doing wrong? I will upload the script also regards andries ... View more

Hi Xander Bakker‌, I have run the Zonal percentile code you wrote on my data for 2 rasters and a snippet of the results is attached as an image here FYI. thanks again for your help in this. much appreciated. Andries  ... View more

Dear Xander Bakker‌ Thanks for the advice. I have changed the code accordingly and it is working now. Will keep you posted when it has finished. regards ap ... View more

Hi Xander Bakker‌ after creating a new geodatabase and copying all relevant rasters and shapefiels to it I now get this error: execfile(r'F:/data/UAV/Lodging_model/Geodatabase/Raster_percentile.py') Processing Raster 'after_raster' - Processing polygon: 1 - ExtractByPolygon... - saved raster as ras0pol1 - fill dict with Value x Count 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 67, in main dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(outname, flds_ras)} RuntimeError: cannot open 'ras0pol1' >>> Would be great if you can perhaps advise me what the likely problem is? cheers andries ... View more

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() ... View more

Hi Xander Bakker‌ is it possible to provide the code line for creating 1 single shapefile or table for the resulting zonal statistics here? thanks andires ... View more

 Magic!! Xander Bakker‌ I just multiplied the raster height by 1000 at your INt convert line ras_pol = arcpy.sa.Int(ras_pol1 * 1000) I would like to run it on 2 rasters and have all data in a table format for each plot. I assume I can copy the second raster to geodatabase and include the INT code transform line in it? much appreciated Andries ... View more

Hi Xander Bakker‌, Apologies for the late reply. The files are too large to upload here, unfortunately. here is the link that has the plot grid data shapefile and the raster. It would be good if you can make it run on more than 1 rater by maybe duplicating the raster. the files are copied to this dropbox link: Dropbox - data_Xander - Simplify your life  I hope you can find the problem. cheers andries ... View more

Hi Zander Thanks for your reply. I will upload the two files during the next two days or can i send you a dropbox link? Andries ... View more

Hi Xander Bakker‌ Thanks for your code on extracting percentiles from raster using polygon. It runs perfectly but all the percentiles have the same value. i.e. Processing polygon: 138 - ExtractByMask... - saved raster as ras138 - fill dict with Value x Count - determine sum - sum=7294 - create percentile dict - iterate perceniles - Perc_dec for is 0.02 - Perc for 2% is 1426 - Store value - Perc_dec for is 0.05 - Perc for 5% is 1426 - Store value - Perc_dec for is 0.1 - Perc for 10% is 1426 - Store value - Perc_dec for is 0.25 - Perc for 25% is 1426 - Store value - Perc_dec for is 0.5 - Perc for 50% is 1426 - Store value - Perc_dec for is 0.75 - Perc for 75% is 1426 - Store value - Perc_dec for is 0.9 - Perc for 90% is 1426 - Store value - Perc_dec for is 0.95 - Perc for 95% is 1426 - Store value - Perc_dec for is 0.98 - Perc for 98% is 1426 - Store value - update row I have followed your advice and used the latest code for extract mask. Still the same. My polygon is not projected and is 1.5mx3m and the raster pixels is 2.5cm.  Big data sets  hope you can advise what I should change to force the percentiles to be calculated uniquely? kind regards andries ... View more