[GRASS-SVN] r63474 - in grass/branches/releasebranch_7_0: . scripts/i.tasscap

svn_grass at osgeo.org svn_grass at osgeo.org
Wed Dec 10 13:17:38 PST 2014


Author: neteler
Date: 2014-12-10 13:17:38 -0800 (Wed, 10 Dec 2014)
New Revision: 63474

Modified:
   grass/branches/releasebranch_7_0/
   grass/branches/releasebranch_7_0/scripts/i.tasscap/i.tasscap.html
   grass/branches/releasebranch_7_0/scripts/i.tasscap/i.tasscap.py
Log:
i.tasscap: MODIS added; partial code rewrite by Leonardo Perathoner, Trento (trunk, r63473)


Property changes on: grass/branches/releasebranch_7_0
___________________________________________________________________
Modified: svn:mergeinfo
   - /grass/trunk:60817,61096,61141,61994,62105,62179-62180,62182,62403,62422,62424,62437,62466,62469,62487,62491,62494,62501,62506,62508-62509,62515,62518-62519,62521,62526,62533,62539,62541,62555,62562,62566,62570,62573,62575,62585,62588,62597,62603,62606,62608-62609,62614,62618,62628,62632,62638,62642,62648-62649,62652,62654-62657,62666,62691,62705,62709,62723,62730,62739,62741,62743,62746,62750,62752,62757,62762,62785,62798,62800-62801,62803,62805,62812,62822,62824,62831,62838,62847,62850,62856,62879,62881,62886,62904,62907-62908,62910,62912,62914,62916,62918,62920,62925,62933,62935,62940,62942,62944-62946,62949,62958,62960,62962,62964,62966-62968,62970,62973,62975,62977,62981,62983,62985,62987,62989,62991,62993,62995,62997,62999-63000,63003,63005,63007,63009,63011,63013,63015,63017,63020,63022,63024,63026,63028-63031,63033,63035,63037,63040,63043-63044,63047,63049,63051,63053,63055,63057,63060,63062-63064,63066,63068,63070-63071,63074,63076,63079,63081,63083,63085,63087,
 63089,63091,63093,63095,63098,63100,63102,63105,63107,63109,63111,63113-63114,63116,63119,63121,63123,63125,63130,63132-63133,63135,63137,63140,63143,63145,63147,63149,63151,63153-63154,63157,63165,63170,63173,63175,63187,63192-63193,63196,63199-63200,63202,63209,63216,63220-63221,63224,63227,63240,63246,63250,63255,63259,63261,63276,63279,63281,63283,63287,63290,63292,63302,63315,63319,63330,63332,63339,63342,63345,63362,63367,63391,63393,63408-63409,63416-63417,63425,63427,63429,63431,63433,63451,63453,63457,63459
   + /grass/trunk:60817,61096,61141,61994,62105,62179-62180,62182,62403,62422,62424,62437,62466,62469,62487,62491,62494,62501,62506,62508-62509,62515,62518-62519,62521,62526,62533,62539,62541,62555,62562,62566,62570,62573,62575,62585,62588,62597,62603,62606,62608-62609,62614,62618,62628,62632,62638,62642,62648-62649,62652,62654-62657,62666,62691,62705,62709,62723,62730,62739,62741,62743,62746,62750,62752,62757,62762,62785,62798,62800-62801,62803,62805,62812,62822,62824,62831,62838,62847,62850,62856,62879,62881,62886,62904,62907-62908,62910,62912,62914,62916,62918,62920,62925,62933,62935,62940,62942,62944-62946,62949,62958,62960,62962,62964,62966-62968,62970,62973,62975,62977,62981,62983,62985,62987,62989,62991,62993,62995,62997,62999-63000,63003,63005,63007,63009,63011,63013,63015,63017,63020,63022,63024,63026,63028-63031,63033,63035,63037,63040,63043-63044,63047,63049,63051,63053,63055,63057,63060,63062-63064,63066,63068,63070-63071,63074,63076,63079,63081,63083,63085,63087,
 63089,63091,63093,63095,63098,63100,63102,63105,63107,63109,63111,63113-63114,63116,63119,63121,63123,63125,63130,63132-63133,63135,63137,63140,63143,63145,63147,63149,63151,63153-63154,63157,63165,63170,63173,63175,63187,63192-63193,63196,63199-63200,63202,63209,63216,63220-63221,63224,63227,63240,63246,63250,63255,63259,63261,63276,63279,63281,63283,63287,63290,63292,63302,63315,63319,63330,63332,63339,63342,63345,63362,63367,63391,63393,63408-63409,63416-63417,63425,63427,63429,63431,63433,63451,63453,63457,63459,63473

