Loading and analysing Digital Earth Australia Sentinel-2 data in R with `rstac` and `gdalcubes`

dea_in_r.R

"""
This code demonstrates how to load Digital Earth Australia Sentinel-2 Analysis Ready Data into R.
It uses `rstac` to search for available data for a time and location using DEA's STAC endpoint,
and `gdalcubes` to load and analyse the data.

Functionality includes:
    * Creating a custom pixel grid to reproject data into
    * Apply a cloud mask using the "s2cloudless" cloud mask
    * Combine data into seasonal composites
    * Create an RGB animation
    * Convert data to NDVI and plot across time as a timeseries plot and animation

Inspired by content in the following tutorial:
https://gdalcubes.github.io/source/tutorials/vignettes/gc02_AWS_Sentinel2.html
"""

library(magrittr)
library(colorspace)
library(rstac)
library(gdalcubes)

# Set parallelisation and AWS credentials to allow anonymous access
Sys.setenv("AWS_NO_SIGN_REQUEST" = "YES")
gdalcubes_options(parallel = 16)

# Set params
time = "2023-01-01/2023-12-31"
bbox = list(left=149.335, top=-35.31, right=149.375, bottom=-35.35)

# Connect to DEA STAC endpoint
items <- stac("https://explorer.dea.ga.gov.au/stac") %>%
  
  # Set up search query for Sentinel-2 data in time period and bounding box
  stac_search(collections = c("ga_s2am_ard_3", "ga_s2bm_ard_3"),
              datetime = time,
              bbox = bbox,
              limit = 1000) %>%

  # Run search query
  get_request()

# Create gdalcube image collection, filtering to specific bands and cloud cover
s2_collection <- stac_image_collection(items$features,
                                       asset_names = c("nbart_red", "nbart_green", "nbart_blue", "nbart_nir_1", "oa_s2cloudless_mask"),
                                       property_filter = function(x) {x[["eo:cloud_cover"]] < 50})

# Creat a custom pixel grid and temporal aggregation
cube_grid <- cube_view(srs = "EPSG:4326", 
                         extent = list(left=bbox$left, 
                                       right=bbox$right, 
                                       top=bbox$top, 
                                       bottom=bbox$bottom,
                                       t0="2023-01-01",
                                       t1="2023-12-31"), 
                         dx = 0.0001,
                         dy = 0.0001, 
                         dt = "P3M",  # create seasonal composites
                         resampling = "nearest", 
                         aggregation = "median")

# Create a cloud mask from the s2cloudless cloud mask layer, treating nodata and cloud as bad
s2_mask = image_mask("oa_s2cloudless_mask", values = c(0, 2))

# Resample satellite data into grid after applying cloud mask
s2_cube <- raster_cube(image_collection = s2_collection, 
                       view = cube_grid, 
                       mask = s2_mask,
                       chunking = c(1, 2024, 2024))

# Plot results in RGB
plot(s2_cube, rgb = c(4, 2, 1), zlim = c(0, 2000))
animate(s2_cube, rgb = c(4, 2, 1), zlim = c(0, 2000), fps = 4)


# Example NDVI analysis
ndvi <- s2_cube %>%
  apply_pixel("(nbart_nir_1-nbart_green)/(nbart_nir_1+nbart_green)", "NDVI")

ndvi %>%
  reduce_space("min(NDVI)", "max(NDVI)", "mean(NDVI)") %>%
  plot(join.timeseries = TRUE)

ndvi_col = function(n) { rev(sequential_hcl(n, "Green-Yellow")) }
ndvi %>%
  animate(key.pos = 1, zlim = c(-0.2,1), col = ndvi_col, nbreaks = 100, fps = 4)