forked from bleutner/RStoolbox
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathspectralIndices.cpp
More file actions
executable file
·157 lines (142 loc) · 5.39 KB
/
spectralIndices.cpp
File metadata and controls
executable file
·157 lines (142 loc) · 5.39 KB
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
#include <Rcpp.h>
using namespace Rcpp;
//[[Rcpp::export]]
NumericMatrix spectralIndicesCpp(NumericMatrix& x, CharacterVector& indices,
const int redBand, const int blueBand, const int greenBand, const int nirBand, const int swir2Band, const int swir1Band,
const double L, const double s, const double G, const double C1, const double C2, const double Levi) {
int nind = indices.size();
int nsamp = x.nrow();
/***
for(int ro = 0; ro < nsamp; ++ro) {
for(int co = 0; co < nind; ++co) {
if(ISNAN(x(ro,co))) x(ro,co) = NA_REAL;
}
}
***/
NumericMatrix out(nsamp, nind);
NumericVector blue, green, red, nir, swir1, swir2;
if(blueBand != NA_INTEGER) blue = x(_,blueBand - 1);
if(greenBand != NA_INTEGER) green = x(_,greenBand - 1);
if(redBand != NA_INTEGER) red = x(_,redBand - 1);
if(nirBand != NA_INTEGER) nir = x(_,nirBand - 1);
if(swir1Band != NA_INTEGER) swir1 = x(_,swir1Band - 1);
if(swir2Band != NA_INTEGER) swir2 = x(_,swir2Band - 1);
for(int j = 0; j < nind; ++j) {
if(indices[j] == "DVI") {
// Difference vegeation index
out(_,j) = (s * nir - red);
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "CTVI") {
// Corrected transformed vegetation index
// Perry and Lautenschlager 1984
NumericVector np = (nir-red)/(nir+red) + 0.5;
out(_,j) = np / sqrt(abs(np));
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "EVI") {
// Enhanced vegetation index
// Huete et al 1990
out(_,j) = G * ((nir - red) / (nir + C1 * red - C2 * blue + Levi));
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "GEMI") {
out(_,j) = (((pow(nir, 2) - pow(red, 2)) * 2.0 + (nir * 1.5) + (red * 0.5) ) /
(nir + red + 0.5)) * (1 - ((((pow(nir,2) - pow(red,2)) * 2 + (nir * 1.5) + (red * 0.5) ) /
(nir + red + 0.5)) * 0.25)) - ((red - 0.125) / (1 - red));
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "LSWI") {
// Land surface water index
out(_,j) = (nir-swir1) / (nir+swir1);
out(_,j) = ifelse((out(_,j) > 1.0) | (out(_,j) < -1.0) | is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "MNDWI") {
// Modified Normalised Difference Water Index
out(_,j) = (green-swir1) / (green+swir1);
out(_,j) = ifelse((out(_,j) > 1.0) | (out(_,j) < -1.0) | is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "MSAVI") {
// Modified soil adjusted vegetation index
out(_,j) = nir + 0.5 - (0.5 * sqrt(pow(2.0 * nir + 1.0, 2) - 8.0 * (nir - (2.0 * red))));
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "MSAVI2") {
// Modified soil adjusted vegetation index 2
out(_,j) = (2.0 * nir + 1.0 - sqrt(pow(2.0 * nir + 1.0, 2) - 8.0 * (nir - red))) / 2.0;
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
/***
else if(indices[j] == "MSI") {
// Moisture stress index
out(_,j) = swir2 / nir;
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
***/
else if(indices[j] == "NBRI"){
// Normalised Burn Ratio Index
out(_,j) = (nir - swir2) / (nir + swir2);
out(_,j) = ifelse((out(_,j) > 1.0) | (out(_,j) < -1.0) | is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "NDVI") {
//Normalized difference vegetation index
out(_,j) = (nir - red) / (nir + red);
out(_,j) = ifelse((out(_,j) > 1.0) | (out(_,j) < -1.0) | is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "NDWI") {
// Normalized difference water index
out(_,j) = (green - nir)/(green + nir);
out(_,j) = ifelse((out(_,j) > 1.0) | (out(_,j) < -1.0) | is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "NRVI") {
// Normalized Ratio Vegetation Index
// Baret and Guyot 1991
NumericVector rvi = red / nir;
out(_,j) = (rvi - 1.0)/(rvi + 1.0);
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "RVI") {
// Ratio Vegetation Index
// Richardson and Wiegand 1977
out(_,j) = red / nir;
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "SATVI"){
// Soil adjusted total vegetation index
out(_,j) = ((swir1 - red) / (swir1 + red + L)) * (1.0 + L) - (swir2 / 2.0);
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "SAVI") {
// Soil adjusted vegetation index
// Huete1988
out(_,j) = (nir - red) * (1.0 + L) / (nir + red + L);
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "SLAVI") {
out(_,j) = nir / (red + swir2);
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "SR") {
// Simple ratio index
out(_,j) = nir / red;
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "TVI") {
// Transformed Vegetation Index
// Deering 1975
out(_,j) = sqrt((nir-red)/(nir+red) + 0.5);
out(_,j) = ifelse((out(_,j) > sqrt(1.5)) | is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "TTVI") {
// Thiam's Transformed Vegetation Index
// Thiam 1997
out(_,j) = sqrt(abs((nir-red)/(nir+red) + 0.5));
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
else if(indices[j] == "WDVI") {
out(_,j) = nir - s * red;
out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j));
}
}
return out;
}