2

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
3

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
  • Thanks! Now it works. The problem was how I wrote the loop (range(len(s))...+s(i)+). – Hannes82 Nov 7 '16 at 12:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.