Nikos A:
> > Is it about trimming Landsat borders?
Anna Z:
> No, it's about selecting an area inside the tile, which I need to work on.
> anyway, since you mentioned it, the issue of cutting away borders is also
> interesting ...
In short
Attached a custom python script (for GRASS GIS 7 I think) which needs a
Landsat scene to be imported based on a GRASS-Wiki Landsat import script [1]
and the corresponding vector tile to be imported in GRASS as well. If you get
into this I'll try to make it more clear how I used it.
Attached as well, a slightly modified import script, based on the GRASS-Wiki
one, to handle, hopefully all: L5, L7 and L8.
Long story
There is this question: <http://gis.stackexchange.com/q/34637/5256> and a
GRASSy answer: <http://gis.stackexchange.com/a/45268/5256>. So there is/was
already a grass-addon named i.landsat.trim in
<http://grasswiki.osgeo.org/wiki/AddOns/GRASS_6#i.landsat.trim>. However, it
is somewhat experimental and not using the exact WRS2 tiles.
Back in May I have tried to give a generic answer to this task by writing a
Python script called "i.landsat.wrs2trim" :-). I imagined that a script could
simply be instructed which tile(s) and row(s) from the official WRS2 tiles [0]
and it should trim all Landsat bands that pertain to an aquisition over the
tile(s) in question.
I ended up, once again, with something kinda hardcoded that would require to
have prepared and have beforehand the vector tile of interest. It is attached.
It works only for "LE7" scenes and only if the bands are imported using the
python script given in GRASS-Wiki [1] (thanks to Martin for this), i.e. they
are named after B1, B2, B3, B4, B5, B61, B62, B8.
Of course (!), one needs to prepare the vector tile for the specific to-be-
trimmed Landsat bands. As far as I remember, this involves getting the WRS2
Shapefile, selecting and saving the tile of interest (in QGIS quicker I think,
load > select > save as...), then importing in GRASS.
The logic is to have Landsat scenes imported in independent Mapsets, that is
all bands of a scene in one Mapset, and all bands of another scene in another
Mapset, et.c. Each Mapset named after the Landsat scenes' ID.
I didn't know how to have available the WRS2 tiles. Let the user download
them and set the path to it, and let the script do its work? Or, make the
script do all this? Packing the WRS2 tiles along with the script somehow?
Bah... Not a programmer here :-/
Oh, thanks to Alexander Muriy for his time at some point for checking the
python script and giving some hints. I just have had no time to work on this
again. Nonetheless, I consider such a tool as a must-have along the other
landsat tools in GRASS.
Nikos
--
[0] <http://landsat.usgs.gov/documents/wrs2_descending.zip>
[1] <http://grasswiki.osgeo.org/wiki/LANDSAT#Automated_data_import>
#!/usr/bin/python
# source: <http://grasswiki.osgeo.org/wiki/LANDSAT#Automated_data_import>
# ToDo
# # use grass.warning(_("...")?
# 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]) # or: os.path.splitext(file)[0].split('_')[-3] + os.path.splitext(file)[0].split('_')[-1]
else:
name = os.path.splitext(file)[0].split('_')[-1]
# band = int(name[1:2] if (len(name) == 3 and name[-1] == '0') else name[1:3] if (len(name) == 3 and name[-1] != '0') else name[-1:])
# found a wrongly named *MTL.TIF file in LE71610432005160ASN00
# skip it!
if ('MTL') in ffile: # use grass.warning(_("...")?
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()#!/usr/bin/env python
# Scripting started 15. 05. 2013
# Working state reached on 20. 05. 2013
############################################################################
#
# MODULE: i.landsat.wrs2trim
# AUTHOR(S): Nikos Alexandris
# PURPOSE: Trims all per-Mapset Landsat-Scene images based on a given WRS2
# REQUIREMENTS: The user should provide a corresponding, vector or raster,
# WRS2 tile to be used as a trimming mask. Obtain the official
# WRS2 Path/Row scene boundaries from:
# <http://landsat.usgs.gov/documents/wrs2_descending.zip>
# More information to be found at:
# <http://landsat.usgs.gov/tools_wrs-2_shapefile.php>
# Path/Row (raster) tile
# COPYRIGHT: (C) 2013 by the GRASS Development Team
#
# This program is free software under the GNU General Public
# License (>=v2). Read the file COPYING that comes with GRASS
# for details.
#
#############################################################suunto core yellow################
#%module
#% description: Trims all per-Mapset images that pertain to a Landsat scene
#% based on a corresponding WRS2 path/row (vector or raster) tile.
#% keywords: imagery
#% keywords: landsat
#% keywords: trimming
#%end
#%option
#% key: tile
#% description: Name of input WRS2 path/row (vector or raster) tile to trim with
#% required : yes
#%end
#%option
#% key: path
#% type: string
#% description: WRS2 Path of interest
#% required: no
#%end
#%option
#% key: row
#% type: string
#% description: WRS2 Row of interest
#% required : no
#%end
#%option
#% key: prefix
#% type: string
#% description: Prefix for trimmed images
#% answer: B.Trimmed
#% required : yes
#%end
# (1) import required libraries
#import atexit # not used yet!
import grass.script as grass
from grass.pygrass.gis import Mapset
from grass.pygrass.vector import VectorTopo
# (3) global variables
# Landsat Spectral bands
SPECTRAL_BANDS_NR = [1, 2, 3, 4, 5, 7]
SPECTRAL_BANDS_NAMES = ['B1', 'B2', 'B3', 'B4', 'B5', 'B7']
SPECTRAL_BANDS_NR_ALT = [10, 20, 30, 40, 50, 70] # for the time being...
SPECTRAL_BANDS_NAMES_ALT = ['B10', 'B20', 'B30', 'B40', 'B50', 'B70']
# Temperature channel
TEMPERATURE_CHANNELS_NR = [61, 62]
TEMPERATURE_CHANNELS_NAMES = ['B61', 'B62']
# Panchromatic channels
PANCHROMATIC_CHANNEL_NR = [8]
PANCHROMATIC_CHANNEL_NAME = ['B8']
PANCHROMATIC_CHANNEL_NR_ALT = [80]
PANCHROMATIC_CHANNEL_NAME_ALT = ['B80']
# (4) functions
#############################################################################
# cleaning up
def cleanup():
if rastertmp:
grass.run_command('g.remove',
rast = rastertmp,
quiet = True
)
grass.run_command('g.remove',
rast = 'MASK',
quiet = True,
stderr = nuldev
)
if mask_found:
grass.message(_("Restoring previous MASK..."))
grass.run_command('g.rename',
rast = (tmpname + "_origmask", 'MASK'),
quiet = True
)
## check if a MASK is already present ========================================
#if grass.find_file('MASK', mapset = Mapset)['file']:
#UserMASK = "usermask_mask." + unique # this came from another script!
#grass.message(_("A user raster mask (MASK) is present. Saving it..."))
#grass.run_command('g.rename',
#quiet = quiet,
#rast = ('MASK',UserMASK)
#)
# List & Resolution of Landsat Rasters ######################################
def retrieve_landsat_rasters(verbose = False):
landsat_elements = dict()
# retrieve existing Spectral bands --------------------------------------
landsat_elements['Spectral Bands'] = \
Mapset().glist('rast', pattern = 'B[123457]')
# alternative naming...
if not landsat_elements['Spectral Bands']:
landsat_elements['Spectral Bands'] = \
Mapset().glist('rast', pattern = 'B[123457]0')
# no results?
if not landsat_elements['Spectral Bands']:
print("No Spectral bands named after a B? or B?0 pattern found!")
# resolution of Spectral bands
landsat_elements['Resolution Spectral'] = \
int(grass.raster_info(landsat_elements['Spectral Bands'][0])['nsres'])
# retrieve existing Temperature channels --------------------------------
landsat_elements['Temperature Channels'] = \
Mapset().glist('rast', pattern = 'B6[12]')
# no results?
if not landsat_elements['Temperature Channels']:
print("There are no Landsat Temperature channels named B61, B62!")
# temperature channels (resolution = 30 or 60m)
landsat_elements['Resolution Temperature'] = \
int(grass.raster_info(landsat_elements['Temperature Channels'][0])['nsres'])
# retrieve existing Panchromatic channels -------------------------------
landsat_elements['Panchromatic Channel'] = \
Mapset().glist('rast', pattern = 'B8')
# alternative naming
if not landsat_elements['Panchromatic Channel']:
landsat_elements['Panchromatic Channel'] = \
Mapset().glist('rast', pattern = 'B80')
# no results?
if not landsat_elements['Panchromatic Channel']:
print("There is no Landsat Panchromatic channel named B8 or B80!")
# panchromatic channels (resolution = 15m)
landsat_elements['Resolution Panchromatic'] = \
int(grass.raster_info(landsat_elements['Panchromatic Channel'])['nsres'])
# verbosity -------------------------------------------------------------
if verbose:
for key in sorted(landsat_elements.keys()):
print("%s: %r") % (key, landsat_elements[key])
return landsat_elements
# check type of queried raster: spectral, temperature or panchromatic #######
def landsat_raster_type(landsat_rasters):
# S controllers: if existing Spectral bands names are "regular"
if any(X in landsat_rasters for X in SPECTRAL_BANDS_NAMES):
landsat_type = 'Sx'
#landsat_raster_name = "Spectral (Bx)"
elif any(X in landsat_rasters for X in SPECTRAL_BANDS_NAMES_ALT):
landsat_type = 'Sxx'
#landsat_raster_name = "Spectral (Bxx)"
# Txx controller
elif any(X in landsat_rasters for X in TEMPERATURE_CHANNELS_NAMES):
landsat_type = 'Txx'
#landsat_raster_name = "Temperature (B6[1,2])"
# P controllers
elif any(X in landsat_rasters for X in PANCHROMATIC_CHANNEL_NAME):
landsat_type = 'Px'
#landsat_raster_name = "Panchromatic (B8)"
elif any(X in landsat_rasters for X in PANCHROMATIC_CHANNEL_NAME_ALT):
landsat_type = 'Pxx'
#landsat_raster_name = "Panchromatic (B80)"
else:
print("Something went wrong... !")
return landsat_type # , landsat_raster_name
# define a trim(map) function ###############################################
def trim(landsat_rasters, landsat_elements):
# region settings -------------------------------------------------------
current_region_settings = grass.parse_command("g.region", flags='g')
resolution_ew = int(current_region_settings['ewres'])
resolution_ns = int(current_region_settings['nsres'])
# prefix for trimmed rasters --------------------------------------------
trimmed_prefix = options['prefix']
# loop over all raster maps ---------------------------------------------
for Landsat_Raster in landsat_rasters:
# set the Resolution
if landsat_raster_type(landsat_rasters) == 'Sx' or \
landsat_raster_type(landsat_rasters) == 'Sxx': # (30m)
landsat_raster_name = 'spectral band'
#print landsat_raster_name
if not resolution_ew == 30 or not resolution_ns == 30:
Resolution = landsat_elements['Resolution Spectral']
#grass.message(_("Resolution set to %s" % Resolution))
if resolution_ew == 30 and resolution_ns == 30:
Resolution = resolution_ns
elif landsat_raster_type(landsat_rasters) == 'Txx': # 30m or 60m
landsat_raster_name = 'temperature channel'
#print landsat_raster_name
Resolution = landsat_elements['Resolution Temperature']
#grass.message(_("Resolution set to %s" % Resolution))
elif landsat_raster_type(landsat_rasters) == 'Px' or \
landsat_raster_type(landsat_rasters) == 'Pxx': # (15m)
landsat_raster_name = 'panchromatic channel'
#print landsat_raster_name
Resolution = landsat_elements['Resolution Panchromatic']
#grass.message(_("Resolution set to %s" % Resolution))
else:
print("Can't identify a resolution value!")
# match region extent to PathRowTile and set resolution of interest
grass.run_command('g.region',
vect = wrs2_Tile,
res = Resolution,
flags = 'a',
verbose = False
)
# expression for grass.mapcalc
expression = "%s.%s = %s" % (trimmed_prefix,
Landsat_Raster[1:],
Landsat_Raster
)
# communicate
grass.message(_("Trimming %s %s (%im)" % (landsat_raster_name,
Landsat_Raster,
Resolution
)
)
)
# r.mapcalc-it
grass.mapcalc(expression,
quiet = False,
verbose = True,
overwrite = True
)
# the main function #########################################################
def main():
# global variables -- better in the beginning of the code?
global wrs2_Tile, \
SPECTRAL_BANDS_NAMES, SPECTRAL_BANDS_NAMES_ALT, \
TEMPERATURE_CHANNELS_NAMES, \
PANCHROMATIC_CHANNEL_NAME, PANCHROMATIC_CHANNEL_NAME_ALT
# Obtain from Options: Path & Row of WRS2 trimming tile -----------------
wrs2_Tile = options['tile']
# check if the tile in question exists
if not grass.find_file(wrs2_Tile, 'vector')['file']:
grass.fatal(_("Tile <%s> not found") % wrs2_Tile)
wrs2_Path = options['path']
wrs2_Row = options['row']
print "Path, Row are set to %s, %s" % (wrs2_Path, wrs2_Row)
#wrs2_PathRow = wrs2_Path + wrs2_Row
# better get it from the Mapset's name?
grass.message(_("Trimming Mask based on <%s>" % wrs2_Tile))
print ''
# obtain list with per-Landsat-Scene Mapsets & loop over them -----------
mapsets = grass.mapsets()
for scene_id in mapsets:
#print scene_id
# ToDo: How to make it work with *any* Landsat Scene ID!
# ToDo: Give option for manual definition of Mapsets to be processed
# -- check at i.histo.match code, for example?
if 'LE7' in scene_id:
# enter in (a) Mapset of interest
Mapset(scene_id).current()
print "\nOperating inside Scene-Mapset <%s>\n" % scene_id
# grant access to current and PERMANENT Mapsets only
grass.run_command("g.mapsets",
operation = 'set',
mapset= '.,PERMANENT',
verbose = False
)
# get WRS2 Path/Row from the Scene Identifier == name of Mapset
wrspr = scene_id[3:9]
# path
#PATH = scene_id[3:6]
# row
#ROW = scene_id[6:9]
# retrieve rasters
print "\n Retrieving information about existing Landsat data in \
Mapset <%s>..." % scene_id
landsat_elements = retrieve_landsat_rasters()
print "\n", landsat_elements
# query the tile's bounding box coordinates (pygrass)
with VectorTopo(wrs2_Tile) as Query_Vector_Tile:
North, South, East, West = Query_Vector_Tile.bbox().nsewtb()
# match region extent to the tile's bounding box coordinates
# Resolution = ?
grass.run_command('g.region',
n = North,
s = South,
e = East,
w = West,
res = 30,
flags = 'a',
verbose = False
)
# set PathRowTile as MASK -- why not use the vector tile?
grass.run_command('r.mask',
vector = wrs2_Tile,
where ='WRSPR = %s' % wrspr,
overwrite = True,
verbose = False
)
print ''
# trim'em
trim(landsat_elements['Spectral Bands'], landsat_elements)
trim(landsat_elements['Temperature Channels'], landsat_elements)
trim(landsat_elements['Panchromatic Channel'], landsat_elements)
# remove PathRowTile-MASK
grass.run_command('r.mask',
flags = 'r',
verbose = True
)
else:
print "This Location will not be processed!\n"
# finally, remove the raster MASK, set to previous region settings
if __name__ == "__main__":
options, flags = grass.parser()
main()
_______________________________________________
grass-user mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/grass-user
