[GRASS-SVN] r33623 - grass/trunk/scripts/i.tasscap
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Sep 30 14:09:19 EDT 2008
Author: glynn
Date: 2008-09-30 14:09:19 -0400 (Tue, 30 Sep 2008)
New Revision: 33623
Added:
grass/trunk/scripts/i.tasscap/i.tasscap.py
Modified:
grass/trunk/scripts/i.tasscap/i.tasscap
Log:
Convert i.tasscap to Python (and convert shell version to work with newer r.mapcalc)
Modified: grass/trunk/scripts/i.tasscap/i.tasscap
===================================================================
--- grass/trunk/scripts/i.tasscap/i.tasscap 2008-09-30 17:27:36 UTC (rev 33622)
+++ grass/trunk/scripts/i.tasscap/i.tasscap 2008-09-30 18:09:19 UTC (rev 33623)
@@ -120,30 +120,30 @@
if [ $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"
+ 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"
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"
+ 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"
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"
+ 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"
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"
+ 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"
else
if [ $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"
+ 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"
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"
+ 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"
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"
+ 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"
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"
+ 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"
else
g.message -e 'Select LANDSAT satellite by flag!'
exit 1
Added: grass/trunk/scripts/i.tasscap/i.tasscap.py
===================================================================
--- grass/trunk/scripts/i.tasscap/i.tasscap.py (rev 0)
+++ grass/trunk/scripts/i.tasscap/i.tasscap.py 2008-09-30 18:09:19 UTC (rev 33623)
@@ -0,0 +1,160 @@
+#!/usr/bin/env python
+#
+############################################################################
+#
+# MODULE: i.tasscap
+# AUTHOR(S): Markus Neteler. neteler itc.it
+# Converted to Python by Glynn Clements
+# PURPOSE: At-satellite reflectance based tasseled cap transformation.
+# COPYRIGHT: (C) 1997-2004,2008 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 at 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
+
+import sys
+import os
+import string
+import grass
+
+parms = [[( 0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863),
+ (-0.2848,-0.2435,-0.5435, 0.7243, 0.0840,-0.1800),
+ ( 0.1509, 0.1973, 0.3279, 0.3406,-0.7112,-0.4572)],
+ [( 0.2909, 0.2493, 0.4806, 0.5568, 0.4438, 0.1706, 10.3695),
+ (-0.2728,-0.2174,-0.5508, 0.7221, 0.0733,-0.1648, -0.7310),
+ ( 0.1446, 0.1761, 0.3322, 0.3396,-0.6210,-0.4186, -3.3828),
+ ( 0.8461,-0.0731,-0.4640,-0.0032,-0.0492,-0.0119, 0.7879)],
+ [( 0.3561, 0.3972, 0.3904, 0.6966, 0.2286, 0.1596),
+ (-0.3344,-0.3544,-0.4556, 0.6966,-0.0242,-0.2630),
+ ( 0.2626, 0.2141, 0.0926, 0.0656,-0.7629,-0.5388),
+ ( 0.0805,-0.0498, 0.1950,-0.1327, 0.5752,-0.7775)]]
+ordinals = ["first", "second", "third", "fourth"]
+names = ["Brightness", "Greenness", "Wetness", "Haze"]
+
+def calc1(out, bands, k1, k2, k3, k4, k5, k7, k0 = 0):
+ t = string.Template("$out = $k1 * $band1 + $k2 * $band2 + $k3 * $band3 + $k4 * $band4 + $k5 * $band5 + $k7 * $band7 + $k0")
+ e = t.substitute(out = out, k1 = k1, k2 = k2, k3 = k3, k4 = k4, k5 = k5, k7 = k7, k0 = k0, **bands)
+ grass.run_command('r.mapcalc', expression = e)
+ grass.run_command('r.colors', map = out, color = 'grey')
+
+def calcN(options, i, n):
+ outpre = options['outprefix']
+ grass.message("LANDSAT-%d..." % n)
+ for j, p in enumerate(parms[i]):
+ out = "%s.%d" % (outpre, j + 1)
+ ord = ordinals[j]
+ if n == 4:
+ name = ''
+ else:
+ name = " (%s)" % names[j]
+ grass.message("Calculating %s TC component %s%s ..." % (ord, out, name))
+ calc1(out, options, *p)
+
+def main():
+ flag4 = flags['4']
+ flag5 = flags['5']
+ flag7 = flags['7']
+
+ if (flag4 and flag5) or (flag4 and flag7) or (flag5 and flag7):
+ grass.fatal("Select only one flag")
+
+ if flag4:
+ calcN(options, 0, 4)
+ elif flag5:
+ calcN(options, 1, 5)
+ elif flag7:
+ calcN(options, 2, 7)
+ else:
+ grass.fatal("Select LANDSAT satellite by flag!")
+
+ grass.message("Tasseled Cap components calculated.")
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main()
Property changes on: grass/trunk/scripts/i.tasscap/i.tasscap.py
___________________________________________________________________
Name: svn:executable
+ *
More information about the grass-commit
mailing list