Code
library(ipumsr)
library(sf)
library(dplyr)
library(tidyr)
library(zoo)
library(lubridate)
library(terra)
library(fs)
library(stringr)
library(purrr)
library(rhdf5)
library(ggplot2)
library(ggspatial)
# Read in H5 files ----
h5_files <- list.files("data_local/evi", full.names = TRUE, pattern = "\\.h5$")
tile_codes <- unique(str_extract(h5_files, "h[0-9]{2}v[0-9]{2}"))
tiles <- map(
tile_codes,
function(code) h5_files[str_detect(h5_files, code)]
)
# Read in metadata from HDF files
tile_metadata <- map(
tiles,
function(t) h5read(t[1], name = "//HDFEOS INFORMATION/StructMetadata.0")
)
# Assign sinusoidal projection information from object tile_metadata
sinu_proj <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"
# Gather upper and lower limits for downloaded tiles
ul <- map(
tile_metadata,
~ str_match(.x, "UpperLeftPointMtrs=\\((.*?)\\)")[, 2]
)
ul <- map(str_split(ul, ","), as.numeric)
lr <- map_chr(
tile_metadata,
~ str_match(.x, "LowerRightMtrs=\\((.*?)\\)")[, 2]
)
lr <- map(str_split(lr, ","), as.numeric)
# Save coordinate extents for tiles
tile_ext <- map2(ul, lr, function(u, l) ext(u[1], l[1], l[2], u[2]))
evi_tiles <- map2(
tiles,
tile_ext,
function(tile, ext) {
# Load raster for the input tile. We select the EVI subdataset
r <- rast(
tile,
subds = "//HDFEOS/GRIDS/VIIRS_Grid_16Day_VI_500m/Data_Fields/500_m_16_days_EVI",
noflip = TRUE
)
crs(r) <- sinu_proj # Attach sinusoidal projection defined above
ext(r) <- ext # Attach this tile's extent
r # Return the updated raster for this tile
}
)
# Mosaic tiles together
evi_mosaic <- reduce(evi_tiles, mosaic)
# Replace EVI values of less than -1 to missing
m <- matrix(c(-Inf, -1, NA), nrow = 1)
evi_mosaic <- classify(evi_mosaic, m)
# Read in Uganda shapefiles and transform the projection to match the EVI data
ug_borders <- ipumsr::read_ipums_sf("shapefiles/geo1_ug1991_2014.zip") |>
st_make_valid() |> # Fix minor border inconsistencies
st_union() |>
st_transform(crs(evi_mosaic))
ug_evi_crop <- crop(evi_mosaic, ug_borders)
ug_evi_mask <- mask(ug_evi_crop, vect(ug_borders))
# Create color pallete for EVI
evi_pal <- list(
pal = c(
"#fdfbdc",
"#f1f4b7",
"#d3ef9f",
"#a5da8d",
"#6cc275",
"#51a55b",
"#397e43",
"#2d673a",
"#1d472e"
),
values = c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 1)
)
#Create plot of EVI for Uganda
uganda_evi_plot <- ggplot() +
layer_spatial(ug_evi_mask[[1]]) +
scale_fill_gradientn(
colors = evi_pal$pal,
values = evi_pal$values,
limits = c(0, 1),
na.value = "transparent"
) +
labs(title = "Uganda Enhanced Vegetative Index (EVI) March 29-April 6, 2024", subtitle = "NASA VIIRS, 500m resolution", fill = "EVI")
uganda_evi_plot
