1. Introduction

This notebook illustrates how to access Sentinel-2 imagery using R. It is based on this vignette from Féret & de Boissieu (2026).

We will use preprocS2, a package dedicated to Sentinel-2 data handling, which provides access to Spatio Temporal Asset Catalogs (STAC) via R scripts. It includes wrapper functions for the rstac package. It allows downloading and preparing Sentinel-2 images for further processings using biodivMapR.

The STAC specification is a common language to describe geospatial information. STAC resources are available online from various institutional and private providers (such as Copernicus dataspace and Nicrosoft Planetary Computers). STAC catalogs provides access to imagery data via spatio-temporal queries. Hence, the first step is to define an area of interest (aoi), and a period or day of interest.

The default STAC collection is sentinel-2-l2a from the Microsoft Planetary Computer catalog. Before proceeding, it is necessary to understand some things about Sentinel-2 imagery. Make sure to take notes about Sentinel-2 spatial and spectral resolution – on your paper notebook.

We will also use the spinR package that computes spectral indices from raster data.

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

2. Installation of packages

First of all, we need to install several packages:

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("remotes")
#remotes::install_github('jbferet/preprocS2')
#remotes::install_github('cran/dissUtils')
#install.packages("spam", repos = "https://reinhardfurrer.r-universe.dev")
#remotes::install_github('jbferet/biodivMapR')
#install.packages("pak")
#pak::pak("gitlab::jbferet/spinR")
#install.packages("tidyterra")

3. Setup

Clean the memory and next, load the libraries:

# clean workspace
rm(list = ls(all=TRUE))
# load libraries
library(preprocS2)  # install instructions: https://github.com/jbferet/preprocS2
library(biodivMapR) # install instructions: https://github.com/jbferet/biodivMapR
library(spinR)      # install instructions: https://github.com/jbferet/spinR
library(terra)
library(tidyterra)
library(sf)
library(ggplot2)

4. Define your area of interest (AOI)

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

Then, read the geojson file:

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

Now, let’s define input & output directories and files:

input_dir_vect <- './lebrija'
input_dir_rast <- './lebrija/sentinel-2'
aoi_path <- file.path(input_dir_vect, 'map.geojson')

5. Define a data of acquisition

A Sentinel-2 acquisition with minimal cloud cover needs to be identified. The Copernicus Browser is a useful service to explore the availability of Sentinel-2 acquisitions. Search images using your municipality’s name. In addition to spatio-temporal filters, it includes filters such as the maximum cloud cover.

We searched images acquired from January 2026 onwards. An acquisition from 17 January 2026 shows decent conditions.

6. Download Sentinel-2 imagery

The library preprocS2 allows downloading from multiple providers. Microsoft Planetary computer is currently the default option.

The function get_s2_raster downloads Sentinel-2 data corresponding to the spatio-temporal query with these parameters:

Let’s try it:

# Change values as needed: 
datetime <- as.Date('2026-01-17') 
stac_info <- list('provider' = 'mpc')
site_name = 'lebrija'
options <- set_options_preprocS2(fun = 'get_s2_raster')
options$overwrite <- FALSE
list_files <- get_s2_raster(aoi_path = aoi_path, 
                            datetime = datetime, 
                            stac_info = stac_info, 
                            output_dir = input_dir_rast, 
                            site_name = site_name, 
                            options = options)
## get S2 tiles corresponding to aoi
## download S2 collection
rast_path <- list_files$Refl_L2A                # S2 L2A reflectance
mask_path <- list_files$vegetation_mask         # S2 binary mask identifying vegetation, discarding clouds & shadows

7. Check what Sentinel-2 files were obtained

After several minutes, the requested data should be in your output folder, organized on a similar structure as this one:

The collections/plot_001.rds corresponds to the item collection resulting from the spatiotemporal query. Use the R function readRDS to access the content of the file.

