Mercurial > repos > ecology > srs_spectral_indices

comparison prosail-master/R/Lib_SpectralIndices.R @ 0:a8dabbf47e15 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/Marie59/Sentinel_2A/srs_tools commit b32737c1642aa02cc672534e42c5cb4abe0cd3e7
author ecology
date Mon, 09 Jan 2023 13:39:08 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:a8dabbf47e15
1 # ============================================================================== =
2 # prosail
3 # Lib_spectralindices.R
4 # ============================================================================== =
5 # PROGRAMMERS:
6 # Jean-Baptiste FERET <jb.feret@teledetection.fr>
7 # Florian de BOISSIEU <fdeboiss@gmail.com>
8 # Copyright 2019/11 Jean-Baptiste FERET
9 # ============================================================================== =
10 # This Library includes aims at computing spectral indices from reflectance data
11 # ============================================================================== =
12
13 #" This function computes Area under curve for continuum removed reflectances
14 #" See Malenovský et al. (2013) for details
15 #" http://dx.doi.org/10.1016/j.rse.2012.12.015
16 #"
17 #" @param refl RasterBrick, RasterStack or list. Raster bands in the order of sensorbands.
18 #" @param sensorbands numeric. vector containing central wavelength for each spectral band in the image
19 #" @param aucminmax list. wavelengths of lower and upper boundaries ("CRmin" and "CRmax")
20 #" @param reflfactor numeric. multiplying factor used to write reflectance in image (==10000 for S2)
21 #"
22 #" @return aucval raster
23 #" @export
24 auc <- function(refl, sensorbands, aucminmax, reflfactor = 1) {
25
26 aucbands <- list()
27 aucbands[["CRmin"]] <- sensorbands[get_closest_bands(sensorbands, aucminmax[["CRmin"]])]
28 aucbands[["CRmax"]] <- sensorbands[get_closest_bands(sensorbands, aucminmax[["CRmax"]])]
29 bands <- get_closest_bands(sensorbands, aucbands)
30 for (i in bands[["CRmin"]]:bands[["CRmax"]]) {
31 if (is.na(match(i, bands))) {
32 aucbands[[paste("B", i, sep = "")]] <- sensorbands[i]
33 }
34 }
35 # compute continuum removal for all spectral bands
36 cr <- cr_wl(refl = refl, sensorbands = sensorbands,
37 crbands = aucbands, reflfactor = reflfactor)
38
39 wl <- sort(unlist(aucbands), decreasing = FALSE)
40 aucval <- 0.5 * (1 - cr[[1]]) * (wl[2] - wl[1])
41 for (i in 2:length(cr)) {
42 aucval <- aucval + 0.5 * (2 - cr[[i - 1]] - cr[[i]]) * (wl[i + 1] - wl[i])
43 }
44 aucval <- aucval + 0.5 * (1 - cr[[length(cr)]]) * (wl[i + 2] - wl[i + 1])
45 return(aucval)
46 }
47
48 #" This function extracts boundaries to be used to compute continuum from reflectance data
49 #"
50 #" @param refl RasterBrick, RasterStack or list. Raster bands in the order of sensorbands.
51 #" @param sensorbands numeric. vector containing central wavelength for each spectral band in the image
52 #" @param crbands list. list of spectral bands (central wavelength) including CRmin and CRmax
53 #" @param reflfactor numeric. multiplying factor used to write reflectance in image ( == 10000 for S2)
54 #"
55 #" @return crminmax list. list of rasters corresponding to minimum and maximum wavelengths
56 #" @export
57 crbound <- function(refl, sensorbands, crbands, reflfactor = 1) {
58
59 # get closest spectral bands from CR1 and CR2
60 bands <- get_closest_bands(sensorbands, list(crbands[["CRmin"]], crbands[["CRmax"]]))
61 wl <- sensorbands[bands]
62 # get equation for line going from CR1 to CR2
63 crminmax <- readrasterbands(refl = refl, bands = bands, reflfactor = reflfactor)
64 names(crminmax) <- paste("wl_", as.character(wl), sep = "")
65 return(crminmax)
66 }
67
68 #" This function extracts boundaries to be used to compute continuum from reflectance data
69 #"
70 #" @param refl RasterBrick, RasterStack or list. Raster bands in the order of sensorbands.
71 #" @param sensorbands numeric. vector containing central wavelength for each spectral band in the image
72 #" @param crbands list. list of spectral bands (central wavelength) including CRmin and CRmax
73 #" @param reflfactor numeric. multiplying factor used to write reflectance in image ( == 10000 for S2)
74 #"
75 #" @return outlier_iqr numeric. band numbers of original sensor corresponding to S2
76 #" @importFrom progress progress_bar
77 #" @export
78 cr_wl <- function(refl, sensorbands, crbands, reflfactor = 1) {
79
80 # Make sure CRmin and CRmax are correctly defined
81 if (is.na(match("CRmin", names(crbands))) || is.na(match("CRmax", names(crbands)))) {
82 stop("Please define CRmin and CRmax (CRmin<CRmax) as spectral bands in crbands")
83 }
84 if (crbands[["CRmax"]] < crbands[["CRmin"]]) {
85 stop("Please define CRmin < CRmax in crbands")
86 }
87 # extract CRmin and CRmax
88 crminmax <- crbound(refl, sensorbands, crbands, reflfactor = reflfactor)
89 # extract other bands and compute CR
90 crmin <- sensorbands[get_closest_bands(sensorbands, crbands[["CRmin"]])]
91 crmax <- sensorbands[get_closest_bands(sensorbands, crbands[["CRmax"]])]
92 crbands[["CRmin"]] <- NULL
93 crbands[["CRmax"]] <- NULL
94 cr <- list()
95 # initiate progress bar
96 pgbarlength <- length(crbands)
97 pb <- progress_bar$new(
98 format = "Computing continuum removal [:bar] :percent in :elapsedfull, estimated time remaining :eta",
99 total = pgbarlength, clear = FALSE, width = 100)
100 # computation for each band
101 for (band in crbands) {
102 pb$tick()
103 bandrank <- get_closest_bands(sensorbands, band)
104 raster2cr <- readrasterbands(refl = refl, bands = bandrank, reflfactor = reflfactor)
105 cr[[as.character(band)]] <- computecr(wlmin = crmin, wlmax = crmax,
106 wltarget = band, boundaries = crminmax,
107 target = raster2cr)
108 }
109 return(cr)
110 }
111
112 #" This function computes continuum removal value for a spectral band of interest,
113 #" based on lower and upper wavelengths corresponding to boundaries of the continuum
114 #"
115 #" @param wlmin numeric. wavelength of the spectral band corresponding to minimum boundary
116 #" @param wlmax numeric. wavelength of the spectral band corresponding to maximum boundary
117 #" @param wltarget numeric. wavelength of the spectral band for which cr is computed
118 #" @param boundaries list. raster objects corresponding to minimum and maximum wavelengths
119 #" @param target list. raster object corresponding target wavelength
120 #"
121 #" @return cr list. raster object corresponding to continuum removed value
122 #" @export
123 computecr <- function(wlmin, wlmax, wltarget, boundaries, target) {
124
125 cr <- target / (boundaries[[1]] + (wltarget - wlmin) * (boundaries[[2]] - boundaries[[1]]) / (wlmax - wlmin))
126 return(cr)
127 }
128
129 #" this function produces a spectral index from an expression defining a spectral index
130 #"
131 #" @param refl RasterBrick, RasterStack or list. Raster bands in the order of sensorbands.
132 #" @param sensorbands numeric. wavelength in nanometers of the spectral bands of refl.
133 #" @param expressindex character. expression corresponding to the spectral index to compute
134 #" @param listbands list. list of spectral bands defined in the "expressindex" variable
135 #" @param reflfactor numeric. multiplying factor used to write reflectance in image ( == 10000 for S2)
136 #" @param nameindex character. name for the index to be computed, provided in the raster layer
137 #"
138 #" @return numeric. band numbers of original sensor corresponding to S2
139 #" @export
140 spectralindices_fromexpression <- function(refl, sensorbands, expressindex, listbands, reflfactor = 1, nameindex = NULL) {
141
142 # define which bands to be used in the spectral index
143 bands <- get_closest_bands(sensorbands, listbands)
144
145 classraster <- class(refl)
146 if (classraster == "RasterBrick" || classraster == "RasterStack" || classraster == "stars") {
147 # if !reflfactor == 1 then apply a reflectance factor
148 if (classraster == "stars") {
149 refl <- refl[bands]
150 } else {
151 refl <- raster::subset(refl, bands)
152 }
153 if (!reflfactor == 1) {
154 refl <- refl / reflfactor
155 }
156 } else if (is.list(refl)) {
157 refl <- raster::stack(refl[bands]) # checks that all rasters have same crs/extent
158 if (!reflfactor == 1) {
159 refl <- refl / reflfactor
160 }
161 } else {
162 stop("refl is expected to be a RasterStack, RasterBrick, Stars object or a list of rasters")
163 }
164 names(refl) <- gsub(pattern = "B", replacement = "Band", x = names(bands))
165
166 nbbands <- unique(as.numeric(gsub(pattern = "B",
167 replacement = "",
168 x = unlist(regmatches(expressindex,
169 gregexpr("B[[:digit:]] + ",
170 expressindex))))))
171 sortband <- sort(nbbands, index.return = TRUE, decreasing = TRUE)
172 matches <- unique(unlist(regmatches(expressindex, gregexpr("B[[:digit:]] + ", expressindex))))[sortband$ix]
173 replaces <- paste("refl[['Band", gsub(pattern = "B", replacement = "", x = matches), "']]", sep = "")
174 expressindex_final <- expressindex
175 for (bb in 1:seq_along(matches)) {
176 expressindex_final <- gsub(pattern = matches[bb], replacement = replaces[bb], x = expressindex_final)
177 }
178 si <- eval(parse(text = expressindex_final))
179 if (!is.null(nameindex)) {
180 names(si) <- nameindex
181 }
182 return(si)
183 }
184
185 #" this function aims at computing spectral indices from Sensor reflectance data in raster object
186 #" it computes the spectral indices based on their computation with Sentinel-2
187 #" and assumes that the bands of the S2 data follow this order
188 #" wavelength = {496.6, 560.0, 664.5, 703.9, 740.2, 782.5, 835.1, 864.8, 1613.7, 2202.4}
189 #" Full description of the indices can be found here:
190 #" https://www.sentinel-hub.com/eotaxonomy/indices
191 #"
192 #" @param refl RasterBrick, RasterStack or list. Raster bands in the order of sensorbands.
193 #" @param sensorbands numeric. wavelength in nanometers of the spectral bands of refl.
194 #" @param sel_indices list. list of spectral indices to be computed
195 #" @param stackout logical. If TRUE returns a stack, otherwise a list of rasters.
196 #" @param reflfactor numeric. multiplying factor used to write reflectance in image ( == 10000 for S2)
197 #"
198 #" @return list. includes
199 #" - spectralindices: List of spectral indices computed from the reflectance initially provided
200 #" - listindices: list of spectral indices computable with the function
201 #" @importFrom methods is
202 #" @importFrom raster stack brick
203 #" @export
204
205 computespectralindices_raster <- function(refl, sensorbands, sel_indices = "ALL", stackout = TRUE, reflfactor = 1) {
206
207 s2bands <- c("B2" = 496.6, "B3" = 560.0, "B4" = 664.5, "B5" = 703.9, "B6" = 740.2,
208 "B7" = 782.5, "B8" = 835.1, "B8A" = 864.8, "B11" = 1613.7, "B12" = 2202.4)
209
210 spectralindices <- list()
211 sen2s2 <- get_closest_bands(sensorbands, s2bands)
212 classraster <- class(refl)
213 if (classraster == "RasterBrick" || classraster == "RasterStack" || classraster == "stars") {
214 # if !reflfactor == 1 then apply a reflectance factor
215 if (classraster == "stars") {
216 refl <- refl[sen2s2]
217 } else {
218 refl <- raster::subset(refl, sen2s2)
219 }
220 if (!reflfactor == 1) {
221 refl <- refl / reflfactor
222 }
223 } else if (is.list(refl)) {
224 refl <- raster::stack(refl[sen2s2]) # checks that all rasters have same crs/extent
225 if (!reflfactor == 1) {
226 refl <- refl / reflfactor
227 }
228 } else {
229 stop("refl is expected to be a RasterStack, RasterBrick, Stars object or a list of rasters")
230 }
231 names(refl) <- names(sen2s2)
232
233 indexall <- list()
234
235 # inelegant but meeeeh
236 listindices <- list("ARI1", "ARI2", "ARVI", "BAI", "BAIS2", "CCCI", "CHL_RE", "CRI1", "CRI2", "EVI", "EVI2",
237 "GRVI1", "GNDVI", "IRECI", "LAI_SAVI", "MCARI", "mNDVI705", "MSAVI2",
238 "MSI", "mSR705", "MTCI", "nBR_RAW", "NDI_45", "NDII", "NDSI", "NDVI", "NDVI_G",
239 "NDVI705", "NDWI1", "NDWI2", "PSRI", "PSRI_NIR", "RE_NDVI", "RE_NDWI", "S2REP",
240 "SAVI", "SIPI", "SR", "CR_SWIR")
241 if (sel_indices[1] == "ALL") {
242 sel_indices <- listindices
243 }
244 if ("ARI1" %in% sel_indices) {
245 ari1 <- (1 / refl[["B3"]]) - (1 / refl[["B5"]])
246 spectralindices$ARI1 <- ari1
247 }
248 if ("ARI2" %in% sel_indices) {
249 ari2 <- (refl[["B8"]] / refl[["B2"]]) - (refl[["B8"]] / refl[["B3"]])
250 spectralindices$ARI2 <- ari2
251 }
252 if ("ARVI" %in% sel_indices) {
253 arvi <- (refl[["B8"]] - (2 * refl[["B4"]] - refl[["B2"]])) / (refl[["B8"]] + (2 * refl[["B4"]] - refl[["B2"]]))
254 spectralindices$ARVI <- arvi
255 }
256 if ("BAI" %in% sel_indices) {
257 bai <- (1 / ((0.1 - refl[["B4"]])**2 + (0.06 - refl[["B8"]])**2))
258 spectralindices$BAI <- bai
259 }
260 if ("BAIS2" %in% sel_indices) {
261 bais2 <- (1 - ((refl[["B6"]] * refl[["B7"]] * refl[["B8A"]]) / refl[["B4"]])**0.5) * ((refl[["B12"]] - refl[["B8A"]]) / ((refl[["B12"]] + refl[["B8A"]])**0.5) + 1)
262 spectralindices$BAIS2 <- bais2
263 }
264 if ("CCCI" %in% sel_indices) {
265 ccci <- ((refl[["B8"]] - refl[["B5"]]) / (refl[["B8"]] + refl[["B5"]])) / ((refl[["B8"]] - refl[["B4"]]) / (refl[["B8"]] + refl[["B4"]]))
266 spectralindices$CCCI <- ccci
267 }
268 if ("CHL_RE" %in% sel_indices) {
269 chl_re <- refl[["B5"]] / refl[["B8"]]
270 spectralindices$CHL_RE <- chl_re
271 }
272 if ("CRI1" %in% sel_indices) {
273 cri1 <- (1 / refl[["B2"]]) - (1 / refl[["B3"]])
274 spectralindices$CRI1 <- cri1
275 }
276 if ("CRI2" %in% sel_indices) {
277 cri2 <- (1 / refl[["B2"]]) - (1 / refl[["B5"]])
278 spectralindices$CRI2 <- cri2
279 }
280 if ("EVI" %in% sel_indices) {
281 evi <- 2.5 * (refl[["B8"]] - refl[["B4"]]) / ((refl[["B8"]] + 6 * refl[["B4"]] - 7.5 * refl[["B2"]] + 1))
282 spectralindices$EVI <- evi
283 }
284 if ("EVI2" %in% sel_indices) {
285 evi2 <- 2.5 * (refl[["B8"]] - refl[["B4"]]) / (refl[["B8"]] + 2.4 * refl[["B4"]] + 1)
286 spectralindices$EVI2 <- evi2
287 }
288 if ("GRVI1" %in% sel_indices) {
289 grvi1 <- (refl[["B4"]] - refl[["B3"]]) / (refl[["B4"]] + refl[["B3"]])
290 spectralindices$GRVI1 <- grvi1
291 }
292 if ("GNDVI" %in% sel_indices) {
293 gndvi <- (refl[["B8"]] - refl[["B3"]]) / (refl[["B8"]] + refl[["B3"]])
294 spectralindices$GNDVI <- gndvi
295 }
296 if ("IRECI" %in% sel_indices) {
297 ireci <- (refl[["B7"]] - refl[["B4"]]) * (refl[["B6"]] / (refl[["B5"]]))
298 spectralindices$IRECI <- ireci
299 }
300 if ("LAI_SAVI" %in% sel_indices) {
301 lai_savi <- - log(0.371 + 1.5 * (refl[["B8"]] - refl[["B4"]]) / (refl[["B8"]] + refl[["B4"]] + 0.5)) / 2.4
302 spectralindices$LAI_SAVI <- lai_savi
303 }
304 if ("MCARI" %in% sel_indices) {
305 mcari <- (1 - 0.2 * (refl[["B5"]] - refl[["B3"]]) / (refl[["B5"]] - refl[["B4"]]))
306 spectralindices$MCARI <- mcari
307 }
308 if ("mNDVI705" %in% sel_indices) {
309 mndvi705 <- (refl[["B6"]] - refl[["B5"]]) / (refl[["B6"]] + refl[["B5"]] - 2 * refl[["B2"]])
310 spectralindices$mNDVI705 <- mndvi705
311 }
312 if ("MSAVI2" %in% sel_indices) {
313 msavi2 <- ((refl[["B8"]] + 1) - 0.5 * sqrt(((2 * refl[["B8"]]) - 1)**2 + 8 * refl[["B4"]]))
314 spectralindices$MSAVI2 <- msavi2
315 }
316 if ("MSI" %in% sel_indices) {
317 msi <- refl[["B11"]] / refl[["B8A"]]
318 spectralindices$MSI <- msi
319 }
320 if ("mSR705" %in% sel_indices) {
321 msr705 <- (refl[["B6"]] - refl[["B2"]]) / (refl[["B5"]] - refl[["B2"]])
322 spectralindices$mSR705 <- msr705
323 }
324 if ("MTCI" %in% sel_indices) {
325 mtci <- (refl[["B6"]] - refl[["B5"]]) / (refl[["B5"]] + refl[["B4"]])
326 spectralindices$MTCI <- mtci
327 }
328 if ("nBR_RAW" %in% sel_indices) {
329 nbr_raw <- (refl[["B8"]] - refl[["B12"]]) / (refl[["B8"]] + refl[["B12"]])
330 spectralindices$nBR_RAW <- nbr_raw
331 }
332 if ("NDI_45" %in% sel_indices) {
333 ndi_45 <- (refl[["B5"]] - refl[["B4"]]) / (refl[["B5"]] + refl[["B4"]])
334 spectralindices$NDI_45 <- ndi_45
335 }
336 if ("NDII" %in% sel_indices) {
337 ndii <- (refl[["B8A"]] - refl[["B11"]]) / (refl[["B8A"]] + refl[["B11"]])
338 spectralindices$NDII <- ndii
339 }
340 if ("NDSI" %in% sel_indices) {
341 ndsi <- (refl[["B3"]] - refl[["B11"]]) / (refl[["B3"]] + refl[["B11"]])
342 spectralindices$NDSI <- ndsi
343 }
344 if ("NDVI" %in% sel_indices) {
345 ndvi <- (refl[["B8"]] - refl[["B4"]]) / (refl[["B8"]] + refl[["B4"]])
346 spectralindices$NDVI <- ndvi
347 }
348 if ("NDVI_G" %in% sel_indices) {
349 ndvi_g <- refl[["B3"]] * (refl[["B8"]] - refl[["B4"]]) / (refl[["B8"]] + refl[["B4"]])
350 spectralindices$NDVI_G <- ndvi_g
351 }
352 if ("NDVI705" %in% sel_indices) {
353 ndvi705 <- (refl[["B6"]] - refl[["B5"]]) / (refl[["B6"]] + refl[["B5"]])
354 spectralindices$NDVI705 <- ndvi705
355 }
356 if ("NDWI" %in% sel_indices) {
357 ndwi <- (refl[["B3"]] - refl[["B8"]]) / (refl[["B3"]] + refl[["B8"]])
358 spectralindices$NDWI <- ndwi
359 }
360 if ("NDWI1" %in% sel_indices) {
361 ndwi1 <- (refl[["B8A"]] - refl[["B11"]]) / (refl[["B8A"]] + refl[["B11"]])
362 spectralindices$NDWI1 <- ndwi1
363 }
364 if ("NDWI2" %in% sel_indices) {
365 ndwi2 <- (refl[["B8A"]] - refl[["B12"]]) / (refl[["B8A"]] + refl[["B12"]])
366 spectralindices$NDWI2 <- ndwi2
367 }
368 if ("PSRI" %in% sel_indices) {
369 psri <- (refl[["B4"]] - refl[["B2"]]) / (refl[["B5"]])
370 spectralindices$PSRI <- psri
371 }
372 if ("PSRI_NIR" %in% sel_indices) {
373 psri_nir <- (refl[["B4"]] - refl[["B2"]]) / (refl[["B8"]])
374 spectralindices$PSRI_NIR <- psri_nir
375 }
376 if ("RE_NDVI" %in% sel_indices) {
377 re_ndvi <- (refl[["B8"]] - refl[["B6"]]) / (refl[["B8"]] + refl[["B6"]])
378 spectralindices$RE_NDVI <- re_ndvi
379 }
380 if ("RE_NDWI" %in% sel_indices) {
381 re_ndwi <- (refl[["B4"]] - refl[["B6"]]) / (refl[["B4"]] + refl[["B6"]])
382 spectralindices$RE_NDWI <- re_ndwi
383 }
384 if ("S2REP" %in% sel_indices) {
385 s2rep <- 705 + 35 * (0.5 * (refl[["B8"]] + refl[["B5"]]) - refl[["B6"]]) / (refl[["B7"]] - refl[["B6"]])
386 spectralindices$S2REP <- s2rep
387 }
388 if ("SAVI" %in% sel_indices) {
389 savi <- 1.5 * (refl[["B8"]] - refl[["B5"]]) / (refl[["B8"]] + refl[["B5"]] + 0.5)
390 spectralindices$SAVI <- savi
391 }
392 if ("SIPI" %in% sel_indices) {
393 sipi <- (refl[["B8"]] - refl[["B2"]]) / (refl[["B8"]] - refl[["B4"]])
394 spectralindices$SIPI <- sipi
395 }
396 if ("SR" %in% sel_indices) {
397 sr <- refl[["B8"]] / refl[["B4"]]
398 spectralindices$SR <- sr
399 }
400 if ("TCARI" %in% sel_indices) {
401 sr <- refl[["B8"]] / refl[["B4"]]
402 spectralindices$SR <- sr
403 }
404 if ("CR_SWIR" %in% sel_indices) {
405 cr_swir <- refl[["B11"]] / (refl[["B8A"]] + (s2bands["B11"] - s2bands["B8A"]) * (refl[["B12"]] - refl[["B8A"]]) / (s2bands["B12"] - s2bands["B8A"]))
406 spectralindices$CR_SWIR <- cr_swir
407 }
408
409 if (stackout)
410 spectralindices <- raster::stack(spectralindices)
411
412 res <- list("spectralindices" = spectralindices, "listindices" = listindices)
413 return(res)
414 }
415
416 #" this function aims at computing spectral indices from Sensor reflectance data.
417 #" it computes the spectral indices based on their computation with Sentinel-2
418 #" and assumes that the bands of the S2 data follow this order
419 #" wavelength = {496.6, 560.0, 664.5, 703.9, 740.2, 782.5, 835.1, 864.8, 1613.7, 2202.4}
420 #" Full description of the indices can be found here:
421 #" https://www.sentinel-hub.com/eotaxonomy/indices
422 #"
423 #" @param refl numeric. reflectance dataset defined in matrix
424 #" @param sel_indices list. list of spectral indices to be computed
425 #" @param sensorbands numeric. wavelength of the spectral bands corresponding to the spectral index
426 #"
427 #" @return list. includes
428 #" - spectralindices: List of spectral indices computed from the reflectance initially provided
429 #" - listindices: list of spectral indices computable with the function
430 #" @export
431
432 computespectralindices_hs <- function(refl, sensorbands, sel_indices = "ALL") {
433
434 spectralindices <- list()
435 s2bands <- data.frame("B2" = 496.6, "B3" = 560.0, "B4" = 664.5, "B5" = 703.9, "B6" = 740.2,
436 "B7" = 782.5, "B8" = 835.1, "B8A" = 864.8, "B11" = 1613.7, "B12" = 2202.4)
437
438 sen2s2 <- get_closest_bands(sensorbands, s2bands)
439 indexall <- list()
440 # set zero vaues to >0 in order to avoid problems
441 selzero <- which(refl == 0)
442 refl[selzero] <- 0.005
443 if (dim(refl)[1] == length(sensorbands)) {
444 refl <- t(refl)
445 }
446
447 # inelegant but meeeeh
448 listindices <- list("ARI1", "ARI2", "ARVI", "BAI", "BAIS2", "CHL_RE", "CRI1", "CRI2", "EVI", "EVI2",
449 "GRVI1", "GNDVI", "IRECI", "LAI_SAVI", "MCARI", "mNDVI705", "MSAVI2",
450 "MSI", "mSR705", "MTCI", "nBR_RAW", "NDI_45", "NDII", "NDVI", "NDVI_G",
451 "NDVI705", "NDWI1", "NDWI2", "PSRI", "PSRI_NIR", "RE_NDVI", "RE_NDWI", "S2REP",
452 "SAVI", "SIPI", "SR", "CR_SWIR")
453 if (sel_indices == "ALL") {
454 sel_indices <- listindices
455 }
456 if ("ARI1" %in% sel_indices) {
457 ari1 <- (1 / refl[, sen2s2[["B3"]]]) - (1 / refl[, sen2s2[["B5"]]])
458 spectralindices$ARI1 <- ari1
459 }
460 if ("ARI2" %in% sel_indices) {
461 ari2 <- (refl[, sen2s2[["B8"]]] / refl[, sen2s2[["B2"]]]) - (refl[, sen2s2[["B8"]]] / refl[, sen2s2[["B3"]]])
462 spectralindices$ARI2 <- ari2
463 }
464 if ("ARVI" %in% sel_indices) {
465 arvi <- (refl[, sen2s2[["B8"]]] - (2 * refl[, sen2s2[["B4"]]] - refl[, sen2s2[["B2"]]])) / (refl[, sen2s2[["B8"]]] + (2 * refl[, sen2s2[["B4"]]] - refl[, sen2s2[["B2"]]]))
466 spectralindices$ARVI <- arvi
467 }
468 if ("BAI" %in% sel_indices) {
469 bai <- (1 / ((0.1 - refl[, sen2s2[["B4"]]])**2 + (0.06 - refl[, sen2s2[["B8"]]])**2))
470 spectralindices$BAI <- bai
471 }
472 if ("BAIS2" %in% sel_indices) {
473 bais2 <- (1 - ((refl[, sen2s2[["B6"]]] * refl[, sen2s2[["B7"]]] * refl[, sen2s2[["B8A"]]]) / refl[, sen2s2[["B4"]]])**0.5) * ((refl[, sen2s2[["B12"]]] - refl[, sen2s2[["B8A"]]]) / ((refl[, sen2s2[["B12"]]] + refl[, sen2s2[["B8A"]]])**0.5) + 1)
474 spectralindices$BAIS2 <- bais2
475 }
476 if ("CCCI" %in% sel_indices) {
477 ccci <- ((refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B5"]]]) / (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B5"]]])) / ((refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B4"]]]) / (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B4"]]]))
478 spectralindices$CCCI <- ccci
479 }
480 if ("CHL_RE" %in% sel_indices) {
481 chl_re <- refl[, sen2s2[["B5"]]] / refl[, sen2s2[["B8"]]]
482 spectralindices$CHL_RE <- chl_re
483 }
484 if ("CRI1" %in% sel_indices) {
485 cri1 <- (1 / refl[, sen2s2[["B2"]]]) - (1 / refl[, sen2s2[["B3"]]])
486 spectralindices$CRI1 <- cri1
487 }
488 if ("CRI2" %in% sel_indices) {
489 cri2 <- (1 / refl[, sen2s2[["B2"]]]) - (1 / refl[, sen2s2[["B5"]]])
490 spectralindices$CRI2 <- cri2
491 }
492 if ("EVI" %in% sel_indices) {
493 evi <- 2.5 * (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B4"]]]) / ((refl[, sen2s2[["B8"]]] + 6 * refl[, sen2s2[["B4"]]] - 7.5 * refl[, sen2s2[["B2"]]] + 1))
494 spectralindices$EVI <- evi
495 }
496 if ("EVI2" %in% sel_indices) {
497 evi2 <- 2.5 * (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B4"]]]) / (refl[, sen2s2[["B8"]]] + 2.4 * refl[, sen2s2[["B4"]]] + 1)
498 spectralindices$EVI2 <- evi2
499 }
500 if ("GRVI1" %in% sel_indices) {
501 grvi1 <- (refl[, sen2s2[["B4"]]] - refl[, sen2s2[["B3"]]]) / (refl[, sen2s2[["B4"]]] + refl[, sen2s2[["B3"]]])
502 spectralindices$GRVI1 <- grvi1
503 }
504 if ("GNDVI" %in% sel_indices) {
505 gndvi <- (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B3"]]]) / (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B3"]]])
506 spectralindices$GNDVI <- gndvi
507 }
508 if ("IRECI" %in% sel_indices) {
509 ireci <- (refl[, sen2s2[["B7"]]] - refl[, sen2s2[["B4"]]]) * (refl[, sen2s2[["B6"]]] / (refl[, sen2s2[["B5"]]]))
510 spectralindices$IRECI <- ireci
511 }
512 if ("LAI_SAVI" %in% sel_indices) {
513 lai_savi <- - log(0.371 + 1.5 * (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B4"]]]) / (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B4"]]] + 0.5)) / 2.4
514 spectralindices$LAI_SAVI <- lai_savi
515 }
516 if ("MCARI" %in% sel_indices) {
517 mcari <- (1 - 0.2 * (refl[, sen2s2[["B5"]]] - refl[, sen2s2[["B3"]]]) / (refl[, sen2s2[["B5"]]] - refl[, sen2s2[["B4"]]]))
518 spectralindices$MCARI <- mcari
519 }
520 if ("mNDVI705" %in% sel_indices) {
521 mndvi705 <- (refl[, sen2s2[["B6"]]] - refl[, sen2s2[["B5"]]]) / (refl[, sen2s2[["B6"]]] + refl[, sen2s2[["B5"]]] - 2 * refl[, sen2s2[["B2"]]])
522 spectralindices$mNDVI705 <- mndvi705
523 }
524 if ("MSAVI2" %in% sel_indices) {
525 msavi2 <- ((refl[, sen2s2[["B8"]]] + 1) - 0.5 * sqrt(((2 * refl[, sen2s2[["B8"]]]) - 1)**2 + 8 * refl[, sen2s2[["B4"]]]))
526 spectralindices$MSAVI2 <- msavi2
527 }
528 if ("MSI" %in% sel_indices) {
529 msi <- refl[, sen2s2[["B11"]]] / refl[, sen2s2[["B8"]]]
530 spectralindices$MSI <- msi
531 }
532 if ("mSR705" %in% sel_indices) {
533 msr705 <- (refl[, sen2s2[["B6"]]] - refl[, sen2s2[["B2"]]]) / (refl[, sen2s2[["B5"]]] - refl[, sen2s2[["B2"]]])
534 spectralindices$mSR705 <- msr705
535 }
536 if ("MTCI" %in% sel_indices) {
537 mtci <- (refl[, sen2s2[["B6"]]] - refl[, sen2s2[["B5"]]]) / (refl[, sen2s2[["B5"]]] + refl[, sen2s2[["B4"]]])
538 spectralindices$MTCI <- mtci
539 }
540 if ("nBR_RAW" %in% sel_indices) {
541 nbr_raw <- (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B12"]]]) / (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B12"]]])
542 spectralindices$nBR_RAW <- nbr_raw
543 }
544 if ("NDI_45" %in% sel_indices) {
545 ndi_45 <- (refl[, sen2s2[["B5"]]] - refl[, sen2s2[["B4"]]]) / (refl[, sen2s2[["B5"]]] + refl[, sen2s2[["B4"]]])
546 spectralindices$NDI_45 <- ndi_45
547 }
548 if ("NDII" %in% sel_indices) {
549 ndii <- (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B11"]]]) / (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B11"]]])
550 spectralindices$NDII <- ndii
551 }
552 if ("NDSI" %in% sel_indices) {
553 ndisi <- (refl[, sen2s2[["B3"]]] - refl[, sen2s2[["B11"]]]) / (refl[, sen2s2[["B3"]]] + refl[, sen2s2[["B11"]]])
554 spectralindices$NDSI <- ndsi
555 }
556 if ("NDVI" %in% sel_indices) {
557 ndvi <- (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B4"]]]) / (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B4"]]])
558 spectralindices$NDVI <- ndvi
559 }
560 if ("NDVI_G" %in% sel_indices) {
561 ndvi_g <- refl[, sen2s2[["B3"]]] * (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B4"]]]) / (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B4"]]])
562 spectralindices$NDVI_G <- ndvi_g
563 }
564 if ("NDVI705" %in% sel_indices) {
565 ndvi705 <- (refl[, sen2s2[["B6"]]] - refl[, sen2s2[["B5"]]]) / (refl[, sen2s2[["B6"]]] + refl[, sen2s2[["B5"]]])
566 spectralindices$NDVI705 <- ndvi705
567 }
568 if ("NDWI1" %in% sel_indices) {
569 ndwi1 <- (refl[, sen2s2[["B8A"]]] - refl[, sen2s2[["B11"]]]) / (refl[, sen2s2[["B8A"]]] + refl[, sen2s2[["B11"]]])
570 spectralindices$NDWI1 <- ndwi1
571 }
572 if ("NDWI2" %in% sel_indices) {
573 ndwi2 <- (refl[, sen2s2[["B8A"]]] - refl[, sen2s2[["B12"]]]) / (refl[, sen2s2[["B8A"]]] + refl[, sen2s2[["B12"]]])
574 spectralindices$NDWI2 <- ndwi2
575 }
576 if ("PSRI" %in% sel_indices) {
577 psri <- (refl[, sen2s2[["B4"]]] - refl[, sen2s2[["B2"]]]) / (refl[, sen2s2[["B5"]]])
578 spectralindices$PSRI <- psri
579 }
580 if ("PSRI_NIR" %in% sel_indices) {
581 psri_nir <- (refl[, sen2s2[["B4"]]] - refl[, sen2s2[["B2"]]]) / (refl[, sen2s2[["B8"]]])
582 spectralindices$PSRI_NIR <- psri_nir
583 }
584 if ("RE_NDVI" %in% sel_indices) {
585 re_ndvi <- (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B6"]]]) / (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B6"]]])
586 spectralindices$RE_NDVI <- re_ndvi
587 }
588 if ("RE_NDWI" %in% sel_indices) {
589 re_ndwi <- (refl[, sen2s2[["B4"]]] - refl[, sen2s2[["B6"]]]) / (refl[, sen2s2[["B4"]]] + refl[, sen2s2[["B6"]]])
590 spectralindices$RE_NDWI <- re_ndwi
591 }
592 if ("S2REP" %in% sel_indices) {
593 s2rep <- 705 + 35 * (0.5 * (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B5"]]]) - refl[, sen2s2[["B6"]]]) / (refl[, sen2s2[["B7"]]] - refl[, sen2s2[["B6"]]])
594 spectralindices$S2REP <- s2rep
595 }
596 if ("SAVI" %in% sel_indices) {
597 savi <- 1.5 * (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B5"]]]) / (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B5"]]] + 0.5)
598 spectralindices$SAVI <- savi
599 }
600 if ("SIPI" %in% sel_indices) {
601 sipi <- (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B2"]]]) / (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B4"]]])
602 spectralindices$SIPI <- sipi
603 }
604 if ("SR" %in% sel_indices) {
605 sr <- refl[, sen2s2[["B8"]]] / refl[, sen2s2[["B4"]]]
606 spectralindices$SR <- sr
607 }
608 if ("CR_SWIR" %in% sel_indices) {
609 cr_swir <- refl[, sen2s2[["B11"]]] / (refl[, sen2s2[["B8A"]]] + (s2bands$B11 - s2bands$B8A) * (refl[, sen2s2[["B12"]]] - refl[, sen2s2[["B8A"]]]) / (s2bands$B12 - s2bands$B8A))
610 spectralindices$CR_SWIR <- cr_swir
611 }
612 res <- list("spectralindices" = spectralindices, "listindices" = listindices)
613 return(res)
614 }
615
616 #" this function identifies the bands of a given sensor with closest match to its spectral characteristics
617 #"
618 #" @param sensorbands numeric. wavelength in nanometer of the sensor of interest
619 #" @param listbands numeric or list. Named vector or list of spectral bands corresponding to sensor
620 #"
621 #" @return numeric. band numbers of original sensor
622 #" @export
623 get_closest_bands <- function(sensorbands, listbands) {
624 sapply(listbands, function(x) {
625 b <- which.min(abs(sensorbands - x))
626 names(b) <- ""
627 b
628 })
629 }
630
631 #" This function computes interquartile range (IQR) criterion, which can be used
632 #" as a criterion for outlier detection
633 #"
634 #" @param distval numeric. vector of distribution of values
635 #" @param weightirq numeric. weighting factor appplied to IRQ to define lower and upper boudaries for outliers
636 #"
637 #" @return outlier_iqr numeric. band numbers of original sensor corresponding to S2
638 #" @importFrom stats IQR quantile
639 #" @export
640 iqr_outliers <- function(distval, weightirq = 1.5) {
641 iqr <- IQR(distval, na.rm = TRUE)
642 range_iqr <- c(quantile(distval, 0.25, na.rm = TRUE), quantile(distval, 0.75, na.rm = TRUE))
643 outlier_iqr <- c(range_iqr[1] - weightirq * iqr, range_iqr[2] + weightirq * iqr)
644 return(outlier_iqr)
645 }
646
647 #" This function selects bands from a raster or stars object
648 #"
649 #" @param refl RasterBrick, RasterStack or list. Raster bands in the order of sensorbands.
650 #" @param bands numeric. rank of bands to be read in refl
651 #" @param reflfactor numeric. multiplying factor used to write reflectance in image ( == 10000 for S2)
652 #"
653 #" @return robj list. R object (default is raster, stars if refl is stars object)
654 #" @importFrom raster subset stack
655 #" @export
656 readrasterbands <- function(refl, bands, reflfactor = 1) {
657
658 # get equation for line going from CR1 to CR2
659 classraster <- class(refl)
660 if (classraster == "RasterBrick" || classraster == "RasterStack" || classraster == "stars") {
661 # if !reflfactor == 1 then apply a reflectance factor
662 if (classraster == "stars") {
663 robj <- refl[bands]
664 } else {
665 robj <- raster::subset(refl, bands)
666 }
667 if (!reflfactor == 1) {
668 robj <- robj / reflfactor
669 }
670 } else if (is.list(refl)) {
671 robj <- raster::stack(refl[bands]) # checks that all rasters have same crs/extent
672 if (!reflfactor == 1) {
673 robj <- robj / reflfactor
674 }
675 } else {
676 stop("refl is expected to be a RasterStack, RasterBrick, Stars object or a list of rasters")
677 }
678 return(robj)
679 }