1 Project Overview

This RPubs-ready practical analysis uses the Sababi Institute six-region Somaliland shapefile and combines:

  • Administrative regions
  • Average annual temperature from WorldClim
  • Elevation / DEM
  • Next 5-day rainfall forecast from CHIRPS3-GEFS

The original project script specifies that the maps should be displayed in the RStudio Plot Pane. In this R Markdown/RPubs version, the same plots are also explicitly rendered inside the knitted HTML document, so the maps and charts become visible on the published RPubs page.

2 Important Before Knitting

Place the folder:

SL_Six_Regions

inside your Downloads folder. It should contain the shapefile and its companion files (.shp, .dbf, .shx, .prj).

An active internet connection is required because WorldClim, elevation, and CHIRPS3-GEFS data are downloaded when needed.

3 1. Packages and Settings

packages <- c(
  "sf", "terra", "ggplot2", "dplyr", "stringr",
  "scales", "viridis", "geodata", "elevatr",
  "raster", "rvest", "httr"
)

missing_pkgs <- packages[!packages %in% rownames(installed.packages())]

if (length(missing_pkgs) > 0) {
  install.packages(missing_pkgs, dependencies = TRUE)
}

invisible(lapply(packages, library, character.only = TRUE))

options(stringsAsFactors = FALSE)
options(timeout = 1800)

4 2. Run the Complete GIS Analysis

The code below is based on the supplied project script. It keeps the original data-processing workflow and prints the maps during knitting so that they are captured by RPubs.

# ============================================================
# SOMALILAND GIS CLIMATE & 5-DAY RAINFALL EARLY WARNING MAPS
# RStudio | Sababi Institute Six-Region Shapefile
# Author: Yahye S. Rageh
#
# FOUR MAPS:
#   1. Somaliland Six Regions
#   2. Average Annual Temperature (WorldClim)
#   3. Elevation (DEM)
#   4. Next 5-Day Rainfall Forecast (CHIRPS3-GEFS)
#
# FIXES APPLIED IN THIS VERSION (vs. the original script):
#   1. Downloads-folder detection now works on Windows, Mac
#      and Linux (the original only worked on Windows because
#      it relied on Sys.getenv("USERPROFILE"), which is blank
#      everywhere else).
#   2. CHIRPS3-GEFS archive path fixed to "05_day" (not
#      "5_day") -- data.chc.ucsb.edu returns 404 on the old path.
#   3. The archive directory listing is now parsed with rvest
#      (real HTML parsing) instead of a raw regex over
#      readLines(), and requests send a normal browser User-
#      Agent -- UCSB's server occasionally blocks R's default
#      agent / plain readLines() on a URL.
#   4. Package list updated to include "raster" and "rvest",
#      which elevatr / the archive-scraping step need.
#   5. get_elev_raster() output is wrapped safely and checked
#      before conversion to SpatRaster.
#   6. Every download step now fails with a clear, actionable
#      message instead of a cryptic low-level error.
#   7. Small robustness fixes: safer year/date parsing, WGS84
#      fallback only when CRS truly missing, and a guard so the
#      script does not silently continue with 0 regions.
# ============================================================

rm(list = ls())
graphics.off()

options(stringsAsFactors = FALSE)
options(timeout = 1800)

# ============================================================
# 1. PACKAGES
# ============================================================

packages <- c(
  "sf",
  "terra",
  "ggplot2",
  "dplyr",
  "stringr",
  "scales",
  "viridis",
  "geodata",
  "elevatr",
  "raster",   # elevatr's get_elev_raster() depends on this
  "rvest",    # robust HTML parsing of the CHIRPS-GEFS archive
  "httr"      # sends a proper User-Agent when reading the archive
)

missing_pkgs <- packages[!packages %in% rownames(installed.packages())]

if (length(missing_pkgs) > 0) {
  message("Installing missing packages: ", paste(missing_pkgs, collapse = ", "))
  install.packages(missing_pkgs, dependencies = TRUE)
}

invisible(lapply(packages, library, character.only = TRUE))

# ============================================================
# 2. PROJECT FOLDERS
# ============================================================