Modified: grass/branches/releasebranch_7_0/scripts/i.tasscap/i.tasscap.html
===================================================================
--- grass/branches/releasebranch_7_0/scripts/i.tasscap/i.tasscap.html	2014-12-10 20:09:37 UTC (rev 63473)
+++ grass/branches/releasebranch_7_0/scripts/i.tasscap/i.tasscap.html	2014-12-10 21:17:38 UTC (rev 63474)
@@ -1,15 +1,15 @@
 <h2>DESCRIPTION</h2>
 
 <em>i.tasscap</em> calculates Tasseled Cap (Kauth Thomas, TC) transformation
-for Landsat TM data (TM4, TM5, ETM7).<!-- add other sensors here once added -->
+for Landsat TM data (TM4, TM5, ETM7) and MODIS data.
 
 The tasseled cap transformation is effectively a compression method to
 reduce multiple spectral data into a few bands. The method was originally
 developed for understanding important phenomena of crop development in
 spectral space (Kauth and Thomas 1976).
 <p>
-Tasseled cap coefficients for Landsat 7 ETM+ at-satellite reflectance
-(C. Huang et al., 2001), the conversion can be achieved with
+Tasseled cap coefficients for Landsat 7 ETM+ are at-satellite reflectance
+values (C. Huang et al., 2001), the conversion can be achieved with
 <em>i.landsat.toar</em>.
 <p>
 The following TC components are generated:
@@ -18,7 +18,7 @@
 <li> tasscap.1: corresponds to brightness,
 <li> tasscap.2: corresponds to greenness,
 <li> tasscap.3: corresponds to wetness,
-<li> tasscap.4: corresponds to atmospheric haze.
+<li> tasscap.4: corresponds to atmospheric haze (only selected sensors: Landsat 5,7,8).
 </ul>
 
 
@@ -83,15 +83,15 @@
 <li>Crist, E. P., 1985, A TM tasseled cap equivalent transformation for reflectance
     factor data, Remote Sensing of Environment, 17: 301-306.
 <li>LANDSAT-7: TASSCAP factors cited from:
-  DERIVATION OF A TASSELED CAP TRANSFORMATION BASED ON LANDSAT 7 AT-SATELLITE REFLECTANCE
+  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
 <br>
  This is published as well in INT. J. OF RS, 2002, VOL 23, NO. 8, 1741-1748.
-<br>
- Compare discussion:
- http://adis.cesnet.cz/cgi-bin/lwgate/IMAGRS-L/archives/imagrs-l.log0211/date/article-14.html
+<li> MODIS Tasselled Cap coefficients - Ref: Lobser & Cohen (2007). MODIS tasselled cap:
+ land cover characteristics expressed through transformed MODIS data.
+ International Journal of Remote Sensing, Volume 28(22), Table 3
 </ul>
 
 <h2>SEE ALSO</h2>

Modified: grass/branches/releasebranch_7_0/scripts/i.tasscap/i.tasscap.py
===================================================================
--- grass/branches/releasebranch_7_0/scripts/i.tasscap/i.tasscap.py	2014-12-10 20:09:37 UTC (rev 63473)
+++ grass/branches/releasebranch_7_0/scripts/i.tasscap/i.tasscap.py	2014-12-10 21:17:38 UTC (rev 63474)
@@ -3,19 +3,17 @@
 ############################################################################
 #
 # MODULE:	i.tasscap
-# AUTHOR(S):	Markus Neteler. neteler itc.it
+# AUTHOR(S):	Agustin Lobo, Markus Neteler
 #               Converted to Python by Glynn Clements
+#               Code improvements by Leonardo Perathoner
+#
 # PURPOSE:	At-satellite reflectance based tasseled cap transformation.
-# COPYRIGHT:	(C) 1997-2004,2008, 2014 by the GRASS Development Team
+# COPYRIGHT:	(C) 1997-2014 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
-#       Add other, e.g from here: http://www.sjsu.edu/faculty/watkins/tassel.htm
-#
 #############################################################################
 # References:
 # LANDSAT-4/LANDSAT-5:
@@ -34,10 +32,16 @@
 #  Compare discussion:
 #  http://adis.cesnet.cz/cgi-bin/lwgate/IMAGRS-L/archives/imagrs-l.log0211/date/article-14.html
 #
