Module III: Working with Raster Data in R
using the terra Package
Abdisalam Hassan Muse, PhD
2025-07-26
1. Introduction to the terra Package
The terra package is the modern, high-performance
replacement for the raster package. It’s designed to work
with both raster and vector data
efficiently. A raster is a grid of cells, where each cell holds a value
(e.g., elevation, temperature, or rainfall).
In terra, a raster object is called a
SpatRaster.
2. Reading and Inspecting Raster Data
We can download climate data directly using the geodata
package, which works seamlessly with terra. Let’s download
the average monthly precipitation for Somalia from the WorldClim
database.
library(terra)
library(geodata)
# Download average monthly precipitation data for Somalia
# This returns a SpatRaster with 12 layers (one for each month)
som_precip <- worldclim_country(country = "Somalia", var = "prec", path = tempdir())
# Inspect the SpatRaster object
som_precip## class : SpatRaster
## dimensions : 1680, 1320, 12 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 40.5, 51.5, -2, 12 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : SOM_wc2.1_30s_prec.tif
## names : SOM_w~rec_1, SOM_w~rec_2, SOM_w~rec_3, SOM_w~rec_4, SOM_w~rec_5, SOM_w~rec_6, ...
## min values : 0, 0, 0, 3, 0, 0, ...
## max values : 43, 49, 87, 208, 285, 165, ...The output tells us: - class: SpatRaster - dimensions: The number of rows, columns, and layers (nlyr = 12 for 12 months). - resolution: The size of each grid cell in degrees. - extent: The geographic bounding box. - crs: The Coordinate Reference System.
We can plot the raster. By default, terra plots all 12
layers.
3. Raster Operations: Summarizing, Cropping, and Masking
A common workflow involves processing a raster to fit our specific research needs.
A. Summarizing Layers We have 12 monthly layers, but
we might want the annual average precipitation. We can use the
mean() function to calculate the mean value for each cell
across all 12 layers.
# Calculate the mean precipitation across all 12 layers
annual_avg_precip <- mean(som_precip)
# Plot the single-layer result
plot(annual_avg_precip, main = "Average Annual Precipitation in Somalia (mm)")
B. Cropping and Masking The raster we downloaded is
a simple rectangle. We often want to clip it to the exact shape of our
study area (e.g., the administrative boundary of Somalia). This is a
two-step process: 1. crop(): Reduces the
raster to the bounding box of the shape, making it smaller and faster to
process. 2. mask(): Sets all raster cells
that fall outside the shape to NA, giving it the
precise outline of our study area.
# We need the vector boundary of Somalia from Lecture 1
library(sf)
som_adm0 <- st_read("som_admbnda_adm0_ocha_20250108.shp", quiet = TRUE)
# Convert the sf object to a SpatVector for use with terra
som_adm0_vect <- vect(som_adm0)
# 1. Crop the raster to the extent of the Somalia boundary
precip_cropped <- crop(annual_avg_precip, som_adm0_vect)
# 2. Mask the cropped raster with the Somalia polygon
precip_masked <- mask(precip_cropped, som_adm0_vect)
# Plot the final, masked raster
plot(precip_masked, main = "Average Annual Precipitation Masked to Somalia's Boundary")
4. Extracting Raster Values
The real power comes from combining raster and vector data. The
extract() function allows us to “pull” values from a raster
at specific locations (points) or summarize them over areas
(polygons).
Example 1: Extracting Precipitation at Specific Cities Let’s find the average annual precipitation for Mogadishu and Hargeisa.
# 1. Create a SpatVector of our cities
cities_df <- data.frame(
name = c("Mogadishu", "Hargeisa"),
lon = c(45.333, 44.067),
lat = c(2.033, 9.567)
)
cities_vect <- vect(cities_df, geom=c("lon", "lat"), crs="EPSG:4326")
# 2. Extract the raster values at these point locations
city_precip_values <- extract(precip_masked, cities_vect)
# 3. Combine the results with the city names
city_results <- cbind(cities_df, city_precip_values)
print(city_results)## name lon lat ID mean
## 1 Mogadishu 45.333 2.033 1 36.66667
## 2 Hargeisa 44.067 9.567 2 38.66667The result shows Hargeisa receives significantly more rainfall than Mogadishu, which matches our understanding of the local climate.
Example 2: Summarizing Precipitation by Region What
is the average precipitation for each administrative region in Somalia?
We can use extract() with polygons.
library(dplyr)
library(ggplot2)
library(tidyterra) # For plotting terra objects with ggplot
# 1. Load the regions vector data (from Lecture 1)
som_adm1 <- st_read("som_admbnda_adm1_ocha_20250108.shp", quiet = TRUE)
som_adm1_vect <- vect(som_adm1)
# 2. Extract the mean precipitation for each polygon (region)
# We use fun="mean" to specify the summary statistic. na.rm=TRUE ignores NA cells.
regional_precip <- extract(precip_masked, som_adm1_vect, fun = "mean", na.rm = TRUE)
# 3. Add the extracted values back to our vector object
som_adm1_vect$avg_precip <- regional_precip[,2] # The second column has the values
# 4. Map the results
ggplot() +
geom_spatvector(data = som_adm1_vect, aes(fill = avg_precip)) +
scale_fill_viridis_c(name = "Avg. Monthly\nPrecipitation (mm)") +
ggtitle("Average Monthly Precipitation by Region") +
theme_bw()
This map provides a clear and powerful visualization of the spatial distribution of rainfall across Somalia, revealing distinct climatic zones that profoundly influence livelihoods, food security, and potential for conflict. The patterns observed align closely with the country’s known topography and agricultural systems.
The “Sorghum Belt” and Riverine Areas (Highest Precipitation): The two areas with the highest average precipitation (shown in yellow, ~40 mm/month or ~480 mm/year) are the northwestern regions of Awdal and Woqooyi Galbeed (Somaliland) and the southwestern regions of Gedo, Bay, and the Jubas.
The northwestern highlands benefit from higher elevation and orographic lift, making them the most productive zone in Somaliland for rain-fed agriculture and pastoralism.
The southwestern area is the breadbasket of Somalia, benefiting not only from higher rainfall but also from the presence of the Jubba and Shabelle rivers, which support irrigation-based farming.
The Arid Central and Northeastern Regions (Lowest Precipitation): A vast corridor stretching from the central regions (Galgaduud, Mudug) up through the entire northeastern tip of the Horn (Nugaal and Bari regions of Puntland) experiences hyper-arid conditions (dark purple, ~10 mm/month or ~120 mm/year). This area is predominantly suited for nomadic pastoralism, and its inhabitants are highly vulnerable to drought. The low rainfall here makes water a scarce and critical resource.
Transitional Zones: The regions shown in green and teal represent transitional semi-arid zones. These areas, such as Togdheer in the north and Hiraan in the south, often serve as crucial grazing lands but are highly susceptible to climate variability.
In conclusion, this map is not just a climatic summary; it is a foundational map for understanding human geography in Somalia. It helps explain why agricultural activity is concentrated in the southwest, why pastoralist migration routes dominate the central and northeast, and which regions are most ecologically fragile and prone to drought-related humanitarian crises. Any analysis of food security, resource management, or development planning in Somalia must be grounded in this spatial reality of water availability.
Example 3: Summarizing Precipitation by District
What is the average precipitation for each administrative district in
Somalia? We can use extract() with polygons.
library(dplyr)
library(ggplot2)
library(tidyterra) # For plotting terra objects with ggplot
# 1. Load the districts vector data (from Lecture 1)
som_adm2 <- st_read("som_admbnda_adm2_ocha_20250108.shp", quiet = TRUE)
som_adm2_vect <- vect(som_adm2)
# 2. Extract the mean precipitation for each polygon (district)
# We use fun="mean" to specify the summary statistic. na.rm=TRUE ignores NA cells.
district_precip <- extract(precip_masked, som_adm2_vect, fun = "mean", na.rm = TRUE)
# 3. Add the extracted values back to our vector object
som_adm2_vect$avg_precip <- district_precip[,2] # The second column has the values
# 4. Map the results
ggplot() +
geom_spatvector(data = som_adm2_vect, aes(fill = avg_precip)) +
scale_fill_viridis_c(name = "Avg. Monthly\nPrecipitation (mm)") +
ggtitle("Average Monthly Precipitation by District") +
theme_bw()
```