Introduction
Welcome to a new post in my ggplot Map Series! So far, we’ve covered polygon maps, vector maps using shapefiles, and Google Maps. Today, we’re going to explore yet another type of maps: raster maps, which are maps created from raster data. Specifically, we’ll be working with GeoTIFF files. They are special TIFF files embedded with georeferencing information such as map projection, coordinate systems, ellipsoids, and datums. GeoTIFFs are commonly used in satellite imagery, aerial photography, and digitized maps for GIS analysis (e.g., climatic maps). So without further ado, let’s get started!
Spatial raster maps in ggplot
1. Data preparation
The data we’ll be using are the WorldClim 1970-2000 5-km2 global monthly precipitation data, which contain 12 GeoTIFF files for each month. These GeoTIFF files can be read into R via the function rast() from the package terra as “SpatRaster” objects.
library(tidyverse)
library(terra)
### Download the WorldClim data
download.file("https://geodata.ucdavis.edu/climate/worldclim/2_1/base/wc2.1_2.5m_prec.zip",
destfile = "./WorldClim_Precipitation.zip")
unzip("WorldClim_Precipitation.zip", exdir = "WorldClim")
### Read the GeoTIFF files
worldclim_filepath <- list.files("./WorldClim", pattern = ".tif") %>%
str_c("./WorldClim/", .)
worldclim_precipitation_raster <- map(worldclim_filepath, rast) # rast() reads the GeoTIFF files as "SpatRaster" objects2. Create a basic raster map
The object “worldclim_precipitation_raster” contains 12 lists of SpatRaster objects, each of which has x- and y-coordinates as well as a layer for the precipitation. Sometimes, a SpatRaster object may have more than one layer, and we can get the layer names using names() from the package tidyterra.
Let’s make a map of the average global precipitation in January (the first list in “worldclim_precipitation_raster”). To create a map from a SpatRaster object, we can use the function geom_spatraster() from the package tidyterra. This function is a wrapper for geom_raster(). In this example, we only have one layer, so it’s not necessary to specify which layer we would like to map to the “fill” aesthetic (but we need to do so if there are multiple layers). We’ll also draw the country boundaries retrieved from the package rnaturalearth on the map.
### Get country boundaries
library(rnaturalearth)
country_boundary <- ne_countries(type = "countries", scale = "small")
### Create a palette for the precipitation
library(colorspace)
col_pal <- sequential_hcl(n = 7, palette = "BluYl", rev = T)
### Get raster layer name(s)
library(tidyterra)
names(worldclim_precipitation_raster[[1]])[1] "wc2.1_2.5m_prec_01"### Create a map
ggplot() +
geom_spatraster(data = worldclim_precipitation_raster[[1]]) + # January precipitation
geom_sf(data = country_boundary, fill = NA, color = "black", size = 0.8) + # country boundaries
scale_fill_gradientn(colors = col_pal[c(1, 4, 6, 7)],
na.value = "transparent",
name = "January precipitation (mm)") +
coord_sf(xlim = c(-180, 180), ylim = c(-90, 90), expand = F) +
theme_bw(base_size = 13) +
theme(legend.position = "top",
legend.direction = "horizontal",
legend.key.width = unit(0.8, "in"),
legend.title = element_text(hjust = -5, vjust = 0.75, size = 12),
legend.text = element_text(size = 10),
axis.text = element_text(color = "black"))
Change the map projection
To change the map projection, we need to transform the CRS of the precipitation SpatRaster object using the function project() from the package terra as well as the CRS of the country boundary sf object using the function st_transform() from the package sf. In this example, we’ll use the Robinson projection (ESRI:54030). Moreover, because of the new projection, the lon-lat lines and labels get a little messy and hard to control via ggplot functions. To deal with this, we can manually add the grid lines and labels to the map as additional geom layers. Basically, we create a set of lon-lat coordinates where we would like the grid lines to be drawn as well as a dataframe specifying the axis labels and positions. We can then convert these into sf objects and add them to the map using geom_sf(). (This method is inspired by this stack overflow solution!)
library(sf)
### Transform the map projection
terraOptions(progress = 0) # hide the progress bar
january_precipitation_raster_robinson <- project(worldclim_precipitation_raster[[1]], "ESRI:54030") # ESRI:54030 for the Robinson projection
country_boundary_robinson <- st_transform(country_boundary, crs = "ESRI:54030")
### Create sf objects for the long-lat lines and latitude labels
# longitude lines
lon_lines <- lapply(c(-180, -120, -60, 0, 60, 120, 180), function(x) {cbind(x, seq(-90, 90, 1))})
# latitude lines
lat_lines <- lapply(c(-90, -60, -30, 0, 30, 60, 90), function(x) {cbind(seq(-180, 180, 1), x)})
# sf object for the long-lat lines
grid_lines <- st_sf(geometry = st_sfc(st_multilinestring(x = append(lon_lines, lat_lines)), crs = "WGS84"))
# sf objects for the latitude labels
lat <- c(-90, -60, -30, 0, 30, 60, 90)
lat_labels <- data.frame(label = str_glue("{abs(lat)} °{if_else(lat == 0, '', if_else(lat < 0, 'S', 'N'))}"),
lon = -180,
lat = lat) %>%
st_as_sf(., coords = c("lon", "lat"), crs = "WGS84")
### Make the map
ggplot() +
geom_spatraster(data = january_precipitation_raster_robinson) +
geom_sf(data = country_boundary_robinson, fill = NA, color = "black", size = 0.8) +
geom_sf(data = grid_lines, linewidth = 0.1, color = "grey") + # add the grid lines
geom_sf_text(data = lat_labels, aes(label = label), size = 3, # add the latitude labels
nudge_x = c(-2e6, -1.5e6, -1.3e6, -1e6, -1.3e6, -1.5e6, -2e6)) + # need to adjust the amount of nudging depending on the specific map projection used
scale_fill_gradientn(colors = col_pal[c(1, 4, 6, 7)],
na.value = "transparent",
name = "January precipitation (mm)") +
theme_minimal() +
coord_sf(expand = F,
crs = "ESRI:54030",
clip = "off") +
theme(legend.position = "top",
legend.direction = "horizontal",
legend.key.width = unit(0.8, "in"),
legend.title = element_text(hjust = -5, vjust = 0.75, size = 12),
legend.text = element_text(size = 10),
panel.border = element_blank(),
plot.background = element_blank(),
panel.grid = element_blank(),
axis.title = element_blank(),
axis.text = element_text(color = "black"))
Create a choropleth map of average precipitation by country
The package terra has a suite of functions for geospatial computation. A useful one is extract(), which extracts cell values from SpatRaster layers based on a set of vector locations (points, lines, or polygons) and summarizes the extracted values (min, max, mean, length, etc.).
Below, we’ll extract the cell values from the January precipitation SpatRaster based on the country boundary polygons and calculate the average precipitation for each country. We’ll then join the country averages to the country boundary sf object to make a choropleth map.
### Summarize the average precipitation by country
january_precipitation_country_mean <- extract(worldclim_precipitation_raster[[1]], # the SpatRaster to extract cell values from
country_boundary, # the vector polygons to define the zones
fun = mean, # the summary function
na.rm = T) %>%
rename(mean_precipitation = `wc2.1_2.5m_prec_01`)
### Add country averages to the country boundary sf object
january_precipitation_country_mean_sf <- bind_cols(country_boundary_robinson, january_precipitation_country_mean)
### Make a choropleth map
ggplot() +
geom_sf(data = january_precipitation_country_mean_sf, aes(fill = mean_precipitation), color = "black", size = 0.8) +
geom_sf(data = grid_lines, linewidth = 0.1, color = "grey") +
geom_sf_text(data = lat_labels, aes(label = label), size = 3,
nudge_x = c(-2e6, -1.5e6, -1.3e6, -1e6, -1.3e6, -1.5e6, -2e6)) +
scale_fill_gradientn(colors = col_pal[c(1, 4, 6, 7)],
na.value = "transparent",
name = "January precipitation (mm)") +
theme_minimal() +
coord_sf(expand = F,
crs = "ESRI:54030",
clip = "off") +
theme(legend.position = "top",
legend.direction = "horizontal",
legend.key.width = unit(0.8, "in"),
legend.title = element_text(hjust = -5, vjust = 0.75, size = 12),
legend.text = element_text(size = 10),
panel.border = element_blank(),
plot.background = element_blank(),
panel.grid = element_blank(),
axis.title = element_blank(),
axis.text = element_text(color = "black"))
Create an animation of monthly global precipitation
Remember we have 12 SpatRaster objects in “worldclim_precipitation_raster”. We can create a map for each and turn them into an animation to explore the temporal patterns of world precipitation. We’ll simply loop through each month and make a map for it, just like what we did for January. This gives us a list of 12 maps. We can then use the functions from the package
magick to make an animation of those maps.
### Create a map for each month
p_world_precipitation_by_month <- map(1:12, function(x) {
# SpatRaster for each month
monthly_precipitation_raster <- worldclim_precipitation_raster[[x]]
# name of precipitation layer
monthly_precipitation_raster_layer <- names(monthly_precipitation_raster)
# set the max precipitation to 1,000 mm (for better visualization)
monthly_precipitation_raster <- monthly_precipitation_raster %>%
mutate(!!monthly_precipitation_raster_layer := if_else(!!sym(monthly_precipitation_raster_layer) >= 1000,
1000,
!!sym(monthly_precipitation_raster_layer)))
# change the map projection
terraOptions(progress = 0)
monthly_precipitation_raster_robinson <- project(monthly_precipitation_raster, "ESRI:54030")
# name of month
month_name <- month.name[x]
# map
p <- ggplot() +
geom_spatraster(data = monthly_precipitation_raster_robinson) +
geom_sf(data = country_boundary_robinson, fill = NA, color = "black", size = 0.8) +
geom_sf(data = grid_lines, linewidth = 0.1, color = "grey") +
geom_sf_text(data = lat_labels, aes(label = label), size = 3,
nudge_x = c(-2e6, -1.5e6, -1.3e6, -1e6, -1.3e6, -1.5e6, -2e6)) +
scale_fill_gradientn(colors = c("#ffffd9", "#238443", "#41b6c4", "#225ea8", "#253494", "#54278f", "#3f007d", "#980043", "#980043"),
limits = c(0, 1000),
breaks = c(0, 250, 500, 750, 1000),
labels = c("0 mm", "250 mm", "500 mm", "750 mm", "> 1,000 mm"),
na.value = "transparent",
name = NULL) +
theme_minimal() +
labs(title = month_name) +
coord_sf(expand = F,
crs = "ESRI:54030",
clip = "off") +
theme(legend.position = "top",
legend.direction = "horizontal",
legend.key.width = unit(0.8, "in"),
legend.title = element_text(hjust = -5, vjust = 0.75, size = 12),
legend.text = element_text(size = 10),
panel.border = element_blank(),
plot.background = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(hjust = 0.51, size = 13, vjust = 0),
axis.title = element_blank(),
axis.text = element_text(color = "black"))
return(p)
})
### Make an animation
library(magick)
# create a Magick image object
image <- image_graph(width = 1350, height = 900, res = 200)
# display the maps on the plotting device
out <- lapply(1:12, function(x) {
p <- p_world_precipitation_by_month[[x]]
print(p)
})
# close the plotting device
dev.off()png
2 # animate the images
animation <- image_animate(image, fps = 2, optimize = T)
animation
Summary
To recap, we learned how to create raster maps from GeoTIFF files in ggplot. We first created a basic raster map of world precipitation using the packages terra and tidyterra. Next, we modified the map by changing the map projection and customizing the grid lines and axis labels. We also did some basic manipulations of the raster data by calculating the average precipitation in each country and created a choropleth map for it. Finally, we made a world precipitation map for each of the 12 months and turned them into an animation using the package magick. Raster data are quite common, especially when it comes to geospatial analysis. Having some “raster” skills under your belt will certainly go a long way!
Hope you learn something useful from this post and don’t forget to leave your comments and suggestions below if you have any!
