User:NikosA/Landsat

From GRASS-Wiki
Revision as of 21:37, 30 June 2013 by ⚠️NikosA (talk | contribs) (Modified existing Landsat import script to take care about Landsat8)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.

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()