project_dir <- getwd()

root_dir <- file.path(project_dir, "Somaliland_Climate_Early_Warning")

dir.create(root_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(file.path(root_dir, "data", "climate"), recursive = TRUE, showWarnings = FALSE)
dir.create(file.path(root_dir, "data", "dem"), recursive = TRUE, showWarnings = FALSE)
dir.create(file.path(root_dir, "data", "forecast"), recursive = TRUE, showWarnings = FALSE)
dir.create(file.path(root_dir, "outputs", "maps"), recursive = TRUE, showWarnings = FALSE)
dir.create(file.path(root_dir, "outputs", "tables"), recursive = TRUE, showWarnings = FALSE)
dir.create(file.path(root_dir, "outputs", "rasters"), recursive = TRUE, showWarnings = FALSE)

out_map <- file.path(root_dir, "outputs", "maps")
out_tab <- file.path(root_dir, "outputs", "tables")
out_ras <- file.path(root_dir, "outputs", "rasters")

# ============================================================
# 3. SABABI INSTITUTE SHAPEFILE
# ============================================================
# FIX: USERPROFILE only exists on Windows. On Mac/Linux it is
# blank, so the original script could never find the folder on
# those systems. This now works on all three.

downloads_dir <- if (.Platform$OS.type == "windows") {
  win_home <- Sys.getenv("USERPROFILE")
  if (!nzchar(win_home)) win_home <- path.expand("~")
  file.path(win_home, "Downloads")
} else {
  file.path(path.expand("~"), "Downloads")
}

sababi_folder <- file.path(downloads_dir, "SL_Six_Regions")

if (!dir.exists(sababi_folder)) {
  stop(
    "\nSABABI FOLDER NOT FOUND.\n\n",
    "Expected:\n", sababi_folder, "\n\n",
    "Place the 'SL_Six_Regions' folder (containing SL_SixRegions.shp\n",
    "and its .dbf/.shx/.prj siblings) inside your Downloads folder,\n",
    "or edit 'sababi_folder' above to point at the correct location."
  )
}

message("\nSababi folder found:")
## 
## Sababi folder found:
message(sababi_folder)
## C:\Users\hp/Downloads/SL_Six_Regions
shp_files <- list.files(
  sababi_folder,
  pattern = "\\.shp$",
  full.names = TRUE,
  recursive = TRUE,
  ignore.case = TRUE
)

if (length(shp_files) == 0) {
  stop("\nNo .shp file was found inside:\n", sababi_folder)
}

message("\nShapefile(s) found:")
## 
## Shapefile(s) found:
print(shp_files)

preferred <- shp_files[grepl("SL_SixRegions", basename(shp_files), ignore.case = TRUE)]
boundary_shp <- if (length(preferred) > 0) preferred[1] else shp_files[1]

message("\nUsing Sababi shapefile:")
## 
## Using Sababi shapefile:
message(boundary_shp)
## C:\Users\hp/Downloads/SL_Six_Regions/SL_SixRegions.shp
# ============================================================
# 4. READ SABABI SIX REGIONS
# ============================================================

somaliland <- sf::st_read(boundary_shp, quiet = TRUE)

if (nrow(somaliland) == 0) {
  stop("\nThe shapefile loaded with 0 features. Check the file is not corrupted.")
}

if (nrow(somaliland) != 6) {
  warning("\nExpected 6 Sababi regions, but found ", nrow(somaliland), " features.")
}

message("\nNumber of Sababi regions: ", nrow(somaliland))
## 
## Number of Sababi regions: 6
message("\nAvailable fields:")
## 
## Available fields:
print(names(somaliland))

if (!"admin1Name" %in% names(somaliland)) {
  stop(
    "\nThe field 'admin1Name' was not found.\n",
    "Available fields:\n", paste(names(somaliland), collapse = ", ")
  )
}

somaliland$Region_raw <- stringr::str_squish(as.character(somaliland$admin1Name))

message("\nOriginal Sababi region names:")
## 
## Original Sababi region names:
print(somaliland$Region_raw)

# ============================================================
# 5. STANDARDIZE REGION NAMES
# ============================================================

somaliland$Region <- dplyr::case_when(
  str_detect(str_to_lower(somaliland$Region_raw), "awdal") ~ "Awdal",
  str_detect(str_to_lower(somaliland$Region_raw), "maroodi|woqooyi|galbeed") ~ "Maroodi Jeex",
  str_detect(str_to_lower(somaliland$Region_raw), "sahil|saxil") ~ "Sahil",
  str_detect(str_to_lower(somaliland$Region_raw), "sanaag") ~ "Sanaag",
  str_detect(str_to_lower(somaliland$Region_raw), "sool") ~ "Sool",
  str_detect(str_to_lower(somaliland$Region_raw), "togdheer|togdher") ~ "Togdheer",
  TRUE ~ somaliland$Region_raw
)

message("\nFinal region names:")
## 
## Final region names:
print(somaliland$Region)

# ============================================================
# 6. GEOMETRY CLEANING
# ============================================================

if (is.na(sf::st_crs(somaliland))) {
  warning("Shapefile CRS is missing. Assuming EPSG:4326.")
  sf::st_crs(somaliland) <- 4326
}

somaliland <- somaliland |>
  sf::st_make_valid() |>
  sf::st_transform(4326)

somaliland <- somaliland[!sf::st_is_empty(somaliland), ]

if (nrow(somaliland) == 0) {
  stop("No valid Somaliland geometries remain after cleaning.")
}

som_vect <- terra::vect(somaliland)

# ============================================================
# 7. SAFE REGION LABELS
# ============================================================
# Project to an equal-area CRS first, then compute point-on-
# surface, then transform back -- avoids the lon/lat warning.

label_points <- somaliland |>
  sf::st_transform(6933) |>
  sf::st_point_on_surface() |>
  sf::st_transform(4326)
## Warning: st_point_on_surface assumes attributes are constant over geometries
# ============================================================
# 8. MAP EXTENT
# ============================================================

bb <- sf::st_bbox(somaliland)
pad_x <- (bb["xmax"] - bb["xmin"]) * 0.025
pad_y <- (bb["ymax"] - bb["ymin"]) * 0.025

map_coord <- ggplot2::coord_sf(
  xlim = c(bb["xmin"] - pad_x, bb["xmax"] + pad_x),
  ylim = c(bb["ymin"] - pad_y, bb["ymax"] + pad_y),
  expand = FALSE,
  datum = NA
)

# ============================================================
# 9. PROFESSIONAL GIS THEME
# ============================================================

gis_theme <- ggplot2::theme_minimal(base_size = 12) +
  ggplot2::theme(
    plot.title = element_text(size = 20, face = "bold", hjust = 0.5, colour = "#08306B"),
    plot.subtitle = element_text(size = 11, hjust = 0.5, colour = "#4D4D4D"),
    plot.caption = element_text(size = 8.5, colour = "grey35", hjust = 0),
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 9),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank(),
    panel.background = element_rect(fill = "#F7FAFC", colour = NA),
    plot.background = element_rect(fill = "white", colour = NA),
    plot.margin = margin(10, 10, 10, 10)
  )