readRDS('./lebrija/sentinel-2/collections/plot_001.rds')
## ###Items
## - features (1 item(s)):
##   - S2C_MSIL2A_20260117T152651_R025_T18NXN_20260117T183212
## - assets: 
## AOT, B01, B02, B03, B04, B05, B06, B07, B08, B09, B11, B12, B8A, datastrip-metadata, granule-metadata, inspire-metadata, product-metadata, rendered_preview, safe-manifest, SCL, tilejson, visual, WVP
## - item's fields: 
## assets, bbox, collection, geometry, id, links, properties, stac_extensions, stac_version, type

The raster_samples directory includes:

─ lebrija_001_2024-01-17.tiff: L2A (surface reflectance) Sentinel-2 data

─ lebrija_001_2026-01-17_BIN.tiff: binary mask derived from lebrija_001_2026-01-17_SCL.tiff and focusing on vegetation class.

─ lebrija_001_2026-01-17_BIN_v2.tiff: binary mask produced from alternative radiometric criteria, including NDVI mask to remove non vegetated pixels, Blue mask to remove hazy/cloudy pixels, and NIR mask to remove shaded pixels. You can adjust threshold for cloudMask (B02), shadeMask (B08), and NDVIMask. See here for additional information on how these masks work.

─ lebrija_001_2024-01-17_SCL.tiff: scene classification provided by sen2cor atmospheric correction.

8. Visualization time

First, let’s read the Sentinel-2 imagery and know its metadata:

# 1. Load the Sentinel-2 bands (e.g., stacked into a multi-layer .tif)
# B4 = Red, B3 = Green, B2 = Blue for Sentinel-2
(s2_image <- rast("./lebrija/sentinel-2/raster_samples/lebrija_001_2026-01-17.tiff"))
## class       : SpatRaster
## size        : 1343, 1826, 10  (nrow, ncol, nlyr)
## resolution  : 10, 10  (x, y)
## extent      : 680190, 698450, 775440, 788870  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source      : lebrija_001_2026-01-17.tiff
## names       :   B02,   B03,   B04,   B05,   B06,   B07, ...
## min values  :    12,    38,  -128,    22,   -81,   -82, ...
## max values  : 14200, 14056, 14104, 14553, 14527, 14594, ...
# Calculate the mean for each layer
(layer_means <- global(s2_image, fun = c("mean", "min", "max"), na.rm = TRUE))

Note two things: (i) there are ten spectral bands; and (ii) the range of values per band is as higher as 15000.

# Apply a custom 1% linear stretch (ignores bottom 1% and top 1% of outliers)
s2_enhanced <- terra::stretch(s2_image, 
                        minq = 0.01, 
                        maxq = 0.99, 
                        minv = 100, 
                        maxv = 12500)

First, a true color composite:

# 2. True color RGB 432
ggplot() +
  geom_spatraster_rgb(data = s2_enhanced, r = 3, g = 2, b = 1, max_col_value= 4800, stretch = "hist") +
  theme_minimal() +
  labs(title = "Sentinel-2 True Color Image")
## <SpatRaster> resampled to 500775 cells.

In your notebook, describe what land cover classes you see.

Now, a false color composite:

# 2. False color RGB B8 – B11 – B3
ggplot() +
  geom_spatraster_rgb(data = s2_enhanced, r = 7, g = 9, b = 2, max_col_value= 6850, stretch = "hist") +
  theme_minimal() +
  labs(title = "Sentinel-2 Agricultural Composite")
## <SpatRaster> resampled to 500775 cells.

In your notebook, describe what land cover classes you see.

9. What have we learned from this notebook

Write down three things you have learned today in relation to:

10. References

Féret, J.-B., de Boissieu, F., 2019. biodivMapR: an R package for α‐ and β‐diversity mapping using remotely‐sensed images. Methods Ecol. Evol. 00:1-7. https://doi.org/10.1111/2041-210X.13310

Féret, J.-B., Asner, G.P., 2014. Mapping tropical forest canopy diversity using high-fidelity imaging spectroscopy. Ecol. Appl. 24, 1289–1296. https://doi.org/10.1890/13-1824.1