1. Introduction

This notebook illustrates how to access time series of Sentinel-2 imagery using the packages rstac and datacubes.

rstac offers features that simplify the process of acquiring spatial data from STAC (SpatioTemporal Asset Catalog). More information here.

gdalcubes s an R package and C++ library aiming at making processing collections of satellite imagery easier, faster, and more interactive. More information here.

We will use the freely available Sentinel-2 COG catalog on Amazon Web Services (AWS) and the corresponding STAC-API endpoint at https://earth-search.aws.element84.com/v0/collections/sentinel-s2-l2a.

In addition, we will use terra for reading raster files and tmap for visualization.

2. Installation of packages

First of all, we need to install the rstac package.

Use the console to run the following code – one line at a time –:

## Do not run this code from here - Use the R console!!!
#install.packages("rstac")
#install.packages("tidyterra")
#install.packages("gdalcubes")
#install.packages("stars")
#install.packages("mapview")
#install.packages("tmap")

3. Setup

Clean the memory and next, load the libraries:

# clean workspace
rm(list = ls(all=TRUE))
# load libraries
library(rstac)
library(gdalcubes)
library(stars)
library(tmap)
library(tmaptools)
library(mapview)
library(leaflet)
library(sf)
library(ggplot2)

4. Define your area of interest (AOI)

Follow these instructions to create a GeoJSON file with your AOI:

In this notebook, the AOI is *La Unión in Valle del Cauca.

Then, read the geojson file:

#use the function list.files() to check your directory and filenames
union <- st_read("la_union/map.geojson")
## Reading layer `map' from data source 
##   `/Users/ials/Documents/unal2025/GB2024/data/la_union/map.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.14981 ymin: 4.437113 xmax: -75.9467 ymax: 4.601428
## Geodetic CRS:  WGS 84

5. Search Sentinel-2 imagery

The function stac_search() allows users to search Sentinel-2 data given an AOI with georeference EPSG:4326, a time period. Let’s try it:

s = stac("https://earth-search.aws.element84.com/v0")
items <- s |>
  stac_search(collections = "sentinel-s2-l2a-cogs",
              bbox = c(st_bbox(union)["xmin"],
                       st_bbox(union)["ymin"],
                       st_bbox(union)["xmax"],
                       st_bbox(union)["ymax"]), 
              datetime = "2025-01-01/2025-12-31",
              limit = 60) |>
  post_request() 
items
## ###Items
## - matched feature(s): 224
## - features (60 item(s) / 164 not fetched):
##   - S2B_18NUK_20251213_0_L2A
##   - S2C_18NUK_20251211_0_L2A
##   - S2A_18NUK_20251210_0_L2A
##   - S2C_18NUK_20251208_0_L2A
##   - S2C_18NUL_20251208_0_L2A
##   - S2B_18NUK_20251206_0_L2A
##   - S2B_18NUL_20251206_0_L2A
##   - S2B_18NUK_20251203_0_L2A
##   - S2B_18NUL_20251203_0_L2A
##   - S2C_18NUK_20251201_0_L2A
##   - ... with 50 more feature(s).
## - assets: 
## AOT, B01, B02, B03, B04, B05, B06, B07, B08, B09, B11, B12, B8A, info, metadata, overview, SCL, thumbnail, visual, WVP
## - item's fields: 
## assets, bbox, collection, geometry, id, links, properties, stac_extensions, stac_version, type

Let print date and time of first and last images:

range(sapply(items$features, function(x) {x$properties$datetime}))
## [1] "2025-10-17T15:41:46Z" "2025-12-13T15:32:03Z"

6. Create a Sentinel-2 collection

We can convert the STAC response to a gdalcubes image collections using stac_image_collection(). This function expects a STAC feature list as input and optionally can apply some filters on metadata and bands. Notice that this operation is much faster than the creation of image collections from local files, because all metadata are already available and no image file must be opened. Below, we create a collection for images with less than 50% cloud cover.

