This practical R Markdown presents a spatial climate analysis workflow for Somalia, Kenya, and Ethiopia.
The supplied script produces:
The source script specifies that maps should appear in the RStudio Plot Pane and should not be saved as PNG, PDF, or JPEG files.
geodataelevatrrequired_packages <- c(
"sf", "terra", "ggplot2", "dplyr",
"geodata", "elevatr", "viridis", "scales"
)
new_packages <- required_packages[
!(required_packages %in% installed.packages()[, "Package"])
]
if(length(new_packages) > 0){
install.packages(new_packages, dependencies = TRUE)
}
library(sf)
library(terra)
library(ggplot2)
library(dplyr)
library(geodata)
library(elevatr)
library(viridis)
library(scales)
options(stringsAsFactors = FALSE)
target_crs <- 4326The following code is the complete supplied workflow, organized as one executable R code section.
############################################################
# ==========================================================
# π EAST AFRICA SPATIAL CLIMATE ANALYSIS
# SOMALIA + KENYA + ETHIOPIA
# ==========================================================
#
# Maps:
# 1. Administrative Regions
# 2. Average Annual Temperature
# 3. Elevation / DEM
# 4. Annual Rainfall
# 5. Rainfall Risk Classification
#
# Countries:
# πΈπ΄ Somalia
# π°πͺ Kenya
# πͺπΉ Ethiopia
#
# Data Sources:
# - Somalia: User-provided shapefile
# - Kenya/Ethiopia: GADM via geodata
# - Temperature: WorldClim
# - Rainfall: WorldClim
# - Elevation: elevatr
#
# Author: Yahye S. Rageh
#
# IMPORTANT:
# ALL MAPS APPEAR IN RSTUDIO PLOT PANE
# NO ggsave()
# NO png()
# NO pdf()
# NO JPEG OUTPUT
#
############################################################
# ==========================================================
# 1. CLEAR ENVIRONMENT
# ==========================================================
rm(list = ls())
graphics.off()
options(stringsAsFactors = FALSE)
# ==========================================================
# 2. INSTALL / LOAD REQUIRED PACKAGES
# ==========================================================
required_packages <- c(
"sf",
"terra",
"ggplot2",
"dplyr",
"geodata",
"elevatr",
"viridis",
"scales"
)
new_packages <- required_packages[
!(required_packages %in%
installed.packages()[, "Package"])
]
if(length(new_packages) > 0){
install.packages(
new_packages,
dependencies = TRUE
)
}
library(sf)
library(terra)
library(ggplot2)
library(dplyr)
library(geodata)
library(elevatr)
library(viridis)
library(scales)
# ==========================================================
# 3. GENERAL SETTINGS
# ==========================================================
# WGS84
target_crs <- 4326
# ==========================================================
# 4. PROFESSIONAL GIS THEME
# ==========================================================
gis_theme <-
theme_minimal() +
theme(
plot.title =
element_text(
size = 20,
face = "bold",
colour = "#08306B",
hjust = 0.5
),
plot.subtitle =
element_text(
size = 11,
colour = "#4D4D4D",
hjust = 0.5
),
plot.caption =
element_text(
size = 8,
colour = "grey40",
hjust = 0
),
legend.position = "right",
legend.title =
element_text(
face = "bold",
size = 10
),
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 = "#F8FAFC",
colour = NA
),
plot.background =
element_rect(
fill = "white",
colour = NA
),
plot.margin =
margin(
10,
10,
10,
10
)
)
# ==========================================================
# 5. FUNCTION β PREPARE ADMINISTRATIVE DATA
# ==========================================================
prepare_admin <- function(
shp,
country_name
){
shp <-
st_make_valid(shp)
shp <-
st_transform(
shp,
target_crs
)
shp$Country <-
country_name
return(shp)
}
# ==========================================================
# 6. SOMALIA β LOAD YOUR SHAPEFILE
# ==========================================================
cat("\n")
cat("============================================================\n")
cat(" πΈπ΄ SOMALIA\n")
cat("============================================================\n")
cat(
"Reading Somalia administrative shapefile...\n"
)
som_adm1 <-
st_read(
"som_admin_boundaries.shp",
layer =
"som_admin1",
quiet =
TRUE
)
som_adm1 <-
prepare_admin(
som_adm1,
"Somalia"
)
cat(
"Somalia ADM1 regions:",
nrow(som_adm1),
"\n"
)
# ==========================================================
# 7. KENYA β AUTOMATIC ADMINISTRATIVE SHAPEFILE
# ==========================================================
cat("\n")
cat("============================================================\n")
cat(" π°πͺ KENYA\n")
cat("============================================================\n")
cat(
"Downloading/loading Kenya ADM1 boundaries...\n"
)
ken_adm1 <-
geodata::gadm(
country = "KEN",
level = 1,
path = tempdir()
)
ken_adm1 <-
st_as_sf(
ken_adm1
)
ken_adm1 <-
prepare_admin(
ken_adm1,
"Kenya"
)
cat(
"Kenya ADM1 regions:",
nrow(ken_adm1),
"\n"
)
# ==========================================================
# 8. ETHIOPIA β AUTOMATIC ADMINISTRATIVE SHAPEFILE
# ==========================================================
cat("\n")
cat("============================================================\n")
cat(" πͺπΉ ETHIOPIA\n")
cat("============================================================\n")
cat(
"Downloading/loading Ethiopia ADM1 boundaries...\n"
)
eth_adm1 <-
geodata::gadm(
country = "ETH",
level = 1,
path = tempdir()
)
eth_adm1 <-
st_as_sf(
eth_adm1
)
eth_adm1 <-
prepare_admin(
eth_adm1,
"Ethiopia"
)
cat(
"Ethiopia ADM1 regions:",
nrow(eth_adm1),
"\n"
)
# ==========================================================
# 9. COMBINE COUNTRIES
# ==========================================================
all_countries <-
list(
Somalia =
som_adm1,
Kenya =
ken_adm1,
Ethiopia =
eth_adm1
)
# ==========================================================
# 10. FUNCTION β ADMINISTRATIVE MAP
# ==========================================================
plot_admin_map <-
function(
admin,
country
){
p <-
ggplot(admin) +
geom_sf(
aes(
fill =
NAME_1
),
colour =
"white",
linewidth =
0.45
) +
geom_sf_text(
aes(
label =
NAME_1
),
size =
2.5,
fontface =
"bold",
colour =
"black"
) +
scale_fill_viridis_d(
option =
"turbo",
guide =
"none"
) +
labs(
title =
paste(
toupper(country),
"ADMINISTRATIVE REGIONS"
),
subtitle =
"Administrative Level 1 Boundaries",
caption =
"Source: Administrative Boundary Dataset"
) +
gis_theme
# Somalia uses adm1_name
if(country == "Somalia"){
p <-
ggplot(admin) +
geom_sf(
aes(
fill =
adm1_name
),
colour =
"white",
linewidth =
0.45
) +
geom_sf_text(
aes(
label =
adm1_name
),
size =
2.5,
fontface =
"bold"
) +
scale_fill_viridis_d(
option =
"turbo",
guide =
"none"
) +
labs(
title =
"SOMALIA ADMINISTRATIVE REGIONS",
subtitle =
"Administrative Level 1 Boundaries",
caption =
"Source: Somalia Administrative Boundary Dataset"
) +
gis_theme
}
return(p)
}
# ==========================================================
# 11. DISPLAY ADMINISTRATIVE MAPS
# ==========================================================
cat("\nDisplaying administrative maps...\n")
print(
plot_admin_map(
som_adm1,
"Somalia"
)
)
print(
plot_admin_map(
ken_adm1,
"Kenya"
)
)
print(
plot_admin_map(
eth_adm1,
"Ethiopia"
)
)
# ==========================================================
# 12. FUNCTION β WORLDCLIM TEMPERATURE
# ==========================================================
get_temperature <- function(
country_code,
country_name,
admin
){
cat(
"\nDownloading WorldClim Tmin:",
country_name,
"\n"
)
tmin <-
geodata::worldclim_country(
country =
country_name,
var =
"tmin",
path =
tempdir()
)
cat(
"Downloading WorldClim Tmax:",
country_name,
"\n"
)
tmax <-
geodata::worldclim_country(
country =
country_name,
var =
"tmax",
path =
tempdir()
)
# --------------------------------------------------------
# CONVERT 0.1 Β°C TO Β°C
# --------------------------------------------------------
tmin <-
tmin / 10
tmax <-
tmax / 10
# --------------------------------------------------------
# IMPORTANT:
# PIXEL-BY-PIXEL MEAN TEMPERATURE
# --------------------------------------------------------
mean_temp <-
(tmin + tmax) / 2
# --------------------------------------------------------
# COUNTRY SPATIAL MASK
# --------------------------------------------------------
admin_vect <-
terra::vect(admin)
mean_temp <-
terra::crop(
mean_temp,
admin_vect
)
mean_temp <-
terra::mask(
mean_temp,
admin_vect
)
return(mean_temp)
}
# ==========================================================
# 13. FUNCTION β ANNUAL RAINFALL
# ==========================================================
get_rainfall <- function(
country_name,
admin
){
cat(
"\nDownloading WorldClim precipitation:",
country_name,
"\n"
)
precipitation <-
geodata::worldclim_country(
country =
country_name,
var =
"prec",
path =
tempdir()
)
# --------------------------------------------------------
# SUM 12 MONTHS = ANNUAL RAINFALL
# --------------------------------------------------------
annual_rainfall <-
sum(
precipitation
)
admin_vect <-
terra::vect(admin)
annual_rainfall <-
terra::crop(
annual_rainfall,
admin_vect
)
annual_rainfall <-
terra::mask(
annual_rainfall,
admin_vect
)
return(
annual_rainfall
)
}
# ==========================================================
# 14. FUNCTION β ELEVATION
# ==========================================================
get_elevation <- function(
admin,
country_name
){
cat(
"\nDownloading elevation:",
country_name,
"\n"
)
# Lower zoom keeps the analysis practical
elev <-
elevatr::get_elev_raster(
locations =
admin,
z =
6,
clip =
"locations"
)
elev <-
terra::rast(
elev
)
admin_vect <-
terra::vect(admin)
elev <-
terra::crop(
elev,
admin_vect
)
elev <-
terra::mask(
elev,
admin_vect
)
return(elev)
}
# ==========================================================
# 15. FUNCTION β TEMPERATURE MAP
# ==========================================================
plot_temperature_map <-
function(
temperature,
admin,
country
){
temp_df <-
as.data.frame(
temperature,
xy =
TRUE,
na.rm =
TRUE
)
colnames(temp_df)[3] <-
"Temperature"
# Somalia boundary label
if(country == "Somalia"){
label_field <-
"adm1_name"
} else {
label_field <-
"NAME_1"
}
p <-
ggplot() +
geom_raster(
data =
temp_df,
aes(
x =
x,
y =
y,
fill =
Temperature
)
) +
geom_sf(
data =
admin,
fill =
NA,
colour =
"white",
linewidth =
0.45
) +
geom_sf_text(
data =
admin,
aes_string(
label =
label_field
),
size =
2.4,
fontface =
"bold"
) +
scale_fill_gradientn(
colours =
c(
"#313695",
"#4575B4",
"#74ADD1",
"#ABD9E9",
"#FFFFBF",
"#FDAE61",
"#F46D43",
"#D73027"
),
name =
"Temperature\n(Β°C)",
labels =
function(x){
paste0(
round(
x,
1
),
"Β°C"
)
}
) +
labs(
title =
paste(
toupper(country),
"AVERAGE ANNUAL TEMPERATURE"
),
subtitle =
"Mean Temperature Derived from WorldClim Tmin & Tmax",
caption =
"Source: WorldClim v2"
) +
gis_theme
return(p)
}
# ==========================================================
# 16. FUNCTION β ELEVATION MAP
# ==========================================================
plot_elevation_map <-
function(
elevation,
admin,
country
){
elev_df <-
as.data.frame(
elevation,
xy =
TRUE,
na.rm =
TRUE
)
colnames(elev_df)[3] <-
"Elevation"
if(country == "Somalia"){
label_field <-
"adm1_name"
} else {
label_field <-
"NAME_1"
}
p <-
ggplot() +
geom_raster(
data =
elev_df,
aes(
x =
x,
y =
y,
fill =
Elevation
)
) +
geom_sf(
data =
admin,
fill =
NA,
colour =
"white",
linewidth =
0.45
) +
geom_sf_text(
data =
admin,
aes_string(
label =
label_field
),
size =
2.4,
fontface =
"bold"
) +
scale_fill_gradientn(
colours =
c(
"#0B3C5D",
"#328CC1",
"#99C24D",
"#E6AF2E",
"#D95D39",
"#7F2704"
),
name =
"Elevation\n(m)"
) +
labs(
title =
paste(
toupper(country),
"DIGITAL ELEVATION MODEL"
),
subtitle =
"Spatial Elevation Distribution",
caption =
"Source: Elevation data via elevatr"
) +
gis_theme
return(p)
}
# ==========================================================
# 17. FUNCTION β RAINFALL MAP
# ==========================================================
plot_rainfall_map <-
function(
rainfall,
admin,
country
){
rain_df <-
as.data.frame(
rainfall,
xy =
TRUE,
na.rm =
TRUE
)
colnames(rain_df)[3] <-
"Rainfall"
if(country == "Somalia"){
label_field <-
"adm1_name"
} else {
label_field <-
"NAME_1"
}
p <-
ggplot() +
geom_raster(
data =
rain_df,
aes(
x =
x,
y =
y,
fill =
Rainfall
)
) +
geom_sf(
data =
admin,
fill =
NA,
colour =
"white",
linewidth =
0.45
) +
geom_sf_text(
data =
admin,
aes_string(
label =
label_field
),
size =
2.4,
fontface =
"bold"
) +
scale_fill_gradientn(
colours =
c(
"#F7FCFD",
"#E0ECF4",
"#BFD3E6",
"#9EBCDA",
"#8C96C6",
"#8856A7",
"#810F7C"
),
name =
"Rainfall\n(mm/year)"
) +
labs(
title =
paste(
toupper(country),
"ANNUAL RAINFALL DISTRIBUTION"
),
subtitle =
"WorldClim Long-Term Precipitation Baseline",
caption =
"Source: WorldClim v2 Precipitation Dataset"
) +
gis_theme
return(p)
}
# ==========================================================
# 18. FUNCTION β RAINFALL RISK MAP
# ==========================================================
plot_risk_map <-
function(
rainfall,
admin,
country
){
rain_df <-
as.data.frame(
rainfall,
xy =
TRUE,
na.rm =
TRUE
)
colnames(rain_df)[3] <-
"Rainfall"
# --------------------------------------------------------
# RAINFALL CLASSES
# --------------------------------------------------------
rain_df <-
rain_df %>%
mutate(
Rainfall_Class =
cut(
Rainfall,
breaks =
c(
-Inf,
200,
400,
600,
800,
Inf
),
labels =
c(
"Very Dry",
"Dry",
"Moderate",
"Wet",
"Very Wet"
),
include.lowest =
TRUE
)
)
if(country == "Somalia"){
label_field <-
"adm1_name"
} else {
label_field <-
"NAME_1"
}
p <-
ggplot() +
geom_raster(
data =
rain_df,
aes(
x =
x,
y =
y,
fill =
Rainfall_Class
)
) +
geom_sf(
data =
admin,
fill =
NA,
colour =
"white",
linewidth =
0.5
) +
geom_sf_text(
data =
admin,
aes_string(
label =
label_field
),
size =
2.4,
fontface =
"bold"
) +
scale_fill_manual(
values =
c(
"Very Dry" =
"#D73027",
"Dry" =
"#FC8D59",
"Moderate" =
"#FEE08B",
"Wet" =
"#91CF60",
"Very Wet" =
"#1A9850"
),
name =
"Rainfall\nRisk"
) +
labs(
title =
paste(
toupper(country),
"RAINFALL RISK CLASSIFICATION"
),
subtitle =
"Based on WorldClim Long-Term Annual Rainfall",
caption =
"Very Dry β Dry β Moderate β Wet β Very Wet"
) +
gis_theme
return(p)
}
# ==========================================================
# 19. ANALYSE SOMALIA
# ==========================================================
cat("\n")
cat("============================================================\n")
cat(" πΈπ΄ ANALYSING SOMALIA\n")
cat("============================================================\n")
som_temperature <-
get_temperature(
country_code =
"SOM",
country_name =
"Somalia",
admin =
som_adm1
)
som_rainfall <-
get_rainfall(
country_name =
"Somalia",
admin =
som_adm1
)
som_elevation <-
get_elevation(
admin =
som_adm1,
country_name =
"Somalia"
)
# ==========================================================
# 20. DISPLAY SOMALIA MAPS
# ==========================================================
print(
plot_temperature_map(
som_temperature,
som_adm1,
"Somalia"
)
)
print(
plot_elevation_map(
som_elevation,
som_adm1,
"Somalia"
)
)
print(
plot_rainfall_map(
som_rainfall,
som_adm1,
"Somalia"
)
)
print(
plot_risk_map(
som_rainfall,
som_adm1,
"Somalia"
)
)
# ==========================================================
# 21. ANALYSE KENYA
# ==========================================================
cat("\n")
cat("============================================================\n")
cat(" π°πͺ ANALYSING KENYA\n")
cat("============================================================\n")
ken_temperature <-
get_temperature(
country_code =
"KEN",
country_name =
"Kenya",
admin =
ken_adm1
)
ken_rainfall <-
get_rainfall(
country_name =
"Kenya",
admin =
ken_adm1
)
ken_elevation <-
get_elevation(
admin =
ken_adm1,
country_name =
"Kenya"
)
# ==========================================================
# 22. DISPLAY KENYA MAPS
# ==========================================================
print(
plot_temperature_map(
ken_temperature,
ken_adm1,
"Kenya"
)
)
print(
plot_elevation_map(
ken_elevation,
ken_adm1,
"Kenya"
)
)
print(
plot_rainfall_map(
ken_rainfall,
ken_adm1,
"Kenya"
)
)
print(
plot_risk_map(
ken_rainfall,
ken_adm1,
"Kenya"
)
)
# ==========================================================
# 23. ANALYSE ETHIOPIA
# ==========================================================
cat("\n")
cat("============================================================\n")
cat(" πͺπΉ ANALYSING ETHIOPIA\n")
cat("============================================================\n")
eth_temperature <-
get_temperature(
country_code =
"ETH",
country_name =
"Ethiopia",
admin =
eth_adm1
)
eth_rainfall <-
get_rainfall(
country_name =
"Ethiopia",
admin =
eth_adm1
)
eth_elevation <-
get_elevation(
admin =
eth_adm1,
country_name =
"Ethiopia"
)
# ==========================================================
# 24. DISPLAY ETHIOPIA MAPS
# ==========================================================
print(
plot_temperature_map(
eth_temperature,
eth_adm1,
"Ethiopia"
)
)
print(
plot_elevation_map(
eth_elevation,
eth_adm1,
"Ethiopia"
)
)
print(
plot_rainfall_map(
eth_rainfall,
eth_adm1,
"Ethiopia"
)
)
print(
plot_risk_map(
eth_rainfall,
eth_adm1,
"Ethiopia"
)
)
# ==========================================================
# 25. FINAL MESSAGE
# ==========================================================
cat("\n")
cat("============================================================\n")
cat(" β
EAST AFRICA SPATIAL CLIMATE ANALYSIS COMPLETED\n")
cat("============================================================\n")
cat("\n")
cat("Countries analysed:\n")
cat("πΈπ΄ Somalia\n")
cat("π°πͺ Kenya\n")
cat("πͺπΉ Ethiopia\n")
cat("\n")
cat("Maps displayed for each country:\n")
cat("1. Administrative Regions\n")
cat("2. Average Annual Temperature\n")
cat("3. Digital Elevation Model\n")
cat("4. Annual Rainfall\n")
cat("5. Rainfall Risk Classification\n")
cat("\n")
cat("Data sources:\n")
cat("- Administrative boundaries\n")
cat("- WorldClim temperature\n")
cat("- WorldClim precipitation\n")
cat("- Elevation DEM\n")
cat("\n")
cat("IMPORTANT:\n")
cat("All maps appear in the RStudio Plot Pane.\n")
cat("No maps are saved to the computer.\n")
cat("No ggsave(), PNG or PDF output is used.\n")
cat("\n")
cat("============================================================\n")For each country, the workflow is designed to display:
The supplied script uses:
target_crs <- 4326for WGS84 mapping, prepares administrative data with
st_make_valid() and st_transform(), and uses
country-specific administrative boundaries. The Somalia shapefile is
expected to contain the som_admin1 layer in
som_admin_boundaries.shp.
All maps are displayed in the RStudio Plot Pane according to the
supplied script. No ggsave(), PNG, PDF, or JPEG export is
used.
.Rmd file in RStudio.