-
Notifications
You must be signed in to change notification settings - Fork 10
/
tct.R
137 lines (131 loc) · 4.79 KB
/
tct.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
#' Tasseled Cap Transformation (TCT)
#'
#' The Tasseled-Cap Transformation is a linear transformation method for various
#' remote sensing data. Not only can it perform volume data compression, but it
#' can also provide parameters associated with the physical characteristics,
#' such as brightness, greenness and wetness indices.
#'
#' @references
#' Crist, E.P., R. Laurin, and R.C. Cicone. 1986. Vegetation and soils information
#' contained in transformed Thematic Mapper data. Pages 1465-1470 Ref. ESA SP-254.
#' European Space Agency, Paris, France. http://www.ciesin.org/docs/005-419/005-419.html.
#'
#' Baig, M.H.A., Shuai, T., Tong, Q., 2014. Derivation of a tasseled cap transformation
#' based on Landsat 8 at-satellite reflectance. Remote Sensing Letters, 5(5), 423-431.
#'
#' Li, B., Ti, C., Zhao, Y., Yan, X., 2016. Estimating Soil Moisture with Landsat Data
#' and Its Application in Extracting the Spatial Distribution of Winter Flooded Paddies.
#' Remote Sensing, 8(1), 38.
#'
#' @note
#' Currently implemented for satellites such as Landsat-4 TM, Landsat-5 TM, Landsat-7 ETM+,
#' Landsat-8 OLI and Sentinel2. The input data must be in top of atmosphere reflectance (toa).
#' Bands required as input must be ordered as:
#'
#' \tabular{ll}{
#' Type of Sensor \tab Name of bands \cr
#' Landsat4TM \tab blue,green,red,nir,swir1,swir2 \cr
#' Landsat5TM \tab blue,green,red,nir,swir1,swir2 \cr
#' Landsat7ETM+ \tab blue,green,red,nir,swir1,swir2 \cr
#' Landsat8OLI \tab blue,green,red,nir,swir1,swir2 \cr
#' Landsat8OLI-Li2016 \tab coastal,blue,green,red,nir,swir1,swir2 \cr
#' Sentinel2MSI \tab coastal,blue,green,red,nir-1,mir-1,mir-2 \cr
#' }
#'
#' @param img It could be RasterStack or RasterBrick.
#' @param sat Specify satellite and sensor type (Landsat5TM, Landsat7ETM
#' or Landsat8OLI).
#'
#' @importFrom raster as.matrix
#'
#' @examples
#' library(ForesToolboxRS)
#' library(raster)
#'
#' # Load an example dataset
#' data(img_l8)
#'
#' # Tasseled-cap using Landsat8OLI
#' sat_tct <- tct(img = img_l8, sat = "Landsat8OLI")
#' plotRGB(sat_tct, 1, 2, 3, stretch = "lin")
#' @export
#'
tct <- function(img, sat = "Landsat8OLI") {
if (!inherits(img, "Raster")) stop("img must be a RasterBrick or RasterStack", call. = TRUE)
if (sat == "Landsat4TM") {
coefc <- matrix(c(
0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863,
-0.2848, -0.2435, -0.5436, 0.7243, 0.0840, -0.1800,
0.1509, 0.1973, 0.3279, 0.3406, -0.7112, -0.4572
), 3, 6,
byrow = TRUE, dimnames = list(
c("Brightness", "Greenness", "Wetness"),
c("B1", "B2", "B3", "B4", "B5", "B7")
)
)
} else if (sat == "Landsat5TM") {
coefc <- matrix(c(
0.2909, 0.2493, 0.4806, 0.5568, 0.4438, 0.1706,
-0.2728, -0.2174, -0.5508, 0.7221, 0.0733, -0.1648,
0.1446, 0.1761, 0.3322, 0.3396, -0.6210, -0.4186
), 3, 6,
byrow = TRUE, dimnames = list(
c("Brightness", "Greenness", "Wetness"),
c("B1", "B2", "B3", "B4", "B5", "B7")
)
)
} else if (sat == "Landsat7ETM") {
coefc <- matrix(c(
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
), 3, 6,
byrow = TRUE, dimnames = list(
c("Brightness", "Greenness", "Wetness"),
c("B1", "B2", "B3", "B4", "B5", "B7")
)
)
} else if (sat == "Landsat8OLI") {
coefc <- matrix(c(
0.3029, 0.2786, 0.4733, 0.5599, 0.5080, 0.1872,
-0.2941, -0.2430, -0.5424, 0.7276, 0.0713, -0.1608,
0.1511, 0.1973, 0.3283, 0.3407, -0.7117, -0.4559
), 3, 6,
byrow = TRUE, dimnames = list(
c("Brightness", "Greenness", "Wetness"),
c("B1", "B2", "B3", "B4", "B5", "B7")
)
)
} else if (sat == "Landsat8OLI-Li2016") {
coefc <- matrix(c(
0.2540, 0.3037, 0.3608, 0.3564, 0.7084, 0.2358, 0.1691,
-0.2578, -0.3064, -0.3300, -0.4325, 0.6860, -0.0383, -0.2674,
0.1877, 0.2097, 0.2038, 0.1017, 0.0685, -0.7460, -0.5548
), 3, 7,
byrow = TRUE, dimnames = list(
c("Brightness", "Greenness", "Wetness"),
c("B1", "B2", "B3", "B4", "B5", "B6", "B7")
)
)
} else if (sat == "Sentinel2MSI") {
coefc <- matrix(c(
0.2381, 0.2569, 0.2934, 0.3020, 0.3580, 0.0896, 0.0780,
-0.2266, -0.2818, -0.3020, -0.4283, 0.3138, -0.1341, -0.2538,
0.1825, 0.1763, 0.1615, 0.0486, -0.0755, -0.7701, -0.5293
), 3, 7,
byrow = TRUE, dimnames = list(
c("Brightness", "Greenness", "Wetness"),
c("B1", "B2", "B3", "B4", "B8", "B11", "B12")
)
)
} else {
stop("Satellite not supported.")
}
if(dim(img)[3] != dim(coefc)[2]) stop(" The number of bands must be equal to the number of coefficients in band.", call. = TRUE)
val <- as.matrix(img) %*% t(coefc)
bgw <- img[[1:3]]
for (i in 1:3) {
values(bgw[[i]]) <- val[, i]
}
return(bgw)
}