save_map <- function(p, filename, width = 13, height = 9, dpi = 600) {
  ggplot2::ggsave(
    filename = file.path(out_map, filename),
    plot = p,
    width = width,
    height = height,
    dpi = dpi,
    bg = "white"
  )
}

# ============================================================
# MAP 1 -- SOMALILAND SIX REGIONS
# ============================================================

message("\n============================================================")
## 
## ============================================================
message("MAP 1 -- SOMALILAND SIX REGIONS")
## MAP 1 -- SOMALILAND SIX REGIONS
message("============================================================")
## ============================================================
regions_map <- ggplot2::ggplot() +
  geom_sf(data = somaliland, aes(fill = Region), colour = "white", linewidth = 0.8) +
  geom_sf_text(data = label_points, aes(label = Region), size = 3.5, fontface = "bold", colour = "#111111") +
  scale_fill_viridis_d(option = "turbo", guide = "none") +
  labs(
    title = "SOMALILAND ADMINISTRATIVE REGIONS",
    subtitle = "Six current Somaliland regions",
    caption = "Boundary source: Sababi Institute"
  ) +
  map_coord +
  gis_theme

print(regions_map)
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
save_map(regions_map, "01_Somaliland_Regions.png")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
# ============================================================
# MAP 2 -- AVERAGE ANNUAL TEMPERATURE
# ============================================================

