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.

closed as off-topic by phloem, xunilk, whyzar, tinlyx, BERA May 16 '18 at 5:04

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "Questions relating to general IT or with no clear GIS component, are off-topic here but can be researched/asked at Stack Overflow (software development), Super User (computing hardware and software), Database Administrators (relational databases) and other SE sites" – phloem, xunilk, whyzar, tinlyx, BERA
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 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.

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