topic Re: Script Zonal statistics: runtime error - cannot open file in Python Questions

Script Zonal statistics: runtime error - cannot open file

FabioCian — Mon, 06 Jul 2015 14:39:26 GMT

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

Re: Script Zonal statistics: runtime error - cannot open file

DanPatterson_Retired — Mon, 06 Jul 2015 15:27:06 GMT

Xander Bakker this one is yours

Re: Script Zonal statistics: runtime error - cannot open file

XanderBakker — Mon, 06 Jul 2015 15:59:06 GMT

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?

Re: Script Zonal statistics: runtime error - cannot open file

FabioCian — Sat, 11 Dec 2021 13:15:06 GMT

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()

Re: Script Zonal statistics: runtime error - cannot open file

XanderBakker — Mon, 06 Jul 2015 18:13:54 GMT

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.

Re: Script Zonal statistics: runtime error - cannot open file

FabioCian — Tue, 07 Jul 2015 08:47:53 GMT

Thanks a lot! That was the problem!

Fabio

Re: Script Zonal statistics: runtime error - cannot open file

XanderBakker — Tue, 07 Jul 2015 20:38:29 GMT

Glad it solved your problem...