assets = c("B01","B02","B03","B04","B05","B06", "B07","B08","B8A","B09","B11","SCL")
s2_collection = stac_image_collection(items$features, asset_names = assets, property_filter = function(x) {x[["eo:cloud_cover"]] < 20})

Let’s check wa we got:

s2_collection
## Image collection object, referencing 3 images with 12 bands
## Images:
##                       name      left      top   bottom     right
## 1 S2C_18NUK_20251208_0_L2A -76.34752 4.523459 3.529607 -75.81212
## 2 S2C_18NUL_20251208_0_L2A -76.14920 5.427668 4.434019 -75.81300
## 3 S2C_18NUL_20251118_0_L2A -76.14245 5.427668 4.434029 -75.81300
##              datetime        srs
## 1 2025-12-08T15:32:13 EPSG:32618
## 2 2025-12-08T15:31:59 EPSG:32618
## 3 2025-11-18T15:32:01 EPSG:32618
## 
## Bands:
##    name offset scale unit nodata image_count
## 1   B01      0     1                       3
## 2   B02      0     1                       3
## 3   B03      0     1                       3
## 4   B04      0     1                       3
## 5   B05      0     1                       3
## 6   B06      0     1                       3
## 7   B07      0     1                       3
## 8   B08      0     1                       3
## 9   B09      0     1                       3
## 10  B11      0     1                       3
## 11  B8A      0     1                       3
## 12  SCL      0     1                       3

In this case, the collection contains 3 images.

7. Create and process data cubes

Having an image collection object, we can now use gdalcubes. Particularly, we define a data cube view and maybe some additional data cube operations. In the following example, we create a coarse resolution (100m) simple median-composite RGB image from a monthly data cube.

(union_box = st_bbox(union))
##       xmin       ymin       xmax       ymax 
## -76.149811   4.437113 -75.946702   4.601428
(union_b9377 <- st_bbox(st_transform(union, "EPSG:9377")))
##    xmin    ymin    xmax    ymax 
## 4650592 2048948 4673070 2067185
gdalcubes_options(parallel = 8)
(v = cube_view(srs="EPSG:9377", dx=10, dy=10, dt="P1M", 
                           aggregation="median", resampling = "average",
                           extent=list(t0 = "2025-01-01", t1 = "2025-12-31",
                                       left=union_b9377["xmin"], right=union_b9377["xmax"],
                                       top=union_b9377["ymax"], bottom=union_b9377["ymin"])))
## A data cube view object
## 
## Dimensions:
##                low             high count pixel_size
## t       2025-01-01       2025-12-31    12        P1M
## y 2048946.65931213 2067186.65931213  1824         10
## x 4650590.60828835 4673070.60828835  2248         10
## 
## SRS: "EPSG:9377"
## Temporal aggregation method: "median"
## Spatial resampling method: "average"

8. Visualizing the data cube

First, natural color (i.e. RGB 432):

raster_cube(s2_collection, v) |>
  select_bands(c("B02","B03","B04")) |>
  reduce_time(c("median(B02)", "median(B03)", "median(B04)")) |>
  plot(rgb = 3:1, zlim = c(0,2500))

Now, false color (i.e. RGB 8-11-4):

raster_cube(s2_collection, v) |>
  select_bands(c("B04","B11","B08")) |>
  reduce_time(c("median(B04)", "median(B11)", "median(B08)")) |>
  plot(rgb = 3:1, zlim = c(0,4500))

9. Download the data

# we "download" the data and save it as a netcdf  or a tif file
tutorialDir = path.expand("./la_union/")
if (!file.exists((file.path(tutorialDir,"Union_10"))))
{
  s2.mask = image_mask("SCL", values = c(3,8,9))
  gdalcubes_options(parallel = 24, ncdf_compression_level = 0)
  raster_cube(s2_collection, v, mask = s2.mask) |>
    write_tif(file.path(tutorialDir,"Union_10.tif"))
}

Check names of downloaded files:

list.files("la_union/Union_10")
##  [1] "cube_12e097724472025-01-01.tif" "cube_12e097724472025-02-01.tif"
##  [3] "cube_12e097724472025-03-01.tif" "cube_12e097724472025-04-01.tif"
##  [5] "cube_12e097724472025-05-01.tif" "cube_12e097724472025-06-01.tif"
##  [7] "cube_12e097724472025-07-01.tif" "cube_12e097724472025-08-01.tif"
##  [9] "cube_12e097724472025-09-01.tif" "cube_12e097724472025-10-01.tif"
## [11] "cube_12e097724472025-11-01.tif" "cube_12e097724472025-12-01.tif"

10. Compute a vegetation index and visualize it

First, read one GTiff file:

f <- "la_union/Union_10/cube_12e097724472025-12-01.tif"
union_s2 <- terra::rast(f)

Then, compute the NDVI index:

# 1. Define the bands
red <- union_s2[[4]] # Band 4
nir <- union_s2[[8]] # Band 8

# 2. Calculate NDVI
ndvi <- (nir - red) / (nir + red)

What we got?

ndvi
## class       : SpatRaster
## size        : 1824, 2248, 1  (nrow, ncol, nlyr)
## resolution  : 10, 10  (x, y)
## extent      : 4650591, 4673071, 2048947, 2067187  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s)   : memory
## varname     : cube_12e097724472025-12-01
## name        :      B08
## min value   : -0.99974
## max value   : 0.999544

Now, visualization time. Let’s start with a simple plot:

# Map NDVI with a continuous Red-Yellow-Green palette
tm_shape(ndvi) +
  tm_raster(
    col_palette = "RdYlGn", 
    style = "cont",
    title = "NDVI Value"
  ) +
  tm_layout(main.title = "Sentinel-2 NDVI Dec. 2025")

As you may know, the results of the NDVI calculation range from -1 to 1. Negative values correspond to areas with water surfaces, man-made structures, rocks, clouds, snow; bare soil usually falls within 0.1- 0.2 range; and plants will always have positive values between 0.2 and 1.

Healthy, dense vegetation canopy should be above 0.5, and sparse vegetation will most likely fall within 0.2 to 0.5. However, it’s only a rule of thumb and you should always take into account the season, type of plant and regional peculiarities to know exactly what NDVI values mean.

Another visualization may help to interpret NDVI values:

# 2. Define the requested Red -> Yellow -> Green color palette
# We wrap the domain values to match the theoretical bounds of NDVI (-1 to 1)
ndvi_palette <- colorNumeric(
  palette = c("red", "yellow", "green"), 
  domain = c(-1, 1),
  na.color = "transparent" # Keeps pixels with missing data invisible
)

# 3. Build the interactive Leaflet Map
map <- leaflet() %>%
  # Add the background Esri Satellite imagery layer
  addProviderTiles("Esri.WorldImagery", group = "Esri Satellite") %>%
  # Define the map view: longitude, latitude, and specific zoom level
  setView(lng = -76.0, lat = 4.50, zoom = 14) %>%
  # Add the custom NDVI raster layer overlay
  addRasterImage(
    x = ndvi, 
    colors = ndvi_palette, 
    opacity = 0.6,          # Adjust opacity so the underlying satellite imagery is partially visible
    group = "NDVI Layer"    # Group identifier for the layer control panel
  ) %>%
  
  # Add a static visual color legend to the map
  addLegend(
    pal = ndvi_palette, 
    values = c(-1, 1), 
    title = "NDVI Index",
    position = "bottomright",
    opacity = 0.8
  ) %>%
  
  # Add structural toggles so you can turn the NDVI map layer on or off
  addLayersControl(
    baseGroups = c("Esri Satellite"),
    overlayGroups = c("NDVI Layer"),
    options = layersControlOptions(collapsed = FALSE)
  )

# Display the map
map

Please note that the date of the satellite image shown in the background does not necessarily match the date of the NDVI.

All you need to know about NDVI may be here

Note: Replicate this notebook for your AOI
Make sure you add your our comments and NDVI interpretation for two different dates

That’s all for today.