⚠️NikosA: Modified existing Landsat import script to take care about Landsat8

2013-06-30T21:37:00Z

Modified existing Landsat import script to take care about Landsat8

New page

= Processing Landsat scenes =

== Importing Landsat scenes in independent Mapsets ==

A modification of the script found at [http://grasswiki.osgeo.org/wiki/LANDSAT#Automated_data_import LANDSAT#Automated_data_import] to take care for newer Landsat 8 scenes.

<source lang="python">
#!/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()
</source>

User:NikosA/Landsat - Revision history

⚠️NikosA: Modified existing Landsat import script to take care about Landsat8