NDVI and NDMI in Tenerife before huge wildire

Download satellite images from Sentinel-2, calculate the NDVI and NDMI indices, and generate a couple of maps.
remote sensing
ggplot2
R
spatial
Author

Adrián Cidre

Published

June 30, 2024

1 Introduction

In this post, we will create a map of the Normalized Difference Vegetation Index (NDVI) and the Normalized Difference Moisture Index (NDMI) for the Tenerife Island, using the average between the 2023-07-15 and 2023-08-15, one month before the wildfire that burned around 15,000 hectares.

We will use the awesome rsi package (rsi?) to access the Sentinel-2 imagery and calculate the indices. This package is used to retrieve data from different STAC servers that provides satellite imagery. It is also possible to use this package to easily get data from Sentinel-1 and Landsat.

Therefore, in this exercise, we will see how was the value of the indices in average one month before the wildfire of Tenerife.

2 Loading packages

We will use the following packages:

## Load pacman
library(pacman)

## Load rest of the packages
p_load(
  giscoR, mapview, patchwork, sf, rsi, terra, tidyterra, tidyverse
)

We will use the giscoR package to download the boundaries of the Tenerife Island, the rsi package to retrieve the Sentinel-2 data and to calculate the indices, the sf package to manipulate vectorial data, the terra package to manipulate raster data, mapview for interactive visualizations, patchwork for plots assembling, and tidyterra for easy visualizations of SpatRasters.

3 Load the data

In this exercise we need to download two sources of data:

  • Study area

  • Satellite image

3.1 Study area

The study area in this case is the Tenerife Island. We can easily load it into R using the giscoR package:

## Get study area
tenerife_sf <- gisco_get_nuts(
  country    = "Spain",
  resolution = "01",
  nuts_level = 3
) |> 
  filter(
    NAME_LATN == "Tenerife"
  ) |> 
  st_transform(25828)

First, we download the Spanish provinces, then we filter only Tenerife, and finally we transform the data to a different coordinates reference system.

mapview(tenerife_sf)

The CRS that you use must be a projected system in meters, since this is required by the rsi package.

3.2 Satellite image

We want to download an image from the Sentinel-2 imagery for the study area that we defined previously. Each Sentinel-2 image contains a total of 12 different bands (variables). However, for our analysis we don’t need all of them. The rsi package makes it easy for us to just download the bands that we need for our analysis. The next dataset contains information about the bands, and it can be used to filter the desired bands:

sentinel2_band_mapping$planetary_computer_v1
An rsi band mapping object with attributes:
names mask_band mask_function stac_source collection_name query_function class scl_name download_function sign_function

  B01   B02   B03   B04   B05   B06   B07   B08   B8A   B09   B11   B12 
  "A"   "B"   "G"   "R" "RE1" "RE2" "RE3"   "N"  "N2"  "WV"  "S1"  "S2"

This is essentially a character vector, so we will filter the RGB bands plus the NIR and SWIR1 that we need to calculate the indices:

sel_bands <- sentinel2_band_mapping$planetary_computer_v1[c("B02", "B03", "B04", "B08", "B11")]

sel_bands
An rsi band mapping object with attributes:
mask_band mask_function stac_source collection_name query_function class scl_name download_function sign_function names

 B02  B03  B04  B08  B11 
 "B"  "G"  "R"  "N" "S1"

Now that we have the bands that we need, we can use the rsi::get_sentinel2_imagery() function to retrieve the satellite image. Here, we will specify the arguments:

  • aoi: area of interest

  • start_date: initial date for retrieving images

  • end_date: end data for retrieving images

  • asset_names: the rsi band mapping object with the desired bands

  • output_filename: output where the image will be downloaded. In this case, I just save it to the temporary directory.

Note that there is one argument called composite_function = "median", which mean that the images between start and end date will be reduced to the median. So now, it’s time to get out image:

tenerife_s2 <- get_sentinel2_imagery(
  aoi             = tenerife_sf,
  start_date      = "2023-07-15",
  end_date        = "2023-08-15",
  asset_names     = sel_bands,
  output_filename = glue::glue("{tempdir()}/sentinel.tif")
)

