Python script for raster calculation using GDAL

Question

I am trying to loop over a raster calculation using GDAL. However, the loop appears not to work and QGIS crashes (without the loop the code which I found somewhere here works). The problem appears to be in the lines where I define input or output file. Any idea what could be the problem?

s = ['F101992','F101993','F101994','F121994','F121995','F121996','F121997','F121998','F121999','F141997','F141998','F141999','F142000','F142001','F142002','F142003','F152000','F152001','F152002','F152003','F152004','F152005','F152006','F152007','F162004','F162005','F162006','F162007','F162008','F162009','F182010','F182011','F182012']
for i in range(len(s)):
nighttime = "{mypath}//"+s(i)+".tif"
GDALnighttime = gdal.Open(nighttime, GA_ReadOnly )
band1 = GDALnighttime.GetRasterBand(1) 
Rdata = BandReadAsArray(band1)
dataOut = Rdata-((Rdata>63)*(Rdata-63))-((Rdata<=5)*Rdata)
outFile = "{mypath}//"+s(i)+"_ic_1.tif"
driver = gdal.GetDriverByName("GTiff") 
dsOut = driver.Create(outFile, GDALnighttime.RasterXSize, GDALnighttime.RasterYSize, 1, gdal.GDT_Float32) 
CopyDatasetInfo(GDALnighttime,dsOut) 
bandOut=dsOut.GetRasterBand(1) 
BandWriteArray(bandOut, dataOut)
dsOut =None
bandOut = None
driver = None

Answer 1

I'm not sure where your script is failing, although you seem to be going about saving the array in a convoluted fashion. It will be easier to use gdal_array.SaveArray, which takes an argument to specify a prototype file to get the projection info from.

from osgeo import gdal, gdal_array
for ext in s:
 # Get file names
 night_img = "{mypath}//"+ext+".tif"
 output = "{mypath}//"+ext+"_ic_1.tif"
 # Open band 1 as array
 ds = gdal.Open(night_img)
 b1 = ds.GetRasterBand(1)
 arr = b1.ReadAsArray()
 # apply equation
 data = arr-((arr>63)*(arr-63))-((arr<=5)*arr)
 # save array, using ds as a prototype
 gdal_array.SaveArray(data.astype("float32"), output, "GTIFF", ds)
 ds = None

• Hannes82: Thanks! Now it works. The problem was how I wrote the loop (range(len(s))...+s(i)+).