message("\n============================================================")
## 
## ============================================================
message("MAP 2 -- AVERAGE ANNUAL TEMPERATURE")
## MAP 2 -- AVERAGE ANNUAL TEMPERATURE
message("============================================================")
## ============================================================
message("\nDownloading WorldClim Tmin...")
## 
## Downloading WorldClim Tmin...
tmin <- geodata::worldclim_country(
  country = "Somalia",
  var = "tmin",
  path = file.path(root_dir, "data", "climate")
)
## Cached as: C:/Users/hp/Downloads/Somaliland_Climate_Early_Warning/data/climate/climate/wc2.1_country/SOM_wc2.1_30s_tmin.tif
message("Downloading WorldClim Tmax...")
## Downloading WorldClim Tmax...
tmax <- geodata::worldclim_country(
  country = "Somalia",
  var = "tmax",
  path = file.path(root_dir, "data", "climate")
)
## Cached as: C:/Users/hp/Downloads/Somaliland_Climate_Early_Warning/data/climate/climate/wc2.1_country/SOM_wc2.1_30s_tmax.tif
# WorldClim monthly temperature is stored in whole degrees C
# (recent releases) -- compute the 12-month mean per cell, then
# average Tmin and Tmax into a single annual mean temperature.
tmin_mean <- terra::app(tmin, mean, na.rm = TRUE)
tmax_mean <- terra::app(tmax, mean, na.rm = TRUE)

temperature_raster <- (tmin_mean + tmax_mean) / 2
temperature_raster <- terra::crop(temperature_raster, som_vect)
temperature_raster <- terra::mask(temperature_raster, som_vect)
names(temperature_raster) <- "Temperature_C"

temperature_df <- as.data.frame(temperature_raster, xy = TRUE, na.rm = TRUE)
names(temperature_df)[3] <- "Temperature_C"

temperature_map <- ggplot2::ggplot() +
  geom_raster(data = temperature_df, aes(x = x, y = y, fill = Temperature_C)) +
  geom_sf(data = somaliland, fill = NA, colour = "white", linewidth = 0.8) +
  geom_sf_text(data = label_points, aes(label = Region), size = 2.8, fontface = "bold", colour = "#111111") +
  scale_fill_gradientn(
    colours = c("#313695", "#4575B4", "#74ADD1", "#ABD9E9", "#FFFFBF", "#FDAE61", "#F46D43", "#D73027"),
    name = "Temperature\n(\u00b0C)",
    labels = function(x) paste0(round(x, 1), "\u00b0C")
  ) +
  labs(
    title = "AVERAGE ANNUAL TEMPERATURE",
    subtitle = "WorldClim spatial temperature pattern",
    caption = "Source: WorldClim v2.1 | Tmin + Tmax cell-by-cell mean"
  ) +
  map_coord +
  gis_theme

print(temperature_map)
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
save_map(temperature_map, "02_Average_Annual_Temperature.png")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not give correct results for longitude/latitude data
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
temp_extract <- terra::extract(temperature_raster, som_vect, fun = mean, na.rm = TRUE)

temperature_summary <- data.frame(
  Region = somaliland$Region,
  Mean_Temperature_C = round(temp_extract[, 2], 1)
)

write.csv(
  temperature_summary,
  file.path(out_tab, "Somaliland_Temperature_By_Region.csv"),
  row.names = FALSE
)

# ============================================================
# MAP 3 -- ELEVATION
# ============================================================

message("\n============================================================")
## 
## ============================================================
message("MAP 3 -- ELEVATION")
## MAP 3 -- ELEVATION
message("============================================================")
## ============================================================
dem_file <- file.path(root_dir, "data", "dem", "Somaliland_DEM.tif")

