Hi everyone,
I'm running a script that it's basically zonal statistic over polylines (I modified one script that I found in this forum Zonal statistics - percentiles ).
The code takes a shape file containing the lines along which I want to analyse the elevation. For every line, it creates a raster (let's say it converts the line into a raster) where in every pixel it writes the corresponding elevation value. It analyses the values of the elevation along the line and returns the max, min and percentiles values. It repeats the same thing for all the lines in the shape file.
The code is running without problems with a certain digital elevation model. Recently I received a new LiDAR file and I tried to run the code but I got an error:
RuntimeError: cannot open "line1.tif"
line1.tif is the line converted to raster. I checked the file and it is exactly the line in tif format with the right values of elevation. The folder in which it is stored, is the same as using the previous DEM.
The DEM and LiDAR have different resolution and therefore line1.tif in case of DEM in something like 2 KB and in case of LIDAR around 60 KB. This is the only difference I can notice.
Anyone has a clue of what is happening?
Thanks a lot
Fabio
Solved! Go to Solution.
mmm, I see what might be going wrong. You are using a raster with floating values. The script uses a search cursor on the attribute table of the raster (see line 82). Since the raster is float, it does not have an attribute table (hence the error).
In this part of the thread it is explained that the code requires an integer raster Re: Zonal statistics - percentiles
If the precision is not too exact want you require, you could convert the raster to integer and try again.
Xander Bakker this one is yours
Thanks Dan, I will try and see if this one can be resolved...
Hi, Fabio Cian ,
Did you make any other changes to the code and if so, could you include the code in this thread? (see: Posting Code blocks in the new GeoNet)
As far as I can deduce from the original code I wrote, it uses polygons to extract the percentiles. To do that, it will create a list of points and use that in the ExtractByPolygon tool to extract the cells inside the polygon. If your geometry is a straight line it probably will create an invalid polygon. If the line has a curve, it will probably include much more than you want in your data. To avoid this, I suppose you could create a a buffer around the lines (using the cellsize as buffer width) and use those polygons in the script.
Is it possible to attach the line1.tif to this thread to have a look of what you are working with, or include a screen shot of the data?
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 = dctperif 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()
mmm, I see what might be going wrong. You are using a raster with floating values. The script uses a search cursor on the attribute table of the raster (see line 82). Since the raster is float, it does not have an attribute table (hence the error).
In this part of the thread it is explained that the code requires an integer raster Re: Zonal statistics - percentiles
If the precision is not too exact want you require, you could convert the raster to integer and try again.
Thanks a lot! That was the problem!
Fabio
Glad it solved your problem...