Case 01 – Brazil (🇧🇷)

Exploring environmental indicators in Brazil

Context

Between 2014 and 2021, the Brazilian Legal Amazon underwent intense environmental and anthropogenic transformations. On one hand, forest loss has accelerated due to logging, agriculture, and infrastructure expansion. At the same time, night-time lights —a proxy for urban growth and human activity— have increased in frontier municipalities, signaling intensified land use. Concurrently, changes in evapotranspiration (SEBAL-MODIS) patterns suggest modifications in local moisture dynamics and energy balance, reflecting the cumulative impact of deforestation and land conversion.

1 Objective

Visualize indicators such as forest loss, human activity and evapotranspiration over time

---

2 Functions Used

Variable	Function	Description
Forest loss	l4h_forest_loss()	Annual forest cover loss
Night-time lights	l4h_night_lights()	Proxy of urbanization and anthropogenic pressure
Evapotranspiration	l4h_sebal_modis()	Surface evapotranspiration from MODIS SEBAL

---

3 Workflow

3.1 Libraries

We need to verify and load all the necessary packages required for this case. It ensures that dependencies for spatial, environmental, and epidemiological analyses are properly installed and available. If any package is missing, install it manually with: install.packages()

for(p in c("geodata","sf","terra","tidyterra","rmapshaper","mapview","tmap",
  "tidyverse","lubridate","purrr","ggrepel",
  "ggplot2","ggpubr","RColorBrewer"
)){
  if(!require(p, character.only = TRUE)){
    install.packages(p)
  }
  library(p, character.only = TRUE)
}

library(land4health)

# ee_Authenticate()   # Run once if not authenticated
ee_Initialize()

── rgee 1.1.7 ─────────────────────────────────────── earthengine-api 0.1.317 ── 
 ✔ user: not_defined 
 ✔ Initializing Google Earth Engine:
 ✔ Initializing Google Earth Engine:  DONE!

 ✔ Earth Engine account: users/land4health_2 

 ✔ Python Path: C:/Users/yomal/AppData/Local/Programs/Python/Python39/python.exe 
──────────────────────────────────────────────────────────────────────────────── 

3.2 Spatial data preparation

We first load and define the administrative boundaries that will serve as the spatial reference for all analyses. The shapefiles are retrieved from the GADM dataset via the geodata package. All departments from Brazil are selected as the study area.

brasil_sf <- geodata::gadm(country = "BRA", level = 1, path = tempdir()) |> 
  st_as_sf()|> filter(NAME_1 %in%  c("Pará","Amapá","Amazonas",
                "Mato Grosso", "Rondônia", "Roraima", "Acre", "Tocantins", "Maranhão"))

brasil_sf <- rmapshaper::ms_simplify(brasil_sf, keep = 0.05)  

rmapshaper::ms_simplify() function simplifies the spatial geometry by removing redundant vertices while maintaining the overall shape and topology of the polygons.

3.3 Forest Loss

We need to use l4h_forest_loss() function which retrieves annual estimates of vegetation loss derived from the Global Forest Change (GFC) dataset developed by Hansen et al. (2013). In this case, the function extracts data for the period 2014–2021 within the administrative boundaries of Brazil (brasil_sf). Setting sf = TRUE returns the results as a spatial object (sf).

# Loss of vegetation cover
forest_loss <-  l4h_forest_loss(
  from = "2014-01-01",
  to = "2021-12-31",
  region = brasil_sf,
  sf = TRUE)%>%
  mutate(year = year(date))

Code

forest_loss%>%
  st_transform(32723) %>%
  mutate(
    X = st_coordinates(st_centroid(geometry))[,1],
    Y = st_coordinates(st_centroid(geometry))[,2]
  ) %>%
  ggplot() +
  geom_sf(aes(fill = value), color = NA) +
  geom_label(
    data = . %>% filter(value >= 8000),
    aes(x = X, y = Y, label = NAME_1, fill = value),
    color = "#333333",        
    size = 2.5,
    fontface = "bold",
    label.size = 0,         
    alpha = 0.8             
  ) +
  scale_fill_gradientn(
    colours = c("#00441B", "#1B7837", "#A6DBA0", "#FEE08B", "#E08214", "#B35806"),
    name = "Forest loss (km²)"
  ) +
  facet_wrap(~year, ncol = 4, nrow = 2) +
  theme_minimal(base_size = 12) +
  theme(
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    strip.text = element_text(size = 11, face = "bold", color = "#333333"),
    axis.text = element_blank(),
    axis.title = element_blank(),
    legend.position = "bottom",
    legend.key.width = unit(1, "cm"),
    legend.title = element_text(size = 10, face = "bold"),
    legend.text = element_text(size = 9)
  )