if (!file.exists(dem_file)) {
  
  message("\nDownloading DEM through elevatr...")
  
  dem_try <- tryCatch(
    elevatr::get_elev_raster(locations = somaliland, z = 7, clip = "locations"),
    error = function(e) {
      message("DEM download failed: ", e$message)
      NULL
    }
  )
  
  if (is.null(dem_try)) {
    stop("\nElevation download failed.\nCheck your internet connection and rerun.")
  }
  
  dem <- terra::rast(dem_try)
  terra::writeRaster(dem, dem_file, overwrite = TRUE)
  
} else {
  message("\nExisting DEM found. Reusing it.")
  dem <- terra::rast(dem_file)
}
## 
## Downloading DEM through elevatr...
## Mosaicing & Projecting
## Clipping DEM to locations
## Note: Elevation units are in meters.
dem <- terra::project(dem, terra::crs(som_vect))
dem <- terra::crop(dem, som_vect)
dem <- terra::mask(dem, som_vect)
names(dem) <- "Elevation_m"

elevation_df <- as.data.frame(dem, xy = TRUE, na.rm = TRUE)
names(elevation_df)[3] <- "Elevation_m"

elevation_map <- ggplot2::ggplot() +
  geom_raster(data = elevation_df, aes(x = x, y = y, fill = Elevation_m)) +
  geom_sf(data = somaliland, fill = NA, colour = "white", linewidth = 0.8) +
  geom_sf_text(data = label_points, aes(label = Region), size = 2.8, fontface = "bold", colour = "#111111") +
  scale_fill_gradientn(
    colours = c("#0B3C5D", "#328CC1", "#99C24D", "#E6AF2E", "#D95D39", "#7F2704"),
    name = "Elevation\n(m)",
    labels = scales::label_number(accuracy = 1)
  ) +
  labs(
    title = "DIGITAL ELEVATION MODEL",
    subtitle = "Elevation distribution across Somaliland",
    caption = "Source: Elevation tiles accessed through elevatr (AWS Terrain Tiles)"
  ) +
  map_coord +
  gis_theme

print(elevation_map)
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
save_map(elevation_map, "03_Elevation.png")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
elev_extract <- terra::extract(dem, som_vect, fun = mean, na.rm = TRUE)

elevation_summary <- data.frame(
  Region = somaliland$Region,
  Mean_Elevation_m = round(elev_extract[, 2], 1)
)

write.csv(
  elevation_summary,
  file.path(out_tab, "Somaliland_Elevation_By_Region.csv"),
  row.names = FALSE
)

# ============================================================
# MAP 4 -- NEXT 5-DAY RAINFALL FORECAST
# ============================================================
# CHIRPS3-GEFS 5-day files are 5-day total precipitation
# forecast products at 0.05 degree resolution.
#
# FIX #1: correct archive path is ".../v3/05_day/..." not
#         ".../v3/5_day/..." (the old path 404s).
# FIX #2: the directory listing is now parsed with rvest, and
#         the HTTP request sends a normal browser User-Agent.
#         Plain readLines(url) on this server sometimes returns
#         an empty/blocked page with no informative error, which
#         made the old regex-based parser fail silently.

