Spectral indices time series
The function calcIndices
is used to calculate spectral indices from multispectral data. The spectral indices calculation is based on the RStoolbox
package. The list of supported spectral indices and their formula is described in the documentation of the spectralIndices
function of the RStoolbox
package (available here. The layer index corresponding to the blue
, green
, red
, nir
, swir1
, swir3
bands need to be provided depending on the indices calculated.
Tasseled Cap brightness (TCB), greenness (TCG) and wetness (TCW) can also be calculated. Calculation of tasseled cap indices is based on the function tasseledCap
from the RSToolbox
package. The name of the satellite that acquired the data has to be provided under sat
and the bands have to be ordered in a specific order as explained in the documentation of RStoolbox::tasseledCap
(bands 1, 2, 3, 4, 5, 7 for Landsat5TM).
For this example, we calculate the Tasseled Cap Brightness (TCB), Greenness (TCG), Wetness (TCW) and NDVI. NDVI requires red
and nir
bands only.
indices_list <- c("NDVI","TCB","TCW","TCG")
# Example for one year
VI_2006 <- calcIndices(spectral_2006,
indices = indices_list,
sat="Landsat5TM",
red=3,
nir=4)
plot(VI_2006[["TCG"]])
To process multiple years at once, all RasterStack objects (one per year) can be placed in a list.
# List all filenames in time-series order: 2006 to 2008
spectral_ts_files <- list(system.file("extdata/vignette/Landsat_BAP/Landsat_BAP_2006.tif",package = "foster"),
system.file("extdata/vignette/Landsat_BAP/Landsat_BAP_2007.tif",package = "foster"),
system.file("extdata/vignette/Landsat_BAP/Landsat_BAP_2008.tif",package = "foster"))
# Open all Landsat images and place them in a list
spectral_ts <- lapply(spectral_ts_files, stack)
spectral_ts
#> [[1]]
#> class : RasterStack
#> dimensions : 133, 134, 17822, 6 (nrow, ncol, ncell, nlayers)
#> resolution : 30, 30 (x, y)
#> extent : 584955, 588975, 5811729, 5815719 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
#> names : Landsat_BAP_2006.1, Landsat_BAP_2006.2, Landsat_BAP_2006.3, Landsat_BAP_2006.4, Landsat_BAP_2006.5, Landsat_BAP_2006.6
#> min values : 7, 89, 59, 42, 13, 18
#> max values : 988, 1133, 1421, 4364, 3285, 2258
#>
#>
#> [[2]]
#> class : RasterStack
#> dimensions : 133, 134, 17822, 6 (nrow, ncol, ncell, nlayers)
#> resolution : 30, 30 (x, y)
#> extent : 584955, 588975, 5811729, 5815719 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
#> names : Landsat_BAP_2007.1, Landsat_BAP_2007.2, Landsat_BAP_2007.3, Landsat_BAP_2007.4, Landsat_BAP_2007.5, Landsat_BAP_2007.6
#> min values : 33, 93, 64, 125, 20, 17
#> max values : 685, 1065, 1211, 4494, 3108, 2135
#>
#>
#> [[3]]
#> class : RasterStack
#> dimensions : 133, 134, 17822, 6 (nrow, ncol, ncell, nlayers)
#> resolution : 30, 30 (x, y)
#> extent : 584955, 588975, 5811729, 5815719 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
#> names : Landsat_BAP_2008.1, Landsat_BAP_2008.2, Landsat_BAP_2008.3, Landsat_BAP_2008.4, Landsat_BAP_2008.5, Landsat_BAP_2008.6
#> min values : 19, 98, 39, 131, 17, 16
#> max values : 629, 965, 1082, 4251, 2845, 1954
VI_ts <- calcIndices(spectral_ts,
indices = indices_list,
sat="Landsat5TM",
red=3,
nir=4)
# The output is a list with VI for each year
VI_ts
#> [[1]]
#> class : RasterStack
#> dimensions : 133, 134, 17822, 4 (nrow, ncol, ncell, nlayers)
#> resolution : 30, 30 (x, y)
#> extent : 584955, 588975, 5811729, 5815719 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
#> names : NDVI, TCB, TCW, TCG
#> min values : -0.3777778, 158.6998000, 23.4080000, -76.8098000
#> max values : 0.9267045, 4436.6082000, 2004.1029000, 2850.4314000
#>
#>
#> [[2]]
#> class : RasterStack
#> dimensions : 133, 134, 17822, 4 (nrow, ncol, ncell, nlayers)
#> resolution : 30, 30 (x, y)
#> extent : 584955, 588975, 5811729, 5815719 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
#> names : NDVI, TCB, TCW, TCG
#> min values : 0.1261261, 198.2568000, 47.2045000, -7.9697000
#> max values : 0.900304, 4265.392400, 2022.956900, 2845.887100
#>
#>
#> [[3]]
#> class : RasterStack
#> dimensions : 133, 134, 17822, 4 (nrow, ncol, ncell, nlayers)
#> resolution : 30, 30 (x, y)
#> extent : 584955, 588975, 5811729, 5815719 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
#> names : NDVI, TCB, TCW, TCG
#> min values : 0.1155116, 206.6494000, 21.4235000, 2.0134000
#> max values : 0.9105449, 3826.3472000, 1813.1592000, 2978.8278000
We can visualize the vegetation indices of the first year of the time series
Once indices have been calculated we can smooth them with a 3 x 3 mean filter, as it was done with the resampled ALS metrics maps.
VI_ts_smooth <- focalMultiBand(VI_ts,
w=filt,
fun=mean,
na.rm=TRUE,
pad = TRUE,
keepNA = TRUE)
We can visualize the smoothed vegetation indices of the first year of the time series
Spectral indices temporal summary metrics
We will now calculate temporal summary metrics of the vegetation indices time series.
The function temporalMetrics
requires the name of a function that returns summary metrics of a numeric vector (e.g mean, standard deviation). The default function used by temporalMetrics
returns the median, IQR and Theil-Sen slope. However, the user can also define another function that returns metrics specific to its needs. For illustration, we create the function funSummary
that returns the same metrics as the default function. Ideally, gap-free composite multispectral images should be used but when creating your own function you can handle the cases where NA values might occur in the time series (by setting the na.rm
argument to TRUE
or FALSE
).
funSummary <- function(x){
c(
median = median(x,na.rm=TRUE),
IQR = IQR(x,na.rm=TRUE),
slope = theilSen(x)
)
}
We can calculate the temporal summary of the 3 year time series of smoothed vegetation indices created above (VI_ts_smooth
)
VI_ts_metrics <- temporalMetrics(VI_ts,
metrics="funSummary")
#> Working on NDVI
#> Working on TCB
#> Working on TCW
#> Working on TCG
The output is a RasterStack where each band is a temporal summary metric of a vegetation index. In this example, there are 12 bands (4 vegetation indices, 3 temporal summary metrics).
VI_ts_metrics
#> class : RasterStack
#> dimensions : 133, 134, 17822, 12 (nrow, ncol, ncell, nlayers)
#> resolution : 30, 30 (x, y)
#> extent : 584955, 588975, 5811729, 5815719 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
#> names : NDVI_median, NDVI_IQR, NDVI_slope, TCB_median, TCB_IQR, TCB_slope, TCW_median, TCW_IQR, TCW_slope, TCG_median, TCG_IQR, TCG_slope
#> min values : 1.261261e-01, 1.286144e-05, -2.372356e-01, 1.982568e+02, 5.820000e-02, -1.034494e+03, 4.720450e+01, 8.410000e-02, -3.770112e+02, -7.969700e+00, 5.355000e-02, -3.133347e+02
#> max values : 0.9029524, 0.3609146, 0.3609146, 4175.8897000, 1349.0888500, 1100.0643000, 1909.1103000, 635.7759000, 530.8367500, 2794.7694000, 634.9761000, 634.9761000
plot(VI_ts_metrics[["NDVI_median"]])