#!/bin/sh # ############################################################################ # # MODULE: i.tasscap # AUTHOR(S): Markus Neteler. neteler itc.it # PURPOSE: At-satellite reflectance based tasseled cap transformation. # COPYRIGHT: (C) 1997-2004 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. # TODO: Check if MODIS Tasseled Cap makes sense to be added # http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1025776 ############################################################################# # References: # LANDSAT-4/LANDSAT-5: # script based on i.tasscap.tm4 from Dr. Agustin Lobo - alobo@ija.csic.es # TC-factor changed to CRIST et al. 1986, p.1467 (Markus Neteler 1/99) # Proc. IGARSS 1986 # # LANDSAT-7: # TASSCAP factors cited from: # DERIVATION OF A TASSELED CAP TRANSFORMATION BASED ON LANDSAT 7 AT-SATELLITE REFLECTANCE # Chengquan Huang, Bruce Wylie, Limin Yang, Collin Homer and Gregory Zylstra Raytheon ITSS, # USGS EROS Data Center Sioux Falls, SD 57198, USA # http://landcover.usgs.gov/pdf/tasseled.pdf # # This is published as well in INT. J. OF RS, 2002, VOL 23, NO. 8, 1741-1748. # Compare discussion: # http://adis.cesnet.cz/cgi-bin/lwgate/IMAGRS-L/archives/imagrs-l.log0211/date/article-14.html ############################################################################# # #%Module #% description: Tasseled Cap (Kauth Thomas) transformation for LANDSAT-TM data #% keywords: raster, imagery #%End #%flag #% key: 4 #% description: use transformation rules for LANDSAT-4 #%END #%flag #% key: 5 #% description: use transformation rules for LANDSAT-5 #%END #%flag #% key: 7 #% description: use transformation rules for LANDSAT-7 #%END #%option #% key: band1 #% type: string #% gisprompt: old,cell,raster #% description: raster input map (LANDSAT channel 1) #% required : yes #%end #%option #% key: band2 #% type: string #% gisprompt: old,cell,raster #% description: raster input map (LANDSAT channel 2) #% required : yes #%end #%option #% key: band3 #% type: string #% gisprompt: old,cell,raster #% description: raster input map (LANDSAT channel 3) #% required : yes #%end #%option #% key: band4 #% type: string #% gisprompt: old,cell,raster #% description: raster input map (LANDSAT channel 4) #% required : yes #%end #%option #% key: band5 #% type: string #% gisprompt: old,cell,raster #% description: raster input map (LANDSAT channel 5) #% required : yes #%end #%option #% key: band7 #% type: string #% gisprompt: old,cell,raster #% description: raster input map (LANDSAT channel 7) #% required : yes #%end #%option #% key: outprefix #% type: string #% gisprompt: new,cell,raster #% description: raster output TC maps prefix #% required : yes #%end #%flag #% key: s #% description: Process bands serially (default: run in parallel) #%end if [ -z "$GISBASE" ] ; then echo "You must be in GRASS GIS to run this program." >&2 exit 1 fi if [ "$1" != "@ARGS_PARSED@" ] ; then exec g.parser "$0" "$@" fi if [ $GIS_FLAG_4 -eq 1 ] ; then g.message "LANDSAT-4..." g.message "Calculating first TC component $GIS_OPT_OUTPREFIX.1 ..." r.mapcalc "\"$GIS_OPT_OUTPREFIX.1\" = 0.3037 * \"$GIS_OPT_BAND1\" \ + 0.2793 * \"$GIS_OPT_BAND2\" + 0.4743 * \"$GIS_OPT_BAND3\" \ + 0.5585 * \"$GIS_OPT_BAND4\" + 0.5082 * \"$GIS_OPT_BAND5\" \ + 0.1863 * \"$GIS_OPT_BAND7\"" & if [ "$GIS_FLAG_S" -eq 1 ] ; then # An alternate approach would be to add a BG='&' enviro var then # eval r.mapcalc "$CMD" $BG # but then the $CMD string would need triple quoting vs. "just" the # double quoting it needs now, and that would make the code even harder # to read. Backgrounding then waiting for processes incurs a little more # overhead, but for just three or four jobs the cumulative cost is small. # Life's a compromise. wait fi g.message "Calculating second TC component $GIS_OPT_OUTPREFIX.2 ..." r.mapcalc "\"$GIS_OPT_OUTPREFIX.2\" = -0.2848 * \"$GIS_OPT_BAND1\" \ - 0.2435 * \"$GIS_OPT_BAND2\" - 0.5435 * \"$GIS_OPT_BAND3\" \ + 0.7243 * \"$GIS_OPT_BAND4\" + 0.0840 * \"$GIS_OPT_BAND5\" \ - 0.1800 * \"$GIS_OPT_BAND7\"" & if [ "$GIS_FLAG_S" -eq 1 ] ; then wait fi g.message "Calculating third TC component $GIS_OPT_OUTPREFIX.3 ..." r.mapcalc "\"$GIS_OPT_OUTPREFIX.3\" = 0.1509 * \"$GIS_OPT_BAND1\" \ + 0.1973 * \"$GIS_OPT_BAND2\" + 0.3279 * \"$GIS_OPT_BAND3\" \ + 0.3406 * \"$GIS_OPT_BAND4\" - 0.7112 * \"$GIS_OPT_BAND5\" \ - 0.4572 * \"$GIS_OPT_BAND7\"" & wait elif [ $GIS_FLAG_5 -eq 1 ] ; then g.message "LANDSAT-5..." g.message "Calculating first TC component $GIS_OPT_OUTPREFIX.1 (Brightness) ..." r.mapcalc "\"$GIS_OPT_OUTPREFIX.1\" = 0.2909 * \"$GIS_OPT_BAND1\" \ + 0.2493 * \"$GIS_OPT_BAND2\" + 0.4806 * \"$GIS_OPT_BAND3\" \ + 0.5568 * \"$GIS_OPT_BAND4\" + 0.4438 * \"$GIS_OPT_BAND5\" \ + 0.1706 * \"$GIS_OPT_BAND7\" + 10.3695" & if [ "$GIS_FLAG_S" -eq 1 ] ; then wait fi g.message "Calculating second TC component $GIS_OPT_OUTPREFIX.2 (Greenness) ..." r.mapcalc "\"$GIS_OPT_OUTPREFIX.2\" = -0.2728 * \"$GIS_OPT_BAND1\" \ - 0.2174 * \"$GIS_OPT_BAND2\" - 0.5508 * \"$GIS_OPT_BAND3\" \ + 0.7221 * \"$GIS_OPT_BAND4\" + 0.0733 * \"$GIS_OPT_BAND5\" \ - 0.1648 * \"$GIS_OPT_BAND7\" - 0.7310" & if [ "$GIS_FLAG_S" -eq 1 ] ; then wait fi g.message "Calculating third TC component $GIS_OPT_OUTPREFIX.3 (Wetness) ..." r.mapcalc "\"$GIS_OPT_OUTPREFIX.3\" = 0.1446 * \"$GIS_OPT_BAND1\" \ + 0.1761 * \"$GIS_OPT_BAND2\" + 0.3322 * \"$GIS_OPT_BAND3\" \ + 0.3396 * \"$GIS_OPT_BAND4\" - 0.6210 * \"$GIS_OPT_BAND5\" \ - 0.4186 * \"$GIS_OPT_BAND7\" - 3.3828" & if [ "$GIS_FLAG_S" -eq 1 ] ; then wait fi g.message "Calculating fourth TC component $GIS_OPT_OUTPREFIX.4. (Haze) ..." r.mapcalc "\"$GIS_OPT_OUTPREFIX.4\" = 0.8461 * \"$GIS_OPT_BAND1\" \ - 0.0731 * \"$GIS_OPT_BAND2\" - 0.4640 * \"$GIS_OPT_BAND3\" \ - 0.0032 * \"$GIS_OPT_BAND4\" - 0.0492 * \"$GIS_OPT_BAND5\" \ - 0.0119 * \"$GIS_OPT_BAND7\" + 0.7879" & wait elif [ $GIS_FLAG_7 -eq 1 ] ; then g.message "LANDSAT-7..." g.message "Calculating first TC component $GIS_OPT_OUTPREFIX.1 (Brightness) ..." r.mapcalc "\"$GIS_OPT_OUTPREFIX.1\" = 0.3561 * \"$GIS_OPT_BAND1\" \ + 0.3972 * \"$GIS_OPT_BAND2\" + 0.3904 * \"$GIS_OPT_BAND3\" \ + 0.6966 * \"$GIS_OPT_BAND4\" + 0.2286 * \"$GIS_OPT_BAND5\" \ + 0.1596 * \"$GIS_OPT_BAND7\"" & if [ "$GIS_FLAG_S" -eq 1 ] ; then wait fi g.message "Calculating second TC component $GIS_OPT_OUTPREFIX.2 (Greenness) ..." r.mapcalc "\"$GIS_OPT_OUTPREFIX.2\" = -0.3344 * \"$GIS_OPT_BAND1\" \ - 0.3544 * \"$GIS_OPT_BAND2\" - 0.4556 * \"$GIS_OPT_BAND3\" \ + 0.6966 * \"$GIS_OPT_BAND4\" - 0.0242 * \"$GIS_OPT_BAND5\" \ - 0.2630 * \"$GIS_OPT_BAND7\"" & if [ "$GIS_FLAG_S" -eq 1 ] ; then wait fi g.message "Calculating third TC component $GIS_OPT_OUTPREFIX.3 (Wetness) ..." r.mapcalc "\"$GIS_OPT_OUTPREFIX.3\" = 0.2626 * \"$GIS_OPT_BAND1\" \ + 0.2141 * \"$GIS_OPT_BAND2\" + 0.0926 * \"$GIS_OPT_BAND3\" \ + 0.0656 * \"$GIS_OPT_BAND4\" - 0.7629 * \"$GIS_OPT_BAND5\" \ - 0.5388 * \"$GIS_OPT_BAND7\"" & if [ "$GIS_FLAG_S" -eq 1 ] ; then wait fi g.message "Calculating fourth TC component $GIS_OPT_OUTPREFIX.4. (Haze) ..." r.mapcalc "\"$GIS_OPT_OUTPREFIX.4\" = 0.0805 * \"$GIS_OPT_BAND1\" \ - 0.0498 * \"$GIS_OPT_BAND2\" + 0.1950 * \"$GIS_OPT_BAND3\" \ - 0.1327 * \"$GIS_OPT_BAND4\" + 0.5752 * \"$GIS_OPT_BAND5\" \ - 0.7775 * \"$GIS_OPT_BAND7\"" & wait else g.message -e 'Select LANDSAT satellite by flag!' exit 1 fi r.colors map="$GIS_OPT_OUTPREFIX.1" color=grey r.colors map="$GIS_OPT_OUTPREFIX.2" color=grey r.colors map="$GIS_OPT_OUTPREFIX.3" color=grey if [ $GIS_FLAG_4 -ne 1 ] ; then r.colors map="$GIS_OPT_OUTPREFIX.4" color=grey fi g.message "Tasseled Cap components calculated."