User:NikosA/Landsat

Processing Landsat scenes

Importing Landsat scenes in independent Mapsets

A modification of the script found at LANDSAT#Automated_data_import to take care for newer Landsat 8 scenes.

#!/usr/bin/python

# source: <http://grasswiki.osgeo.org/wiki/LANDSAT#Automated_data_import>

# imports
import os
import shutil
import sys
import glob
import grass.script as grass

def copy_metafile(mapset):
    """!Copies the *MTL.txt metadata file in the cell_misc directory inside the
    Landsat scene's independent Mapset
    """

    # get the metadata file
    try:
        metafile = glob.glob(mapset + '/*MTL.txt')[0]
        print '\nThe identified metadata file is <%s>.' % metafile.split('/')[1]

    except IndexError:
        return

    # get environment variables & define path to "cell_misc"
    gisenv=grass.gisenv()
    CELL_MISC_DIR = gisenv['GISDBASE'] + '/' + gisenv['LOCATION_NAME'] + '/' + gisenv['MAPSET'] + '/cell_misc'
    print 'It will be copied at: <%s>.\n' % CELL_MISC_DIR # better check if really copied!

    # copy the metadata file
    shutil.copy (metafile, CELL_MISC_DIR)

def get_timestamp(mapset):
    """!Gets the timestamp for each band of a Landsat scene from its metadata

    Returns the date of acquisition for each band of a Landsat scene from 
    its metadata file (*MTL.txt)
    """
    try:
        metafile = glob.glob(mapset + '/*MTL.txt')[0]

    except IndexError:
        return

    result = dict()

    try:
        fd = open(metafile)
        for line in fd.readlines():
            line = line.rstrip('\n')
            if len(line) == 0:
                continue
            if any(x in line for x in ('DATE_ACQUIRED', 'ACQUISITION_DATE')):
                result['date'] = line.strip().split('=')[1].strip()
    finally:
        fd.close()
 
    return result
 
def import_tifs(mapset):
    """!Imports all bands (TIF format) of a Landsat scene

    Imports all bands (TIF format) of a Landsat scene, be it Landsat 5, 7 or 8.
    All known naming conventions are respected, such as "VCID" and "QA" found
    in newer Landsat scenes for temperature channels and Quality... respectively.
    """
  
    #grass.message ('Importing... \n') # why is it printed twice?
    
    meta = get_timestamp(mapset)
    
    for file in os.listdir(mapset):
      
        if os.path.splitext(file)[-1] != '.TIF':
            continue
        
        ffile = os.path.join(mapset, file)
        # if correctly handled below, usr the "any" instruction!
        if ('QA') in ffile:
            name = "".join((os.path.splitext(file)[0].split('_'))[1::2])
        elif ('VCID') in ffile:
            name = "".join((os.path.splitext(file)[0].split('_'))[1::2])
        else:
            name = os.path.splitext(file)[0].split('_')[-1]
        
        # found a wrongly named *MTL.TIF file in LE71610432005160ASN00 > issue Warning!
        # really break the import process?
        if ('MTL') in ffile: 
            print("Found a wrongly named *MTL.TIF file!")
            print("Please rename the extension to .txt")
            break
        elif ('QA') in ffile:
            band = name
        elif len(name) == 3 and name[-1] == '0':
            band = int(name[1:2])
        elif len(name) == 3 and name[-1] != '0':
            band = int(name[1:3])
        else:
            band = int(name[-1:])
        
        grass.message('%s > %s@%s...' % (file, name, mapset))

        # create mapset of interest
        grass.run_command('g.mapset',
                          flags = 'c',
                          mapset = mapset,
                          quiet = True,
                          stderr = open(os.devnull, 'w'))

        # import bands
        if isinstance(band, str):
            grass.run_command('r.in.gdal',
                              input = ffile,
                              output = name,
                              quiet = True,
                              overwrite = True,
                              title = 'band %s' % band)
        else:
            grass.run_command('r.in.gdal',
                              input = ffile,
                              output = name,
                              quiet = True,
                              overwrite = True,
                              title = 'band %d' % band)

        if meta:
            year, month, day = meta['date'].split('-')
            if month == '01':
                month = 'jan'
            elif month == '02':
                month = 'feb'
            elif month == '03':
                month = 'mar'
            elif month == '04':
                month = 'apr'
            elif month == '05':
                month = 'may'
            elif month == '06':
                month = 'jun'
            elif month == '07':
                month = 'jul'
            elif month == '08':
                month = 'aug'
            elif month == '09':
                month = 'sep'
            elif month == '10':
                month = 'oct'
            elif month == '11':
                month = 'nov'
            elif month == '12':
                month = 'dec'
 
            grass.run_command('r.timestamp',
                              map = name,
                              date = ' '.join((day, month, year)))
    # copy MTL file
    copy_metafile(mapset)
 
def main():
    if len(sys.argv) == 1:
        for directory in filter(os.path.isdir, os.listdir(os.getcwd())):
            import_tifs(directory)
    else:
        import_tifs(sys.argv[1])
    
 
if __name__ == "__main__":
    main()