LANDSAT

2012-09-11T08:10:23Z

⚠️Ljbade: Fixed bug in python code

== LANDSAT Data Availability==

* data download from [http://landsat.usgs.gov/ USGS LANDSAT archive] and [http://glcfapp.umiacs.umd.edu GLCF] and [http://edcsns17.cr.usgs.gov/NewEarthExplorer/ EarthExplorer] and [http://glovis.usgs.gov USGS GloVis]
* see also the [[Global datasets]] wiki page

== Modules overview ==

* {{cmd|d.rgb}} - display 3-band data
* {{cmd|i.landsat.rgb}} - auto-enhance colors
* {{cmd|i.atcorr}} - correct top of atmosphere to surface reflectance - see also the [[Atmospheric correction]] wiki page
* {{cmd|i.oif}} - calculate the 3 bands showing the greatest difference (for use as R,G,B bands)
* {{cmd|r.composite}} - flatten 3-bands of data into a single image (lossy)

* {{AddonCmd|i.landsat.toar}} (addon, included in GRASS 6.5 and 7) - convert DN to top of atmosphere radiance
* {{cmd|i.landsat.acca}} - cloud identification
* {{AddonCmd|i.landsat.dehaze}} (addon) - haze removal
* {{cmd|i.topo.corr}} -used to topographically correct reflectance from imagery files, e.g. obtained with {{cmd|i.landsat.toar}}, using a sun illumination terrain model.

== LANDSAT Pre-Processing ==

=== Import data ===

# Open GRASS GIS, select "Location Wizard" in order to create a new location from georeferenced file
# Use {{cmd|r.in.gdal}} or the "Import file tool" to import the GeoTIFF files into GRASS GIS. See also "Automated data import" below.
## Most Landsat scenes are delivered in "north-is-up" orientation, hence import is straightforward
## If you get "ERROR: Input map is rotated - cannot import.", the image must be first externally rotated, applying the rotation info stored in the metadata field of the raster image file. For example, the gdalwarp software can be used to transform the map to North-up (note, there are several gdalwarp parameters to select the resampling algorithm):
gdalwarp rotated.tif northup.tif

==== Automated data import ====

An example of Python script bellow imports raster data into GRASS. For each image creates its own mapset and imports bands as <tt>B<id></tt>, e.g. B10, B20, ..., B80. This script also sets up timestamp based on MTL file.

The {{AddonCmd|i.landsat.toar}} and {{AddonCmd|i.landsat.acca}} modules want the maps to be named such as "scenename.1, .2, .3", etc. for the different bands. GloVis names LANDSAT-7 like _B10, _B20, _B30, etc. Save the following script as "import_landsat.py" file:

<source lang=python>
#!/usr/bin/python

import os
import sys
import glob
import grass.script as grass

def get_timestamp(mapset):
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 'ACQUISITION_DATE' in line:
result['date'] = line.strip().split('=')[1].strip()
finally:
fd.close()

return result

def import_tifs(mapset):
meta = get_timestamp(mapset)
for file in os.listdir(mapset):
if os.path.splitext(file)[-1] != '.TIF':
continue
ffile = os.path.join(mapset, file)
name = os.path.splitext(file)[0].split('_')[-1]
band = int(name[-1])
grass.message('Importing %s -> %s@%s...' % (file, name, mapset))
grass.run_command('g.mapset',
flags = 'c',
mapset = mapset,
quiet = True,
stderr = open(os.devnull, 'w'))
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)))

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>

Example of usage:

* <tt>./import_landsat.py</tt> walk through current directory and import all found satellite images.
* <tt>./import_landsat.py LM41890261983200FFF03</tt> imports images only from given directory.

=== Hint: Minimal disk space copies ===

Here's a little trick with {{cmd|r.reclass}} to rename maps without touching the data or wasting disk space:

<source lang="bash">
BASE=L71074092_09220040924

for BAND in 10 20 30 40 50 61 70 80; do
BAND1st=`echo $BAND | sed -e 's/0$//'`
r.reclass in="${BASE}_B$BAND" out=$BASE.$BAND1st << EOF
* = *
EOF
done
</source>

=== Natural color composites ===

Solution A) use {{cmd|i.landsat.rgb}}

Solution B) equalize colors on each R,G,B band with

<source lang="bash">
r.colors -e map=band1 color=grey
</source>

Composite: then run {{cmd|r.composite}}

=== Reset color tables ===

<source lang="bash">
BASE=L71074092_09220040924

for map in `g.mlist pat="$BASE.[0-8]*"` ; do
r.colors $map color=grey255
done
</source>

=== Create an MASK to only show data where coverage exists for all bands ===

<source lang="bash">
BASE=L71074092_09220040924

g.region rast=$BASE.1
r.series in=`g.mlist pat="$BASE.[0-8]*" sep=,` -n out=$BASE.thresh method=threshold thresh=1
r.mapcalc "$BASE.mask = if(isnull($BASE.thresh))"
g.remove $BASE.thresh
</source>

=== Calculate top-of-atmosphere reflectance and band-6 temperature ===

Calculate top-of-atmosphere reflectance and band-6 temperature (-t only if MTL, not MET file) with {{AddonCmd|i.landsat.toar}}:
i.landsat.toar input_prefix=$BASE output_prefix=${BASE}_toar sensor=tm7 metfile=${BASE}_MTL.txt -t

Note: the resulting temperature map is in Kelvin:
# convert to degree Celsius
r.mapcalc "$BASE.temp_celsius = ${BASE}_toar.6 - 273.15"
r.info -r $BASE.temp_celsius

=== Haze removal ===

Simple haze removal can be done with {{AddonCmd|i.landsat.dehaze}}. i.landsat.dehaze applies a bandwise haze correction using tasscap4 (haze) and linear regression.