-#   Landsat8: Baig, M.H.A., Zhang, L., Shuai, T., Tong, Q., 2014. Derivation of a tasselled cap transformation
+#  Landsat8: Baig, M.H.A., Zhang, L., Shuai, T., Tong, Q., 2014. Derivation of a tasselled cap transformation
 #              based on Landsat 8 at-satellite reflectance. Remote Sensing Letters 5, 423-431. 
 #              doi:10.1080/2150704X.2014.915434
 #
+#  MODIS Tasselled Cap coefficients
+#  https://gis.stackexchange.com/questions/116107/tasseled-cap-transformation-on-modis-in-grass/116110
+#  Ref: Lobser & Cohen (2007). MODIS tasselled cap: land cover characteristics 
+#                              expressed through transformed MODIS data.
+#       International Journal of Remote Sensing, Volume 28(22), Table 3
+#
 #############################################################################
 #
 #%Module
@@ -45,10 +49,11 @@
 #% keywords: imagery
 #% keywords: transformation
 #% keywords: Landsat
+#% keywords: MODIS
 #% keywords: Tasseled Cap transformation
 #%end
 #%option G_OPT_R_INPUTS
-#% description: For Landsat4-7: bands 1, 2, 3, 4, 5, and 7; for Landsat8: bands 2, 3, 4, 5, 6, and 7
+#% description: For Landsat4-7: bands 1, 2, 3, 4, 5, 7; for Landsat8: bands 2, 3, 4, 5, 6, 7; for MODIS: bands 1, 2, 3, 4, 5, 6, 7
 #%end
 #%option G_OPT_R_BASENAME_OUTPUT
 #% label: Name for output basename raster map(s)
@@ -59,13 +64,14 @@
 #% description: Satellite sensor
 #% required: yes
 #% multiple: no
-#% options: landsat4_tm,landsat5_tm,landsat7_etm,landsat8
-#% descriptions: landsat4_tm;Use transformation rules for Landsat 4 TM;landsat5_tm;Use transformation rules for Landsat 5 TM;landsat7_etm;Use transformation rules for Landsat 7 ETM;landsat8;Use transformation rules for Landsat 8
+#% options: landsat4_tm,landsat5_tm,landsat7_etm,landsat8_oli,modis
+#% descriptions: landsat4_tm;Use transformation rules for Landsat 4 TM;landsat5_tm;Use transformation rules for Landsat 5 TM;landsat7_etm;Use transformation rules for Landsat 7 ETM;landsat8_oli;Use transformation rules for Landsat 8 OLI;modis;Use transformation rules for MODIS
 #%end
 
 import grass.script as grass
 
 # weights for 6 Landsat bands: TM4, TM5, TM7, OLI
+# MODIS: Red, NIR1, Blue, Green, NIR2, SWIR1, SWIR2
 parms = [[( 0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863), # Landsat TM4
 	  (-0.2848,-0.2435,-0.5435, 0.7243, 0.0840,-0.1800),
 	  ( 0.1509, 0.1973, 0.3279, 0.3406,-0.7112,-0.4572)],
@@ -80,27 +86,36 @@
 	 [( 0.3029, 0.2786, 0.4733, 0.5599, 0.5080, 0.1872), # Landsat TM8
 	  (-0.2941,-0.2430,-0.5424, 0.7276, 0.0713,-0.1608),
 	  ( 0.1511, 0.1973, 0.3283, 0.3407,-0.7117,-0.4559),
-	  (-0.8239, 0.0849, 0.4396, -0.058, 0.2013,-0.2773)]]
+	  (-0.8239, 0.0849, 0.4396, -0.058, 0.2013,-0.2773)],
+	 [( 0.4395, 0.5945, 0.2460, 0.3918, 0.3506, 0.2136, 0.2678), # MODIS
+	  (-0.4064, 0.5129,-0.2744,-0.2893, 0.4882,-0.0036,-0.4169),
+	  ( 0.1147, 0.2489, 0.2408, 0.3132,-0.3122,-0.6416,-0.5087)]]
+#satellite information
+satellites = ["landsat4_tm", 'landsat5_tm', 'landsat7_etm', 'landsat8_oli', 'modis']
+used_bands = [6,6,6,6,7]
+#components information
 ordinals = ["first", "second", "third", "fourth"]
 names = ["Brightness", "Greenness", "Wetness", "Haze"]
 
