rterra

How to calculate NDVI from a four-dimensional raster?


I have several Sentinel-2 satellite scenes for the same location from different dates. Each scene contains basic spectral bands (including red "B04" and near infrared "B08"). I loaded the data as a SpatRasterDataset using the terra package. How I can calculate the NDVI defined as NDVI = (NIR - RED) / (NIR + RED) for each date?

library(terra)
r <- rast(system.file("ex/logo.tif", package="terra")) / c(1,2,3)
r <- c(r, r[[1]] + r[[2]])
names(r) <- paste0("B0", c(2:4, 8)) 
dates <- c("2023-02-08", "2023-03-02", "2023-04-06")
s <- sds(r, r, r)
names(s) <- dates

Solution

  • Example data

    library(terra)
    r <- rast(system.file("ex/logo.tif", package="terra")) / c(1,2,3)
    r <- c(r, r[[1]] + r[[2]])
    names(r) <- paste0("B0", c(2:4, 8)) 
    dates <- c("2023-02-08", "2023-03-02", "2023-04-06")
    s <- sds(r, r, r)
    names(s) <- dates
    

    Solution (here assuming "BO2" is "NIR" and "B04" is "RED")

    ndvi <- function(NIR, RED) (NIR - RED) / (NIR + RED)
    v <- lapply(s, \(x) lapp(x[[c("B02", "B04")]], ndvi))
    v <- rast(v)
    panel(v, nr=1)
    

    enter image description here

    You could use a list instead of a SpatRasterDataset

    ss <- list(r, r, r)
    names(ss) <- dates
    vv <- lapply(ss, \(x) lapp(x[[c("B02", "B04")]], ndvi)) |> rast()