[GIS] Python script for raster calculation using GDAL

gdalpyqgisraster-calculator

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

Best Answer

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
Related Question