0

I have a list of landsat 8 rasters in tif format and i want to calculate an index with python for each escene, using red and NIR bands.

The individual (simplified) script is:

from osgeo import gdal
from osgeo.gdalnumeric import *
from osgeo.gdalconst import *
import numpy as np
import os
import glob
import sys

L8_red = 'Z:/bands/LC08_L1TP_233087_20180325_20180404_01_T1_sr_band4.tif'
L8_nir = 'Z:/bands/LC08_L1TP_233087_20180325_20180404_01_T1_sr_band5.tif'

outFile = "out.tif"

# open
redds = gdal.Open(L8_red)
nirds = gdal.Open(L8_nir)

redband = redds.GetRasterBand(1)
nirband = nirds.GetRasterBand(1)

red = redband.ReadAsArray()
nir = nirband.ReadAsArray()

red = red.astype(np.float64)
nir = nir.astype(np.float64)

#The actual calculation
def SR_pqa(red, nir):
    SR = (nir / red)
    return SR

SR = SR_pqa(red, nir)

#Write the out file
driver = gdal.GetDriverByName("GTiff")
dsOut = driver.Create("out.tif", redds.RasterXSize, redds.RasterYSize, 1, gdal.GDT_Float32)
CopyDatasetInfo(redds,dsOut)
bandOut=dsOut.GetRasterBand(1)
BandWriteArray(bandOut, SR)

I need to do the same with all the scenes in the folder.

They match until the date (LC08_L1TP_233087_*_01_T1_sr_bandX.tif).

There's 50 different dates in the path, so is very slow to do individualy in QGIS raster calculator.

  • So what's the question? Aren't landsat bands 16 bit integer? Why convert to float64? Shouldn't you test in SR_pqa if red==0: return 0 else: do the calc? otherwise you have a division by 0 error. – Michael Stimson May 15 '18 at 22:40
  • 4
    If I understand correctly, your question boils down to "how to iterate through each file (or similarly named files) in a folder and do something to each," which is a pretty general coding problem, best researched on a more general site like Stack Overflow. – phloem May 15 '18 at 22:44
  • See: stackoverflow.com/questions/19587118/… – BERA May 16 '18 at 5:04
  • with int16 the output values results without decimals. Your right about the division by 0 error, but i simplified the code for the question, so these pixels are no data. Thanks for the comments, no doubt i should edit. – OSCAR_P May 16 '18 at 16:16
1

You could try something like this:

# assuming all files have the format LC08_L1TP_233087_20180325_20180404_01_T1_sr_band4.tif
import os, sys
dates   = []                        # empty lists to start
Rasters = []

for thisFile in os.listdir(r'Z:\bands'):
    fN, fE = os.path.splitext(thisFile) # break into name and extension
    if fE.upper() == '.TIF':
        Rasters.append(thisFile)

for ThisRaster in Rasters:
    RasterNameSplit = ThisRaster.split('_')
    if len(RasterNameSplit) > 5: # omit non-conforming names
        ThisDateString  = '_'.join(RasterNameSplit[2:5]) # join the 3rd to 5th element like 233087_20180325_20180404
        if not ThisDateString in dates:
            dates.append(ThisDateString)

# validate the dates
ValidDates=[]
for ThisDateString in dates:
    L8_red = 'LC08_L1TP_{}_01_T1_sr_band4.tif'.format(ThisDateString)
    L8_nir = 'LC08_L1TP_{}_01_T1_sr_band5.tif'.format(ThisDateString)
    if os.path.exists(os.path.join(r'Z:\bands',L8_red)):
        if os.path.exists(os.path.join(r'Z:\bands',L8_nir)):
            ValidDates.append(ThisDateString) # both red and nir bands exist

for ThisDateString in ValidDates:
    L8_red = 'LC08_L1TP_{}_01_T1_sr_band4.tif'.format(ThisDateString)
    L8_nir = 'LC08_L1TP_{}_01_T1_sr_band5.tif'.format(ThisDateString)
    outFile = "out_{}.tif".format(ThisDateString) # make the name unique

    # open
    redds = gdal.Open(os.path.join(r'Z:\bands',L8_red))
    nirds = gdal.Open(os.path.join(r'Z:\bands',L8_nir))
    # rest of your script

This uses a few steps:

  1. List all the files in your folder, find the ones that end with '.tif' and add to a list.
  2. Go through the list of TIFF files and find the date string matching, add each unique value to the list.
  3. Iterate the date strings and find the ones that have a band4 and band5 existing appending to a list of valid rasters.
  4. Iterate the valid rasters to call your script. I have changed the name of the output file to include the date string otherwise it would be overwritten.
| improve this answer | |

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