← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to scripts/i.tasscap/i.tasscap.py

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/i.tasscap/i.tasscap.py
 
1
#!/usr/bin/env python
 
2
#
 
3
############################################################################
 
4
#
 
5
# MODULE:       i.tasscap
 
6
# AUTHOR(S):    Agustin Lobo, Markus Neteler
 
7
#               Converted to Python by Glynn Clements
 
8
#               Code improvements by Leonardo Perathoner
 
9
#
 
10
# PURPOSE:      At-satellite reflectance based tasseled cap transformation.
 
11
# COPYRIGHT:    (C) 1997-2014 by the GRASS Development Team
 
12
#
 
13
#               This program is free software under the GNU General Public
 
14
#               License (>=v2). Read the file COPYING that comes with GRASS
 
15
#               for details.
 
16
#
 
17
#############################################################################
 
18
# References:
 
19
# LANDSAT-4/LANDSAT-5:
 
20
#  script based on i.tasscap.tm4 from Dr. Agustin Lobo - alobo@ija.csic.es
 
21
#  TC-factor changed to CRIST et al. 1986, p.1467 (Markus Neteler 1/99)
 
22
#                       Proc. IGARSS 1986
 
23
#
 
24
# LANDSAT-7:
 
25
#  TASSCAP factors cited from:
 
26
#   DERIVATION OF A TASSELED CAP TRANSFORMATION BASED ON LANDSAT 7 AT-SATELLITE REFLECTANCE
 
27
#   Chengquan Huang, Bruce Wylie, Limin Yang, Collin Homer and Gregory Zylstra Raytheon ITSS, 
 
28
#   USGS EROS Data Center Sioux Falls, SD 57198, USA
 
29
#   http://landcover.usgs.gov/pdf/tasseled.pdf
 
30
#
 
31
#  This is published as well in INT. J. OF RS, 2002, VOL 23, NO. 8, 1741-1748.
 
32
#  Compare discussion:
 
33
#  http://adis.cesnet.cz/cgi-bin/lwgate/IMAGRS-L/archives/imagrs-l.log0211/date/article-14.html
 
34
#
 
35
#  Landsat8: Baig, M.H.A., Zhang, L., Shuai, T., Tong, Q., 2014. Derivation of a tasselled cap transformation
 
36
#              based on Landsat 8 at-satellite reflectance. Remote Sensing Letters 5, 423-431. 
 
37
#              doi:10.1080/2150704X.2014.915434
 
38
#
 
39
#  MODIS Tasselled Cap coefficients
 
40
#  https://gis.stackexchange.com/questions/116107/tasseled-cap-transformation-on-modis-in-grass/116110
 
41
#  Ref: Lobser & Cohen (2007). MODIS tasselled cap: land cover characteristics 
 
42
#                              expressed through transformed MODIS data.
 
43
#       International Journal of Remote Sensing, Volume 28(22), Table 3
 
44
#
 
45
#############################################################################
 
46
#
 
47
#%Module
 
48
#% description: Performs Tasseled Cap (Kauth Thomas) transformation.
 
49
#% keyword: imagery
 
50
#% keyword: transformation
 
51
#% keyword: Landsat
 
52
#% keyword: MODIS
 
53
#% keyword: Tasseled Cap transformation
 
54
#%end
 
55
#%option G_OPT_R_INPUTS
 
56
#% 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
 
57
#%end
 
58
#%option G_OPT_R_BASENAME_OUTPUT
 
59
#% label: Name for output basename raster map(s)
 
60
#%end
 
61
#%option
 
62
#% key: sensor
 
63
#% type: string
 
64
#% description: Satellite sensor
 
65
#% required: yes
 
66
#% multiple: no
 
67
#% options: landsat4_tm,landsat5_tm,landsat7_etm,landsat8_oli,modis
 
68
#% 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
 
69
#%end
 
70
 
 
71
import grass.script as grass
 
72
 
 
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 TM8
 
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.058, 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', 'modis']
 
95
used_bands = [6,6,6,6,7]
 
96
#components information
 
97
ordinals = ["first", "second", "third", "fourth"]
 
98
names = ["Brightness", "Greenness", "Wetness", "Haze"]
 
99
 
 
100
def calc1bands6(out, bands, k1, k2, k3, k4, k5, k6, k0 = 0):
 
101
    grass.mapcalc("$out = $k1 * $in1band + $k2 * $in2band + $k3 * $in3band + $k4 * $in4band + $k5 * $in5band + $k6 * $in6band + $k0",
 
102
    out = out, k1 = k1, k2 = k2, k3 = k3, k4 = k4, k5 = k5, k6 = k6, k0 = k0, **bands)    
 
103
 
 
104
def calc1bands7(out, bands, k1, k2, k3, k4, k5, k6, k7):
 
105
    grass.mapcalc("$out = $k1 * $in1band + $k2 * $in2band + $k3 * $in3band + $k4 * $in4band + $k5 * $in5band + $k6 * $in6band + $k7 * $in7band",
 
106
    out = out, k1 = k1, k2 = k2, k3 = k3, k4 = k4, k5 = k5, k6 = k6, k7 = k7, **bands)
 
107
 
 
108
def calcN(outpre, bands, satel):
 
109
    i=satellites.index(satel)
 
110
    grass.message(_("Satellite %s...") % satel)        
 
111
    for j, p in enumerate(parms[i]):
 
112
        out = "%s.%d" % (outpre, j + 1)
 
113
        ord = ordinals[j]
 
114
        name = " (%s)" % names[j]        
 
115
        grass.message(_("Calculating %s TC component %s%s ...") % (ord, out, name))
 
116
        bands_num=used_bands[i]
 
117
        eval("calc1bands%d(out, bands, *p)" % bands_num) #use combination function suitable for used number of bands
 
118
        grass.run_command('r.colors', map = out, color = 'grey', quiet=True)
 
119
 
 
120
 
 
121
def main():
 
122
    options, flags = grass.parser()
 
123
    satellite = options['sensor']
 
124
    output_basename = options['output']
 
125
    inputs = options['input'].split(',')
 
126
    num_of_bands = used_bands[satellites.index(satellite)]
 
127
    if len(inputs) != num_of_bands:
 
128
        grass.fatal(_("The number of input raster maps (bands) should be %s") % num_of_bands)
 
129
 
 
130
    bands = {}
 
131
    for i, band in enumerate(inputs):        
 
132
        band_num = i + 1        
 
133
        bands['in' + str(band_num) + 'band'] = band
 
134
    grass.debug(1, bands)
 
135
 
 
136
    calcN(output_basename, bands, satellite) #core tasseled cap components computation
 
137
 
 
138
    #assign "Data Description" field in all four component maps
 
139
    for i, comp in enumerate(names):
 
140
        grass.run_command('r.support', map = "%s.%d" % (output_basename, i+1), description = "Tasseled Cap %d: %s" % (i+1, comp))
 
141
 
 
142
    grass.message(_("Tasseled Cap components calculated"))
 
143
 
 
144
if __name__ == "__main__":
 
145
    main()