-def calc1(out, bands, k1, k2, k3, k4, k5, k7, k0 = 0):
-    grass.mapcalc(
-	"$out = $k1 * $band1 + $k2 * $band2 + $k3 * $band3 + $k4 * $band4 + $k5 * $band5 + $k7 * $band7 + $k0",
-	out = out, k1 = k1, k2 = k2, k3 = k3, k4 = k4, k5 = k5, k7 = k7, k0 = k0, **bands)
-    grass.run_command('r.colors', map = out, color = 'grey')
+def calc1bands6(out, bands, k1, k2, k3, k4, k5, k6, k0 = 0):
+    grass.mapcalc("$out = $k1 * $in1band + $k2 * $in2band + $k3 * $in3band + $k4 * $in4band + $k5 * $in5band + $k6 * $in6band + $k0",
+    out = out, k1 = k1, k2 = k2, k3 = k3, k4 = k4, k5 = k5, k6 = k6, k0 = k0, **bands)    
 
-def calcN(outpre, bands, i, n):
-    grass.message(_("Landsat-%d...") % n)
+def calc1bands7(out, bands, k1, k2, k3, k4, k5, k6, k7):
+    grass.mapcalc("$out = $k1 * $in1band + $k2 * $in2band + $k3 * $in3band + $k4 * $in4band + $k5 * $in5band + $k6 * $in6band + $k7 * $in7band",
+    out = out, k1 = k1, k2 = k2, k3 = k3, k4 = k4, k5 = k5, k6 = k6, k7 = k7, **bands)
+
+def calcN(outpre, bands, satel):
+    i=satellites.index(satel)
+    grass.message(_("Satellite %s...") % satel)        
     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, bands, *p)
+        out = "%s.%d" % (outpre, j + 1)
+        ord = ordinals[j]
+        name = " (%s)" % names[j]        
+        grass.message(_("Calculating %s TC component %s%s ...") % (ord, out, name))
+        bands_num=used_bands[i]
+        eval("calc1bands%d(out, bands, *p)" % bands_num) #use combination function suitable for used number of bands
+        grass.run_command('r.colors', map = out, color = 'grey', quiet=True)
 
 
 def main():
@@ -108,35 +123,21 @@
     satellite = options['sensor']
     output_basename = options['output']
     inputs = options['input'].split(',')
-    num_of_bands = 6
+    num_of_bands = used_bands[satellites.index(satellite)]
     if len(inputs) != num_of_bands:
         grass.fatal(_("The number of input raster maps (bands) should be %s") % num_of_bands)
 
-    # this is here just for the compatibility with r.mapcalc expression
-    # remove this if not really needed in new implementation
     bands = {}
-    for i, band in enumerate(inputs):
-        band_num = i + 1
-        if band_num == 6:
-            band_num = 7
-        bands['band' + str(band_num)] = band
+    for i, band in enumerate(inputs):        
+        band_num = i + 1        
+        bands['in' + str(band_num) + 'band'] = band
     grass.debug(1, bands)
 
-    if satellite == 'landsat4_tm':
-        calcN(output_basename, bands, 0, 4)
-    elif satellite == 'landsat5_tm':
-        calcN(output_basename, bands, 1, 5)
-    elif satellite == 'landsat7_etm':
-        calcN(output_basename, bands, 2, 7)
-    elif satellite == 'landsat8':
-        calcN(output_basename, bands, 3, 8)
-    else:
-        raise RuntimeError("Invalid satellite: " + satellite)
+    calcN(output_basename, bands, satellite) #core tasseled cap components computation
 
-    grass.run_command('r.support', map = "%s.%d" % (output_basename, 1), description = "Tasseled Cap 1: brightness")
-    grass.run_command('r.support', map = "%s.%d" % (output_basename, 2), description = "Tasseled Cap 2: greenness")
-    grass.run_command('r.support', map = "%s.%d" % (output_basename, 3), description = "Tasseled Cap 3: wetness")
-    grass.run_command('r.support', map = "%s.%d" % (output_basename, 4), description = "Tasseled Cap 4: atmospheric haze")
+    #assign "Data Description" field in all four component maps
+    for i, comp in enumerate(names):
+        grass.run_command('r.support', map = "%s.%d" % (output_basename, i+1), description = "Tasseled Cap %d: %s" % (i+1, comp))
 
     grass.message(_("Tasseled Cap components calculated"))
 



More information about the grass-commit mailing list