During the 2014–2021 period, Pará consistently exhibits the largest forest loss, with peaks observed in 2016–2017, followed by a gradual decrease in subsequent years. Mato Grosso and Maranhão also show significant deforestation during these peak years. From 2018 onward, a relative reduction in forest loss is observed across most regions.

3.4 Night-time Lights

Night-time light (NTL) intensity serves as a proxy for human activity, urbanization, and energy use, providing insight into the dynamics of development and land transformation. Here we use l4h_night_lights() function, which retrieves annual composites of satellite-derived radiance data.

night_lights<-l4h_night_lights(
  from = "2014-01-01",
  to = "2021-12-31",
  region = brasil_sf,
  stat = "mean"
)%>%
  mutate(year = year(date))

Code

tm_shape(night_lights) +
  tm_polygons(
    col = "value",
    palette = "viridis",
    title = "Radiance\n(nW/cm²/sr)"
  ) +
  tm_facets(by = "year", ncol = 4, nrow = 2) 

States in the southeastern and southern Amazon exhibit a gradual increase in brightness after 2019, indicating intensified human activity and urban expansion, whereas northern areas such as Amazonas and Roraima remain consistently low, reflecting limited electrification and preserved forest cover.

3.5 Evapotranspiration

Evapotranspiration quantifies the total water flux from land surfaces to the atmosphere, integrating soil evaporation and vegetation transpiration.
To estimate this variable, the function l4h_sebal_modis() will be use to derive mean annual evapotranspiration from MODIS-based SEBAL (Surface Energy Balance Algorithm for Land) data between 2014 and 2021.

sebal_annual <- l4h_sebal_modis(
  from   = "2014-01-01",
  to     = "2021-12-31",
  by     = "annual",
  fun    = "mean",
  region = brasil_sf,
  sf     = TRUE
)%>%
  mutate(year = year(date))

Code

sebal_annual %>%
  st_transform(32723) %>%
  mutate(
    X = st_coordinates(st_centroid(geometry))[,1],
    Y = st_coordinates(st_centroid(geometry))[,2]
  ) %>%
  ggplot() +
  geom_sf(aes(fill = value), color = NA) +
  geom_label(
    data = . %>% filter(value >= 3500),
    aes(x = X, y = Y, label = NAME_1, fill = value),
    color = "white",        
    size = 2.5,
    fontface = "bold",
    label.size = 0,         
    alpha = 0.8             
  ) +
  scale_fill_gradientn(
    colours = c("#87CEFA", "#4682B4", "#08306B"),
    name = "Evapotranspiration (mm)"
  ) +
  facet_wrap(~year, ncol = 4, nrow = 2) +
  theme_minimal(base_size = 12) +
  theme(
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    strip.text = element_text(size = 11, face = "bold", color = "#333333"),
    axis.text = element_blank(),
    axis.title = element_blank(),
    legend.position = "bottom",
    legend.key.width = unit(1, "cm"),
    legend.title = element_text(size = 10, face = "bold"),
    legend.text = element_text(size = 9)
  )

Spatial patterns derived from MODIS-based SEBAL data show relatively stable evapotranspiration rates over time, with consistently higher values observed in northern and northeastern states such as Amapá, Pará, and Maranhão. These regions exhibit strong vegetative activity and higher surface energy fluxes due to dense forest cover and humid climatic conditions.

3.6 How can we see variation across variables?

Bivariate maps allow us to visualize how two environmental or socioeconomic variables interact spatially over time.
Instead of analyzing them independently, this approach overlays both indicators using a combined color scale, where each color represents a specific combination of quartiles of the two variables.

In this analysis, we use bivariate color maps to explore spatial associations between indicators such as:

• Evapotranspiration (biophysical factor)
  
• Night-time lights (human activity proxy)
  
• Forest loss, precipitation, or accessibility (contextual drivers)

Code

