Visualizing vegetation phenological trend.

Here we use ridgeplots to display the variation of vegetation parameter (NDVI) related to phenology on an annual scale making use of the ‘ggridges’ R package.

We start by calling the libraries and importing the files:

## Libraries ####

library(raster)
## Caricamento del pacchetto richiesto: sp
library(terra)
## terra 1.7.55
library(purrr)
library(furrr)
## Caricamento del pacchetto richiesto: future
library(ggridges)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks terra::extract(), raster::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ dplyr::select()  masks raster::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Import data

setwd("C:/Tirocinio")

NDVI_300m <- list.files(pattern = "NDVI300_2014", recursive = T) %>%
  map(~rast(.))

# select one raster per each month to simplify calculations

NDVI_300m <- NDVI_300m[seq(1, length(NDVI_300m), by = 3)]

Preprocessing

Here we crop the rasters to visualize the NDVI in Italy

## Pre-processing ####

# crop the rasters on Italian extent

Italy <- getData("GADM", country = "ITA", level = 0)
## Warning in getData("GADM", country = "ITA", level = 0): getData will be removed in a future version of raster
## . Please use the geodata package instead
NDVI_300m_crop <- lapply(NDVI_300m, crop, Italy)

# plot the results

plot(NDVI_300m_crop[[1]])

Now, let’s extract the data from the raster files and let’s tidy it for the final ridgeplot!

# Extract NDVI values in a list

values_NDVI_300m <- list()

for (i in 1:12) {
  
  new_object <- terra::as.data.frame(NDVI_300m_crop[[i]])
  
  values_NDVI_300m[[paste0("NDVI_300m_", i)]] <- new_object
  
  values_NDVI_300m[[i]]$month <- as.numeric(paste(i))  # create a column reporting the month of reference for each NDVI value
  
  rm(new_object)
  
}


# Merge list objects into a single dataframe

values_NDVI300_tot <- do.call(rbind, values_NDVI_300m) %>%
  
  subset(NDVI >= 0.2 | is.na(NDVI)) # set the condition excluding uninteresting cover types and NAs


# order rows according to month

values_NDVI300_tot <- values_NDVI300_tot[order(values_NDVI300_tot$month), ]


# check the varaible levels

unique(values_NDVI300_tot$month)
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12
# turn variable month into a factor for the ridgeplot function requirements

values_NDVI300_tot$month <- as.factor(values_NDVI300_tot$month)


# Correct months names

values_NDVI300_tot$month <- fct_recode(values_NDVI300_tot$month, 
                                       January = "1", 
                                       February = "2",
                                       March = "3",
                                       April = "4",
                                       May = "5",
                                       June = "6",
                                       July = "7",
                                       August = "8",
                                       September ="9",
                                       October = "10",
                                       November = "11",
                                       December = "12")

Ridgeplots

Now that the dataset is ready, let’s go plotting our ridgeplots!

## Ridgeplots ####

# free some memory

rm(NDVI_300m, NDVI_300m_crop, values_NDVI_300m)

gc()
##             used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells  74921351 4001.3  121862690  6508.2   89064258  4756.6
## Vcells 539268755 4114.3 1585959634 12100.0 1561181829 11910.9
# Simple NDVI ridgeplot with quantiles

ridgeplot_NDVI <- ggplot(values_NDVI300_tot, aes(NDVI, month, fill = factor(stat(quantile)))) +
  stat_density_ridges(
    geom = "density_ridges_gradient", calc_ecdf = TRUE,
    quantiles = 4, quantile_lines = TRUE)

print(ridgeplot_NDVI)
## Warning: `stat(quantile)` was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(quantile)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Picking joint bandwidth of 0.00709
## Warning: Using the `size` aesthetic with geom_segment was deprecated in ggplot2 3.4.0.
## ℹ Please use the `linewidth` aesthetic instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Gradient ridgeplot type

ridgeplot_NDVI <- ggplot(values_NDVI300_tot, aes(x = `NDVI`, y = month, fill = after_stat(x))) +
  geom_density_ridges_gradient(scale = 3, rel_min_height = 0.01, gradient_lwd = 1.) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_discrete(expand = expansion(mult = c(0.01, 0.25))) +
  scale_fill_viridis_c(name = "NDVI", option = "C") +
  labs(x = "NDVI",
       y = "Month",
    title = 'Phenological trend in Italy (2017)',
    subtitle = 'NDVI by month for 2017'
  ) +
  theme_ridges(font_size = 13, grid = TRUE)

print(ridgeplot_NDVI)
## Picking joint bandwidth of 0.00709