i.tasscap.py 7.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182
  1. #!/usr/bin/env python
  2. ############################################################################
  3. #
  4. # MODULE: i.tasscap
  5. # AUTHOR(S): Agustin Lobo, Markus Neteler
  6. # Converted to Python by Glynn Clements
  7. # Code improvements by Leonardo Perathoner
  8. #
  9. # PURPOSE: At-satellite reflectance based tasseled cap transformation.
  10. # COPYRIGHT: (C) 1997-2014 by the GRASS Development Team
  11. #
  12. # This program is free software under the GNU General Public
  13. # License (>=v2). Read the file COPYING that comes with GRASS
  14. # for details.
  15. #
  16. #############################################################################
  17. # References:
  18. # LANDSAT-4/LANDSAT-5:
  19. # script based on i.tasscap.tm4 from Dr. Agustin Lobo - alobo@ija.csic.es
  20. # TC-factor changed to CRIST et al. 1986, p.1467 (Markus Neteler 1/99)
  21. # Proc. IGARSS 1986
  22. #
  23. # LANDSAT-7:
  24. # TASSCAP factors cited from:
  25. # DERIVATION OF A TASSELED CAP TRANSFORMATION BASED ON LANDSAT 7 AT-SATELLITE REFLECTANCE
  26. # Chengquan Huang, Bruce Wylie, Limin Yang, Collin Homer and Gregory Zylstra Raytheon ITSS,
  27. # USGS EROS Data Center Sioux Falls, SD 57198, USA
  28. # http://landcover.usgs.gov/pdf/tasseled.pdf
  29. #
  30. # This is published as well in INT. J. OF RS, 2002, VOL 23, NO. 8, 1741-1748.
  31. # Compare discussion:
  32. # http://adis.cesnet.cz/cgi-bin/lwgate/IMAGRS-L/archives/imagrs-l.log0211/date/article-14.html
  33. #
  34. # Landsat8: Baig, M.H.A., Zhang, L., Shuai, T., Tong, Q., 2014. Derivation of a tasselled cap transformation
  35. # based on Landsat 8 at-satellite reflectance. Remote Sensing Letters 5, 423-431.
  36. # doi:10.1080/2150704X.2014.915434
  37. #
  38. # MODIS Tasselled Cap coefficients
  39. # https://gis.stackexchange.com/questions/116107/tasseled-cap-transformation-on-modis-in-grass/116110
  40. # Ref: Lobser & Cohen (2007). MODIS tasselled cap: land cover characteristics
  41. # expressed through transformed MODIS data.
  42. # International Journal of Remote Sensing, Volume 28(22), Table 3
  43. #
  44. #############################################################################
  45. #%Module
  46. #% description: Performs Tasseled Cap (Kauth Thomas) transformation.
  47. #% keyword: imagery
  48. #% keyword: transformation
  49. #% keyword: Landsat
  50. #% keyword: MODIS
  51. #% keyword: Tasseled Cap transformation
  52. #%end
  53. #%option G_OPT_R_INPUTS
  54. #% 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
  55. #%end
  56. #%option G_OPT_R_BASENAME_OUTPUT
  57. #% label: Name for output basename raster map(s)
  58. #%end
  59. #%option
  60. #% key: sensor
  61. #% type: string
  62. #% description: Satellite sensor
  63. #% required: yes
  64. #% multiple: no
  65. #% options: landsat4_tm,landsat5_tm,landsat7_etm,landsat8_oli,modis
  66. #% 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
  67. #%end
  68. import grass.script as grass
  69. # i18N
  70. import os
  71. import gettext
  72. gettext.install('grassmods', os.path.join(os.getenv("GISBASE"), 'locale'))
  73. # weights for 6 Landsat bands: TM4, TM5, TM7, OLI
  74. # MODIS: Red, NIR1, Blue, Green, NIR2, SWIR1, SWIR2
  75. parms = [[(0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863), # Landsat TM4
  76. (-0.2848, -0.2435, -0.5435, 0.7243, 0.0840, -0.1800),
  77. (0.1509, 0.1973, 0.3279, 0.3406, -0.7112, -0.4572)],
  78. [(0.2909, 0.2493, 0.4806, 0.5568, 0.4438, 0.1706, 10.3695), # Landsat TM5
  79. (-0.2728, -0.2174, -0.5508, 0.7221, 0.0733, -0.1648, -0.7310),
  80. (0.1446, 0.1761, 0.3322, 0.3396, -0.6210, -0.4186, -3.3828),
  81. (0.8461, -0.0731, -0.4640, -0.0032, -0.0492, -0.0119, 0.7879)],
  82. [(0.3561, 0.3972, 0.3904, 0.6966, 0.2286, 0.1596), # Landsat TM7
  83. (-0.3344, -0.3544, -0.4556, 0.6966, -0.0242, -0.2630),
  84. (0.2626, 0.2141, 0.0926, 0.0656, -0.7629, -0.5388),
  85. (0.0805, -0.0498, 0.1950, -0.1327, 0.5752, -0.7775)],
  86. [(0.3029, 0.2786, 0.4733, 0.5599, 0.5080, 0.1872), # Landsat OLI
  87. (-0.2941, -0.2430, -0.5424, 0.7276, 0.0713, -0.1608),
  88. (0.1511, 0.1973, 0.3283, 0.3407, -0.7117, -0.4559),
  89. (-0.8239, 0.0849, 0.4396, -0.0580, 0.2013, -0.2773)],
  90. [(0.4395, 0.5945, 0.2460, 0.3918, 0.3506, 0.2136, 0.2678), # MODIS
  91. (-0.4064, 0.5129, -0.2744, -0.2893, 0.4882, -0.0036, -0.4169),
  92. (0.1147, 0.2489, 0.2408, 0.3132, -0.3122, -0.6416, -0.5087)]]
  93. # satellite information
  94. satellites = ['landsat4_tm', 'landsat5_tm', 'landsat7_etm', 'landsat8_oli',
  95. 'modis']
  96. used_bands = [6, 6, 6, 6, 7]
  97. # components information
  98. ordinals = ["first", "second", "third", "fourth"]
  99. names = ["Brightness", "Greenness", "Wetness", "Haze"]
  100. def calc1bands6(out, bands, k1, k2, k3, k4, k5, k6, k0=0):
  101. """
  102. Tasseled cap transformation equation for Landsat bands
  103. """
  104. equation = ('$out = $k1 * $in1band + $k2 * $in2band + $k3 * $in3band + '
  105. '$k4 * $in4band + $k5 * $in5band + $k6 * $in6band + $k0')
  106. grass.mapcalc(equation, out=out, k1=k1, k2=k2, k3=k3, k4=k4, k5=k5,
  107. k6=k6, k0=k0, **bands)
  108. def calc1bands7(out, bands, k1, k2, k3, k4, k5, k6, k7):
  109. """
  110. Tasseled cap transformation equation for MODIS bands
  111. """
  112. equation = ('$out = $k1 * $in1band + $k2 * $in2band + $k3 * $in3band + '
  113. '$k4 * $in4band + $k5 * $in5band + $k6 * $in6band + $k7 * '
  114. '$in7band')
  115. grass.mapcalc(equation, out=out, k1=k1, k2=k2, k3=k3, k4=k4, k5=k5, k6=k6,
  116. k7=k7, **bands)
  117. def calcN(outpre, bands, satel):
  118. """
  119. Calculating Tasseled Cap components
  120. """
  121. i = satellites.index(satel)
  122. grass.message(_("Satellite %s...") % satel)
  123. for j, p in enumerate(parms[i]):
  124. out = "%s.%d" % (outpre, j + 1)
  125. ord = ordinals[j]
  126. name = " (%s)" % names[j]
  127. message = "Calculating {ordinal} TC component {outprefix}{outname} ..."
  128. message = message.format(ordinal=ord, outprefix=out, outname=name)
  129. grass.message(_(message))
  130. bands_num = used_bands[i]
  131. # use combination function suitable for used number of bands
  132. eval("calc1bands%d(out, bands, *p)" % bands_num)
  133. grass.run_command('r.colors', map=out, color='grey', quiet=True)
  134. def main():
  135. options, flags = grass.parser()
  136. satellite = options['sensor']
  137. output_basename = options['output']
  138. inputs = options['input'].split(',')
  139. num_of_bands = used_bands[satellites.index(satellite)]
  140. if len(inputs) != num_of_bands:
  141. grass.fatal(_("The number of input raster maps (bands) should be %s") % num_of_bands)
  142. bands = {}
  143. for i, band in enumerate(inputs):
  144. band_num = i + 1
  145. bands['in' + str(band_num) + 'band'] = band
  146. grass.debug(bands, 1)
  147. # core tasseled cap components computation
  148. calcN(output_basename, bands, satellite)
  149. # assign "Data Description" field in all four component maps
  150. for i, comp in enumerate(names):
  151. grass.run_command('r.support', map="%s.%d" % (output_basename, i + 1),
  152. description="Tasseled Cap %d: %s" % (i + 1, comp))
  153. grass.message(_("Tasseled Cap components calculated"))
  154. if __name__ == "__main__":
  155. main()