## Load pacman
library(pacman)
## Load rest of the packages
p_load(
giscoR, mapview, patchwork, sf, rsi, terra, tidyterra, tidyverse )
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:
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
<- gisco_get_nuts(
tenerife_sf country = "Spain",
resolution = "01",
nuts_level = 3
|>
) filter(
== "Tenerife"
NAME_LATN |>
) 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:
$planetary_computer_v1 sentinel2_band_mapping
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:
<- sentinel2_band_mapping$planetary_computer_v1[c("B02", "B03", "B04", "B08", "B11")]
sel_bands
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:
<- get_sentinel2_imagery(
tenerife_s2 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:
<- rast(tenerife_s2) / 10000 tenerife_sr
Before calculating the indices, we can check how does it look the RGB image:
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:
<- mask(tenerife_sr, tenerife_sf) tenerife_sr
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:
<- spectral_indices() |>
indices_tbl filter(
%in% c("NDMI", "NDVI")
short_name
)
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
<- calculate_indices(
indices_path raster = tenerife_sr,
indices = indices_tbl,
output_filename = glue::glue("{tempdir()}/spectral_indices.tif")
)
## Read indices
<- rast(indices_path)
indices_sr
## 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:
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
<- matrix(
ndmi_mat c(
-Inf, -.2, 1,
-.2, 0, 2,
0, .2, 3,
2, .4, 4,
.4, Inf, 5
.
),ncol = 3,
byrow = TRUE
)
## Classify NDMI
$NDMI <- indices_sr$NDMI |>
indices_srclassify(
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:
<- ggplot() +
ndmi_gg 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
<- ggplot() +
ndvi_gg 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
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"
)
) )
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.