LANDSAT

From GRASS-Wiki
Revision as of 15:16, 10 September 2010 by Neteler (talk | contribs) (more hints)
Jump to navigation Jump to search

LANDSAT Data Availability

Modules

  • d.rgb - display 3-band data
  • i.landsat.rgb - auto-enhance colors
  • i.atcorr - correct top of atmosphere to surface reflectance
  • i.oif - calculate the 3 bands showing the greatest difference (for use as R,G,B bands)
  • r.composite - flatten 3-bands of data into a single image (lossy)

LANDSAT Pre-Processing

Some notes from Yann Chemin:

  1. Open GRASS GIS, select create location from georeferenced file
  2. Use r.in.gdal to import into GRASS GIS.
    1. Most landsat scenes are delivered in "north-is-up" orientation, hence import is straightforward
    2. 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
  1. In grass/add-ons look for any of these to correct from radiance to reflectance at top of atmosphere
    1. i.landsat.toar
    2. i.dn2full.l5
    3. i.dn2full.l7
  2. Use i.atcorr to correct top of atmosphere to surface reflectance

Notes

The i.landsat.toar and 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.

Here's a little trick with r.reclass to rename it without touching the data or wasting disk space:

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

Calculate top-of-atmosphere reflectance:

i.landsat.toar -7 band_pre=$BASE metfile=L71074092_09220040924_MTL.txt

Identify clouds in the image:

i.landsat.acca -2 -f band_prefix=$BASE.toar out=$BASE.acca

Mask out the clouds:

r.mapcalc "MASK = if( isnull($BASE.acca), 1, null() )"