=== Atmospheric correction ===

See [[Atmospheric correction]]

=== Cloud identification ===

Identify clouds in the image with {{cmd|i.landsat.acca|version=65}}:
i.landsat.acca -f input_prefix=226_62_toar. output=226_62.acca_cloudmask

Mask out the clouds:
r.mapcalc "MASK = if(isnull($BASE.acca_cloudmask))"

== Download sample data ==

=== Preprocessed Landsat-7 data for North Carolina ===

The [[Sample_datasets|North Carolina 2008 sample dataset]] comes with 3 different Landsat scenes:
: (''Wake County -- path: 16 row: 35 for various dates'')

The above import efforts are not needed since the data are already in a GRASS location.

=== Landsat-5: Oct 14, 1987 ===

;Glovis download: LT50160351987287XXX08
: http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LT50160351987287XXX08&dataset_name=LANDSAT_TM



Provided metadata:
IMAGE_ID=P016R35_5T871014
PATH=16
ROW=35
DATE=10/14/87
PLATFORM=LANDSAT5

=== Landsat-7: Mar 31, 2000 ===

;Glovis download: LE70160352000091EDC00
: http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LE70160352000091EDC00&dataset_name=LANDSAT_ETM

Values in the metadata below indicate that the version provided with the NC 2008 dataset's production date was after July 1, 2000, and that the channel gains were <tt>HHHLHLHHL</tt>.

Convert DNs to radiance/temperatures:
GRASS> i.landsat.toar -v band=lsat7_2000 sensor=7 date=2000-03-31 \
product_date=2000-07-02 solar_elevation=51.5246529 gain=HHHLHLHHL

Provided metadata:
SPACECRAFT_ID=Landsat7
SENSOR_ID=ETM+
'''ACQUISITION_DATE=2000-03-31'''
WRS_PATH=16
CPF_FILE_NAME=L7CPF20000101_20000331_12
SUN_AZIMUTH=139.6033279
'''SUN_ELEVATION=51.5246529'''
LMAX_BAND1=191.600
LMIN_BAND1=-6.200
LMAX_BAND2=196.500
LMIN_BAND2=-6.400
LMAX_BAND3=152.900
LMIN_BAND3=-5.000
LMAX_BAND4=241.100
LMIN_BAND4=-5.100
LMAX_BAND5=31.060
LMIN_BAND5=-1.000
LMAX_BAND61=17.040
LMIN_BAND61=0.000
LMAX_BAND62=12.650
LMIN_BAND62=3.200
LMAX_BAND7=10.800
LMIN_BAND7=-0.350
LMAX_BAND8=243.100
LMIN_BAND8=-4.700
QCALMAX_BAND1=255.0
QCALMIN_BAND1=1.0
QCALMAX_BAND2=255.0
QCALMIN_BAND2=1.0
QCALMAX_BAND3=255.0
QCALMIN_BAND3=1.0
QCALMAX_BAND4=255.0
QCALMIN_BAND4=1.0
QCALMAX_BAND5=255.0
QCALMIN_BAND5=1.0
QCALMAX_BAND61=255.0
QCALMIN_BAND61=1.0
QCALMAX_BAND62=255.0
QCALMIN_BAND62=1.0
QCALMAX_BAND7=255.0
QCALMIN_BAND7=1.0
QCALMAX_BAND8=255.0
QCALMIN_BAND8=1.0

=== Landsat-7: May 24, 2002 ===

;Glovis download: LE70160352002144EDC00
: http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LE70160352002144EDC00&dataset_name=LANDSAT_ETM



Provided metadata: (`p016r035_7x20020524.met`)
SPACECRAFT_ID=Landsat7
SENSOR_ID=ETM+
'''ACQUISITION_DATE=2002-05-24'''
WRS_PATH=016
WRS_ROW=035
SUN_AZIMUTH=120.8810347
'''SUN_ELEVATION=64.7730999'''
QA_PERCENT_MISSING_DATA=0
CLOUD_COVER=0
CPF_FILE_NAME=L7CPF20020401_20020630_03
LMAX_BAND1=191.600
LMIN_BAND1=-6.200
LMAX_BAND2=196.500
LMIN_BAND2=-6.400
LMAX_BAND3=152.900
LMIN_BAND3=-5.000
LMAX_BAND4=241.100
LMIN_BAND4=-5.100
LMAX_BAND5=31.060
LMIN_BAND5=-1.000
LMAX_BAND61=17.040
LMIN_BAND61=0.000
LMAX_BAND62=12.650
LMIN_BAND62=3.200
LMAX_BAND7=10.800
LMIN_BAND7=-0.350
LMAX_BAND8=243.100
LMIN_BAND8=-4.700
QCALMAX_BAND1=255.0
QCALMIN_BAND1=1.0
QCALMAX_BAND2=255.0
QCALMIN_BAND2=1.0
QCALMAX_BAND3=255.0
QCALMIN_BAND3=1.0
QCALMAX_BAND4=255.0
QCALMIN_BAND4=1.0
QCALMAX_BAND5=255.0
QCALMIN_BAND5=1.0
QCALMAX_BAND61=255.0
QCALMIN_BAND61=1.0
QCALMAX_BAND62=255.0
QCALMIN_BAND62=1.0
QCALMAX_BAND7=255.0
QCALMIN_BAND7=1.0
QCALMAX_BAND8=255.0
QCALMIN_BAND8=1.0

== LANDSAT Image classification ==

See [[Image classification]]

[[Category:Documentation]]
[[Category:FAQ]]
[[Category:Image processing]]

GRASS-Wiki - User contributions [en]

LANDSAT