Case 03 – Peru(🇵🇪)

Exploring environmental and accessibility drivers of malaria risk across the Peruvian Amazon

Context

Between 2010 and 2023, the Peruvian Amazon has experienced accelerated deforestation and recurrent malaria outbreaks. Forest loss in Loreto, Ucayali, and Madre de Dios has modified local temperature regimes and humidity patterns, favoring Anopheles vector proliferation. Accessibility to diagnostic and treatment facilities remains highly heterogeneous, especially along riverine and forest-edge communities.

1 Objective

Visualize the relationship between association between malaria incidence count, temperature, forest loss, and travel time to health centers across the Peruvian Amazon

2 Functions from land4health()

Variable	Function	Description
Malaria	l4h_layers_available_malaria()	Provides a tidy listing of all malaria data layers in the Malaria Atlas Project GeoServer by species and metric.
Forest loss	l4h_forest_loss()	Calculates forest loss within a user-defined region for a specified year range.
Travel time	l4h_travel_time()	Retrieves the travel time raster (in minutes) to the nearest health facility or city using the Oxford Global Accessibility dataset.
Surface temperature	l4h_surface_temp()	Retrieves daily Land Surface Temperature (MODIS MOD11A1) data from Google Earth Engine.

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()

library(land4health)
library(geodata)
library(ggplot2)
library(tidyverse)
library(sf)
library(terra)
library(tidyterra)
library(malariaAtlas)
library(RColorBrewer)
library(lubridate)
library(leaflet.extras2)
library(purrr)
library(mapview)
# 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/yomali_ferreyra 

 ✔ 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. Only the Amazonian departments (Loreto and Ucayali) are selected as the study area.

peru_sf <- geodata::gadm(country = "PE", level = 2, path = tempdir()) |>
st_as_sf()
sel_deps <- c("Loreto", "Ucayali")
amazon_sf <- peru_sf |> filter(NAME_1 %in% sel_deps)

amazon_sf<-rmapshaper::ms_simplify(amazon_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 Malaria

First, we retrieve annual Plasmodium vivax malaria incidence data from the Malaria Atlas Project (MAP) GeoServer using the l4h_layers_available_malaria() and getRaster() functions. Each raster layer represents the estimated number of malaria cases per 5 km × 5 km grid cell.The rasters are then stacked across multiple years (2016–2022) and clipped to the Peruvian Amazon shapefile for subsequent spatial analyses.

layer_malaria<-l4h_layers_available_malaria(year = 2022, 
                                            species = "plasmodium vivax", 
                                            measure = "incidence_count")
dataset<-layer_malaria$dataset_id
years <- 2016:2021

layers <- map(years, function(y) {
  getRaster(
    dataset_id = dataset,  
    year = y,              
    shp = amazon_sf        
  )
})

# Each layer represents a different year (2016–2022)
malaria_incidence<- rast(layers)

malaria_zonal <- terra::extract(
  malaria_incidence,
  amazon_sf[c("NAME_1", "NAME_2")],  
  fun = sum,
  na.rm = TRUE,
  bind = TRUE                        
)

malaria_sf <- sf::st_as_sf(malaria_zonal)
names(malaria_sf)[3:7] <- as.character(2016:2021)

malaria_long <- malaria_sf %>%
  pivot_longer(
    cols = matches("^[0-9]{4}$"),  
    names_to = "year",
    values_to = "incidence"
  ) 

malaria_plot<-malaria_long%>%
  group_by(NAME_1, year) %>%
  summarise(total_cases = sum(incidence, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(year = as.numeric(year))


ggplot(malaria_plot, aes(x = year, y = total_cases / 1000, fill = NAME_1)) +
  geom_col(position = "dodge", color = "black", width = 0.7, alpha = 0.9) +
  scale_fill_manual(
    values = c("Loreto" = "#b2182b", "Ucayali" = "#2166ac"),
    name = NULL
  ) +
  labs(
    x = "Malaria cases (thousands)",
    y = "Year"
  ) +
  theme_minimal(base_size = 8) +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.major.y = element_line(color = "gray85", linewidth = 0.3),
    axis.title = element_text(face = "bold"),
    plot.title = element_text(face = "bold", size = 14),
    legend.position = "top",
    legend.text = element_text(size = 8),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

3.4 Forest loss

In this part, we assess forest cover loss across the Peruvian Amazon using the l4h_forest_loss() function. This function computes the annual area of forest loss within a defined region, based on the Hansen Global Forest Change dataset (Hansen et al., 2013). Forest loss is defined as a stand‐replacement disturbance, meaning the conversion from forest to non‐forest cover. We will extract annual loss data for the period 2016-2022

# Loss of vegetation cover
forest_loss <-  l4h_forest_loss(
  from = "2016-01-01",
  to = "2022-12-31",
  region = amazon_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) +
  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)
  )

3.5 Travel time

In this part, we estimate the travel time to healthcare facilities across the Peruvian Amazon using the l4h_travel_time() function. This function extracts data from the Oxford Global Map of Accessibility, which models the time (in minutes) required to reach the nearest major city.

travel_time<- l4h_travel_time(
  region        = amazon_sf,
  destination   = "cities",
  sf =TRUE)

malaria_avg <- malaria_long %>%
  group_by(NAME_1, NAME_2) %>%
  summarise(mean_incidence = mean(incidence, na.rm = TRUE),
            geometry = st_union(geometry)) %>%
  ungroup()

malaria_centroids <- malaria_avg %>%
  st_centroid() %>%
  st_transform(st_crs(travel_time)) %>%  
  mutate(lon = st_coordinates(.)[,1],
         lat = st_coordinates(.)[,2])

ggplot() +
  geom_sf(
    data = travel_time,
    aes(fill = city_access),
    color = NA
  ) +
  scale_fill_viridis_c(
    option = "magma",
    direction = -1,
    name = "Travel time (minutes)"
  ) +
  coord_sf(xlim = c(-78, -69.5), ylim = c(-11.5, -0.2)) +
  labs(
    x = NULL, y = NULL
  ) +
  theme_minimal(base_size = 13) +
  theme(
    panel.background = element_rect(fill = "white", color = NA),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank(),
    legend.position = "right",
    legend.title = element_text(size = 8, face = "bold"),
    legend.text = element_text(size = 9)
  )