The images are stored in the cloud after multiplying them by a factor of 10,000 to avoid the storage of floats. Therefore, we will apply back this factor to get the real digital values:

tenerife_sr <- rast(tenerife_s2) / 10000

Before calculating the indices, we can check how does it look the RGB image:

plotRGB(
  tenerife_sr,
  3, 2, 1,
  scale   = 2,
  stretch = "lin"
)
Figure 1: RGB image of Tenerife, averaged from 15th July 2023 to 15th July 2023. Source: Sentinel-2

4 Data preparation

Before proceeding with the analysis, I just want to get the indices for the Island, not for the sea. I mean, as you can see in Figure 1, the image is a rectangle where we also have the sea. To eliminate these pixels, we can mask those values as follows:

tenerife_sr <- mask(tenerife_sr, tenerife_sf)

We are basically converting to NA the values of tenerife_sr which are outside of tenerife_sf. Excellent!! Let’s go to calculate the indices.

4.1 Select indices

We could manually calculate the indices, since it’s very simple. However, I want to show you something from the rsi package that really stands out! We have access to the awesome spectral indices database through the {spectral_indices} function.

spectral_indices()
# A tibble: 243 × 9
   application_domain bands     contributor   date_of_addition formula long_name
   <chr>              <list>    <chr>         <chr>            <chr>   <chr>    
 1 vegetation         <chr [2]> https://gith… 2021-11-17       (N - 0… Aerosol …
 2 vegetation         <chr [2]> https://gith… 2021-11-17       (N - 0… Aerosol …
 3 water              <chr [6]> https://gith… 2022-09-22       (B + G… Augmente…
 4 vegetation         <chr [2]> https://gith… 2021-09-20       (1 / G… Anthocya…
 5 vegetation         <chr [3]> https://gith… 2022-04-08       N * ((… Anthocya…
 6 vegetation         <chr [4]> https://gith… 2021-05-11       (N - (… Atmosphe…
 7 vegetation         <chr [4]> https://gith… 2021-05-14       sla * … Adjusted…
 8 vegetation         <chr [2]> https://gith… 2022-04-08       (N * (… Advanced…
 9 water              <chr [4]> https://gith… 2021-09-18       4.0 * … Automate…
10 water              <chr [5]> https://gith… 2021-09-18       B + 2.… Automate…
# ℹ 233 more rows
# ℹ 3 more variables: platforms <list>, reference <chr>, short_name <chr>

A total of 243 indices that can be calculated 😱. Honestly, I don’t know most of them and I think for a tutorial like this one is enough to showcase two of them. So we can filter them from that tibble using a regular filter:

indices_tbl <- spectral_indices() |> 
  filter(
    short_name %in% c("NDMI", "NDVI")
  )

indices_tbl
# A tibble: 2 × 9
  application_domain bands     contributor    date_of_addition formula long_name
  <chr>              <list>    <chr>          <chr>            <chr>   <chr>    
1 vegetation         <chr [2]> https://githu… 2021-12-01       (N - S… Normaliz…
2 vegetation         <chr [2]> https://githu… 2021-04-07       (N - R… Normaliz…
# ℹ 3 more variables: platforms <list>, reference <chr>, short_name <chr>

4.2 Calculate indices

Once we filtered the indices we want in the spectral_indices() table, we can use the calculate_indices() function to obtain them. And this is very straightforward, you just need the image, the filtered indices in the table, and specify the output file name.

## Calculate indices
indices_path <- calculate_indices(
  raster          = tenerife_sr,
  indices         = indices_tbl,
  output_filename = glue::glue("{tempdir()}/spectral_indices.tif")
)

## Read indices
indices_sr <- rast(indices_path)