message("\n============================================================")
## 
## ============================================================
message("MAP 4 -- NEXT 5-DAY RAINFALL FORECAST")
## MAP 4 -- NEXT 5-DAY RAINFALL FORECAST
message("============================================================")
## ============================================================
download_latest_chirps3_gefs_5day <- function() {
  
  base_url <- "https://data.chc.ucsb.edu/products/CHIRPS-GEFS/v3/05_day/global/data/"
  ua <- httr::user_agent("Mozilla/5.0 (Somaliland-GIS-EarlyWarning-Script)")
  
  message("\nCHIRPS3-GEFS archive:")
  message(base_url)
  
  get_links <- function(url) {
    resp <- tryCatch(
      httr::GET(url, ua, httr::timeout(120)),
      error = function(e) {
        stop(
          "\nCould not reach CHIRPS3-GEFS archive.\n\nURL:\n", url,
          "\n\nError:\n", e$message
        )
      }
    )
    
    if (httr::status_code(resp) != 200) {
      stop(
        "\nCHIRPS3-GEFS archive returned HTTP ", httr::status_code(resp),
        " for:\n", url
      )
    }
    
    page <- httr::content(resp, as = "text", encoding = "UTF-8")
    html <- rvest::read_html(page)
    rvest::html_attr(rvest::html_elements(html, "a"), "href")
  }
  
  # ---- find year folders -------------------------------------------------
  links <- get_links(base_url)
  year_links <- links[grepl("^[0-9]{4}/?$", links)]
  years <- as.integer(gsub("/", "", year_links))
  years <- years[!is.na(years)]
  
  if (length(years) == 0) {
    stop("\nNo year folders found in CHIRPS3-GEFS 05-day archive.")
  }
  
  latest_year <- max(years)
  message("Latest archive year: ", latest_year)
  
  # ---- find latest forecast file in that year ----------------------------
  year_url <- paste0(base_url, latest_year, "/")
  year_links <- get_links(year_url)
  
  files <- unique(year_links[grepl("^c3g_[0-9]{4}\\.[0-9]{2}\\.[0-9]{2}\\.tif$", year_links)])
  
  if (length(files) == 0) {
    stop("\nNo c3g 5-day GeoTIFF files found in:\n", year_url)
  }
  
  latest_file <- sort(files, decreasing = TRUE)[1]
  forecast_url <- paste0(year_url, latest_file)
  local_file <- file.path(root_dir, "data", "forecast", latest_file)
  
  message("\nLatest CHIRPS3-GEFS forecast:")
  message(latest_file)
  message("\nForecast URL:")
  message(forecast_url)
  
  # ---- download if necessary ---------------------------------------------
  if (!file.exists(local_file)) {
    
    message("\nDownloading forecast GeoTIFF (may take a while)...")
    
    tryCatch(
      {
        resp <- httr::GET(
          forecast_url,
          ua,
          httr::write_disk(local_file, overwrite = TRUE),
          httr::progress()
        )
        if (httr::status_code(resp) != 200) {
          stop("HTTP status ", httr::status_code(resp))
        }
      },
      error = function(e) {
        if (file.exists(local_file)) unlink(local_file)
        stop("\nCHIRPS3-GEFS download failed.\n\n", e$message)
      }
    )
    
    if (!file.exists(local_file)) {
      stop("\nThe forecast file was not created.")
    }
    
  } else {
    message("\nForecast file already exists locally.")
  }
  
  # ---- validate ------------------------------------------------------------
  file_size_mb <- file.info(local_file)$size / (1024^2)
  message("\nDownloaded file size: ", round(file_size_mb, 1), " MB")
  
  if (is.na(file_size_mb) || file_size_mb < 1) {
    unlink(local_file)
    stop("\nDownloaded file is invalid or incomplete.\nPlease rerun the script.")
  }
  
  list(file = local_file, filename = latest_file, url = forecast_url)
}

forecast_5 <- download_latest_chirps3_gefs_5day()
## 
## CHIRPS3-GEFS archive:
## https://data.chc.ucsb.edu/products/CHIRPS-GEFS/v3/05_day/global/data/
## Latest archive year: 2026
## 
## Latest CHIRPS3-GEFS forecast:
## c3g_2026.09.15.tif
## 
## Forecast URL:
## https://data.chc.ucsb.edu/products/CHIRPS-GEFS/v3/05_day/global/data/2026/c3g_2026.09.15.tif
## 
## Downloading forecast GeoTIFF (may take a while)...
## 
## Downloaded file size: 61.9 MB
# ---- read + reproject + clip forecast raster -------------------------------

message("\nReading rainfall forecast raster...")
## 
## Reading rainfall forecast raster...
forecast_raster <- terra::rast(forecast_5$file)

message("Forecast raster CRS:")
## Forecast raster CRS:
print(terra::crs(forecast_raster))