plot_bivariate_facet <- function(
  data1, data2,
  var1_label = "Variable 1",
  var2_label = "Variable 2",
  years = NULL,
  scale_type = "quartile",     
  pal_var1 = NULL,  
  pal_var2 = NULL,  
  ncol_facets = 4,
  legend_pos = "right"
) {
  library(sf)
  library(dplyr)
  library(ggplot2)
  library(scales)
  library(cowplot)
  
  if (is.null(years)) years <- sort(unique(data1$year))

  biv_data <- data1 %>%
    select(NAME_1, geometry, year, var1 = value) %>%
    left_join(
      data2 %>%
        st_drop_geometry() %>%
        select(NAME_1, year, var2 = value),
      by = c("NAME_1", "year")
    ) %>%
    filter(year %in% years) %>%
    st_as_sf()

  n_breaks <- switch(scale_type,
                     "tertile"  = 3,
                     "quartile" = 4,
                     "quintile" = 5,
                     "decile"   = 10,
                     4)

  probs <- seq(0, 1, 1 / n_breaks)
  v1_qtiles <- quantile(biv_data$var1, probs = probs, na.rm = TRUE)
  v2_qtiles <- quantile(biv_data$var2, probs = probs, na.rm = TRUE)

  biv_data <- biv_data %>%
    mutate(
      v1_class = cut(var1, breaks = v1_qtiles, include.lowest = TRUE, labels = 1:n_breaks),
      v2_class = cut(var2, breaks = v2_qtiles, include.lowest = TRUE, labels = 1:n_breaks),
      bi_class = paste0(v1_class, "-", v2_class)
    )

  make_bi_pal <- function(pal1, pal2, n = n_breaks) {
    col1 <- colorRampPalette(pal1)(n)
    col2 <- colorRampPalette(pal2)(n)
    outer(seq_len(n), seq_len(n), Vectorize(function(i, j) {
      rgb(
        (col2rgb(col1[i])[1] + col2rgb(col2[j])[1]) / 510,
        (col2rgb(col1[i])[2] + col2rgb(col2[j])[2]) / 510,
        (col2rgb(col1[i])[3] + col2rgb(col2[j])[3]) / 510
      )
    }))
  }

  bi_pal <- make_bi_pal(pal_var1, pal_var2)
  bi_colors <- as.vector(t(bi_pal))

  color_map <- expand.grid(v1_class = 1:n_breaks, v2_class = 1:n_breaks) %>%
    mutate(bi_class = paste0(v1_class, "-", v2_class),
           fill = bi_colors)

  biv_data <- left_join(biv_data, color_map, by = "bi_class")

  map_faceted <- ggplot(biv_data) +
    geom_sf(aes(fill = fill), color = "white", size = 0.05) +
    scale_fill_identity() +
    facet_wrap(~year, ncol = ncol_facets) +
    theme_void() +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
      strip.text = element_text(face = "bold", size = 10),
      panel.spacing = unit(0.4, "lines")
    )
  legend_df <- expand.grid(v1_class = 1:n_breaks, v2_class = 1:n_breaks) %>%
    mutate(fill = as.vector(t(bi_pal)))

  legend_plot <- ggplot(legend_df, aes(x = v1_class, y = v2_class, fill = fill)) +
    geom_tile(color = "white") +
    scale_fill_identity() +
    coord_equal() +
    theme_void() +
    theme(plot.margin = margin(10, 10, 10, 10)) +
    annotate("segment", x = 1, xend = n_breaks + 0.5, y = 0.3, yend = 0.3,
             arrow = arrow(length = unit(0.25, "cm")), color = "black", size = 0.6) +
    annotate("segment", x = 0.3, xend = 0.3, y = 1, yend = n_breaks + 0.5,
             arrow = arrow(length = unit(0.25, "cm")), color = "black", size = 0.6) +
    annotate("text", x = (n_breaks + 1) / 2, y = -0.2,
             label = paste0(var1_label, " (Q1–Q", n_breaks, ")"),
             size = 2, fontface = "bold") +
    annotate("text", x = 0, y = (n_breaks + 1) / 2,
             label = paste0(var2_label, " (Q1–Q", n_breaks, ")"),
             size = 2, fontface = "bold", angle = 90)
  if (legend_pos == "right") {
    out <- cowplot::plot_grid(map_faceted, legend_plot,
                              ncol = 2,
                              rel_widths = c(0.8, 0.22),
                              align = "h")
  } else {
    out <- cowplot::ggdraw() +
      draw_plot(map_faceted, 0, 0, 1, 1) +
      draw_plot(legend_plot, 0.75, 0.05, 0.25, 0.25)
  }

  return(out)
}

plot_bivariate_facet(
  sebal_annual, forest_loss,
  var1_label = "Evapotranspiration",
  var2_label = "Forest loss",
  scale_type = "quartil",
  pal_var1 = c("#87CEFA", "#4682B4", "#08306B"),     # azul
  pal_var2 = c("#00441B", "#1B7837", "#A6DBA0", "#FEE08B", "#E08214", "#B35806")      
)