## Print
indices_sr
class       : SpatRaster 
dimensions  : 6556, 7895, 2  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 311453.7, 390403.7, 3098149, 3163709  (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 28N (EPSG:25828) 
source      : 007_indices.tif 
names       :       NDMI,       NDVI 
min values  : -0.5467158, -0.5689685 
max values  :  0.6188195,  0.7418026

As you can see, we have a SpatRaster with two bands, one per index.

4.3 Classify NDMI

Instead of showing the NDMI as a continuous variable, we can classify it using some thresholds. In this case, I followed the recommendations given my EOS Data Analytics slightly modified. We can have a look at the histograms of these two variables:

hist(indices_sr)

There are almost no values in the tails, so I decided to reduce the classification for this exercise to the following 5 classes:

## Classification matrix for NDMI
ndmi_mat <- matrix(
  c(
    -Inf, -.2, 1,
    -.2, 0, 2,
    0, .2, 3,
    .2, .4, 4,
    .4, Inf, 5
  ),
  ncol  = 3,
  byrow = TRUE
)

## Classify NDMI
indices_sr$NDMI <- indices_sr$NDMI |> 
  classify(
    rcl = ndmi_mat
  ) |> 
  as.factor()

## Label values
levels(indices_sr$NDMI)[[1]][, 2] <- c(
  "Low or very low canopy cover, dry or very low canopy cover, wet",
  "Mid-low canopy cover, high water stress or low canopy cover, low water stress",
  "Average canopy cover, high water stress or mid-low canopy cover, low water stress",
  "Mid-high canopy cover, high water stress or average canopy cover, low water stress",
  "No water stress"
)

5 Visualize

The last step in this analysis is to create a map. We create two maps independently, and then assemble them together using the patchwork package:

ndmi_gg <- ggplot() +
  geom_spatraster(
    data = indices_sr$NDMI
  ) +
  scale_fill_whitebox_d(
    palette   = "bl_yl_rd",
    direction = -1
  ) +
  guides(
    fill = guide_legend(
      position = "inside",
      title    = NULL
    )
  ) +
  theme_void(
    base_size   = 8,
    base_family = "Roboto"
  ) +
  theme(
    legend.position.inside = c(.3, .8),
    legend.key.spacing.y   = unit(2, "mm"),
    legend.key.width       = unit(5, "mm"),
    legend.key.height      = unit(1, "mm"),
    legend.key             = element_rect(colour = "black", linewidth = .2),
    legend.text            = element_text(size = 5),
  )

ndmi_gg
Figure 2: Normalized Difference Moisture Index reduced into 5 classes 
ndvi_gg <- ggplot() +
  geom_spatraster(
    data = indices_sr$NDVI
  ) +
  scale_fill_whitebox_c(
    palette   = "muted",
    direction = -1
  ) +
  guides(
    fill = guide_colorbar(
      position       = "inside",
      title          = "NDVI",
      title.position = "top",
      title.hjust    = .5,
      direction      = "horizontal"
    )
  ) +
  theme_void(
    base_size   = 8,
    base_family = "Roboto"
  ) +
  theme(
    legend.position.inside = c(.3, .8),
    legend.key.width       = unit(1, "cm"),
    legend.key.height      = unit(1.5, "mm"),
    legend.text            = element_text(size = 5),
    legend.title           = element_text(size = 8)
  )

ndvi_gg
Figure 3: Normalized Difference Vegetation Index

So the final visualization is created with patchwork as follows:

ndvi_gg + 
  ndmi_gg +
    plot_annotation(
      title    = "Pre-Wildfire Vegetation and Moisture Analysis of Tenerife",
      subtitle = "NDVI and NDMI averaged from 15-July 2023 to 15-August 2023",
      caption  = "Author: Adrián Cidre | https://adrian-cidre.com | Data source: Sentinel-2",
      theme    = theme(
        plot.title = element_text(
          family     = "Merriweather",
          face       = "bold",
          lineheight = 1.2,
          margin     = margin(t = 5, l = 5, b = 10),
          hjust = .5
        ),
        plot.subtitle = element_text(hjust = .5),
        plot.caption  = element_text(
          hjust  = .5,
          family = "Roboto"
        )
      )
    )
Figure 4: Final visualization with NDMI and NDVI assembled with the {patchwork} package

6 Conclusions

Using the rsi package is very easy to retrieve satellite images from different sources. In this exercise, we saw how to obtain a Sentinel-2 image for a given area, and we calculated two spectral indices using the same package.

If you found this post interesting and helpful, please share it with others or create your own, and tag me in linkedin. You can also watch the video in YouTube to learn more.

References