if (!terra::same.crs(forecast_raster, som_vect)) {
  message("\nReprojecting forecast raster to Somaliland CRS...")
  forecast_raster <- terra::project(forecast_raster, terra::crs(som_vect))
}

forecast_raster <- terra::crop(forecast_raster, som_vect)
forecast_raster <- terra::mask(forecast_raster, som_vect)
names(forecast_raster) <- "Rainfall_5Day_mm"

terra::writeRaster(
  forecast_raster,
  file.path(out_ras, "Somaliland_CHIRPS3_GEFS_5Day_Forecast.tif"),
  overwrite = TRUE
)

forecast_df <- as.data.frame(forecast_raster, xy = TRUE, na.rm = TRUE)
names(forecast_df)[3] <- "Rainfall_5Day_mm"

forecast_df$Rainfall_Category <- cut(
  forecast_df$Rainfall_5Day_mm,
  breaks = c(-Inf, 5, 20, 50, 100, Inf),
  labels = c("Very Low", "Low", "Moderate", "High", "Very High"),
  include.lowest = TRUE
)

rainfall_map <- ggplot2::ggplot() +
  geom_raster(data = forecast_df, aes(x = x, y = y, fill = Rainfall_5Day_mm)) +
  geom_sf(data = somaliland, fill = NA, colour = "white", linewidth = 0.8) +
  geom_sf_text(data = label_points, aes(label = Region), size = 2.8, fontface = "bold", colour = "#111111") +
  scale_fill_gradientn(
    colours = c("#FFF7BC", "#FEC44F", "#FE9929", "#EC7014", "#CC4C02", "#8C2D04", "#54278F"),
    name = "Rainfall\n(mm)",
    labels = scales::label_number(accuracy = 1)
  ) +
  labs(
    title = "5-DAY RAINFALL FORECAST -- SOMALILAND",
    subtitle = paste0("CHIRPS3-GEFS | Latest forecast: ", forecast_5$filename),
    caption = "Source: Climate Hazards Center, UC Santa Barbara | Bias-corrected GEFS precipitation forecast | 0.05\u00b0"
  ) +
  map_coord +
  gis_theme

print(rainfall_map)
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
save_map(rainfall_map, "04_5_Day_Rainfall_Forecast.png", width = 13, height = 9, dpi = 600)
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
# ---- regional table ---------------------------------------------------------

forecast_extract <- terra::extract(forecast_raster, som_vect, fun = mean, na.rm = TRUE, ID = TRUE)

if (nrow(forecast_extract) == 0) {
  stop("\nRegional rainfall extraction returned zero rows.")
}

value_col <- names(forecast_extract)[2]

regional_forecast <- data.frame(
  ID = forecast_extract$ID,
  Forecast_Rainfall_mm = as.numeric(forecast_extract[[value_col]])
) |>
  dplyr::filter(!is.na(Forecast_Rainfall_mm)) |>
  dplyr::mutate(
    Region = somaliland$Region[ID],
    Forecast_Rainfall_mm = round(Forecast_Rainfall_mm, 1),
    Rainfall_Category = cut(
      Forecast_Rainfall_mm,
      breaks = c(-Inf, 5, 20, 50, 100, Inf),
      labels = c("Very Low", "Low", "Moderate", "High", "Very High"),
      include.lowest = TRUE
    ),
    Heavy_Rainfall_Watch = Forecast_Rainfall_mm >= 50
  ) |>
  dplyr::select(Region, Forecast_Rainfall_mm, Rainfall_Category, Heavy_Rainfall_Watch) |>
  dplyr::arrange(dplyr::desc(Forecast_Rainfall_mm))

message("\n============================================================")
## 
## ============================================================
message("REGIONAL 5-DAY RAINFALL FORECAST")
## REGIONAL 5-DAY RAINFALL FORECAST
message("============================================================")
## ============================================================
print(regional_forecast)

write.csv(
  regional_forecast,
  file.path(out_tab, "Somaliland_5_Day_Rainfall_Forecast_By_Region.csv"),
  row.names = FALSE
)

# ============================================================
# FINAL SUMMARY
# ============================================================

cat("\n\n============================================================\n")
cat(" SOMALILAND GIS CLIMATE EARLY WARNING -- COMPLETED\n")
cat("============================================================\n")
cat("Sababi shapefile: ", basename(boundary_shp), "\n")
cat("Regions: ", nrow(somaliland), "\n")
cat("Latest 5-day forecast: ", forecast_5$filename, "\n")
cat("\nMAPS:\n")
cat("1. 01_Somaliland_Regions.png\n")
cat("2. 02_Average_Annual_Temperature.png\n")
cat("3. 03_Elevation.png\n")
cat("4. 04_5_Day_Rainfall_Forecast.png\n")
cat("\nTABLES:\n")
cat("Somaliland_Temperature_By_Region.csv\n")
cat("Somaliland_Elevation_By_Region.csv\n")
cat("Somaliland_5_Day_Rainfall_Forecast_By_Region.csv\n")
cat("\nOUTPUT FOLDER:\n", root_dir, "\n")
cat("============================================================\n")

message("\nALL FOUR MAPS HAVE BEEN CREATED SUCCESSFULLY.")
## 
## ALL FOUR MAPS HAVE BEEN CREATED SUCCESSFULLY.

6 4. Regional Summary Tables

6.1 Temperature by Region

##         Region Mean_Temperature_C
## 1        Sahil               26.5
## 2 Maroodi Jeex               22.9
## 3     Togdheer               23.8
## 4        Awdal               27.3
## 5         Sool               25.0
## 6       Sanaag               23.3

6.2 Elevation by Region

##         Region Mean_Elevation_m
## 1        Sahil            488.5
## 2 Maroodi Jeex           1147.4
## 3     Togdheer            938.2
## 4        Awdal            516.0
## 5         Sool            708.0
## 6       Sanaag            952.0

6.3 5-Day Rainfall Forecast by Region

##         Region Forecast_Rainfall_mm Rainfall_Category Heavy_Rainfall_Watch
## 1       Sanaag                 10.4               Low                FALSE
## 2 Maroodi Jeex                  6.7               Low                FALSE
## 3        Sahil                  5.9               Low                FALSE
## 4     Togdheer                  3.6          Very Low                FALSE
## 5        Awdal                  2.1          Very Low                FALSE
## 6         Sool                  1.4          Very Low                FALSE

7 5. Regional Graphs

These charts are added specifically for the RPubs report so that the publication contains both maps and graphs.

7.1 Average Annual Temperature by Region

7.2 Mean Elevation by Region

7.3 Next 5-Day Rainfall Forecast by Region

8 6. Rainfall Risk Categories

9 7. Interpretation Guide

9.0.1 Administrative map

Shows the spatial organization of the six Somaliland regions.

9.0.2 Temperature map

Shows the spatial pattern of average annual temperature derived from WorldClim minimum and maximum temperature layers.

9.0.3 Elevation map

Shows the spatial distribution of elevation using the DEM accessed through elevatr.

9.0.4 5-day rainfall map

Shows the latest available CHIRPS3-GEFS 5-day precipitation forecast clipped to Somaliland.

9.0.5 Regional graphs

The bar charts provide a direct region-to-region comparison of temperature, elevation, and forecast rainfall.

10 8. RPubs Publishing Steps

  1. Open this .Rmd file in RStudio.
  2. Confirm the SL_Six_Regions shapefile folder is in Downloads.
  3. Click Knit → Knit to HTML.
  4. Wait for the maps and charts to finish rendering.
  5. Check that all four maps and the regional charts appear in the HTML.
  6. Click Publish → Publish to RPubs.
  7. Sign in to RPubs and publish the report.

11 9. Expected RPubs Output

The published report should contain:

Map 1: Somaliland Six Regions
Map 2: Average Annual Temperature
Map 3: Elevation / DEM
Map 4: Next 5-Day Rainfall Forecast

and regional graphs for:

Temperature → Elevation → Rainfall Forecast → Rainfall Categories

12 Author

Yahye S. Rageh

Somaliland GIS Climate & 5-Day Rainfall Early Warning Analysis