Benjamin Bell: Blog

Today I presented at the British Ecological Society Palaeoecology Group "Masterclass in R" seminar series on using R for mapping and spatial analysis.

If you couldn't make the presentation and would like to see what I spoke about, as well as access the full R code, read more for details!

Using R for mapping and spatial analysis

Below you can find the presentation and the full R code that was used to create the maps and figures in the presentation. Please note, this is provided as code only, rather than the usual guides that I write. However, I have added lots of #comments to the code so that it should be easy to follow along, but please feel free to get in touch if you have any queries.

You can also check out some of my existing R guides for creating maps - although these use the older sp, and raster packages, the concepts are generally the same. I shall be updating these guides to use sf and terra soon.

The presentation is available on Google Drive, and should hopefully appear below! Please feel free to share.

The code is split into two parts - the first is for creating maps from existing data, and the second is for creating raster data and simple spatial analysis. The code has #comments to indicate the relevant slide.

This code was tested with R version 4.2.1 on Windows 11, and the latest version of terra and sf packages at the time this guide was published. Compatibility with other systems or versions of R or the packages cannot be guaranteed.

It is recommended that your system has at least 16 GB of RAM to follow this code.

Part 1: Creating and plotting maps

Part 2: Creating your own raster data and spatial analysis

Part 1: Creating and plotting maps - Full R code

####################################
### BES R Masterclass 2022
# Part 1: Creating and plotting maps - Full R code (with notes!)
# This code follows the presentation, available on my blog.
# You may find additional help/info for creating maps and working with spatial data on my blog.
# Please feel free to get in touch with any queries (contact details on my blog)
# Benjamin Bell: https://www.benjaminbell.co.uk
####################################

####################################
# Recommended 16 GB ram minimum to run this code!
# Tested with R v4.2.1 on Windows 11 
####################################

#### Slide 3 ####

### Install required packages
install.packages(c("sf", "terra"))

### Load packages (You must do this each time you start R!)
library(sf)
library(terra)

############
# Download GEBCO bathymetry data (Global data) from: https://www.gebco.net/data_and_products/gridded_bathymetry_data/
# Download: GEBCO_2022 Grid (ice surface elevation) - netCDF file
# (Put file in your dem folder)

# Important! If your system does not have enough RAM - use the downloader application to define a smaller area to download: https://download.gebco.net/

############

############
# Download vector data from Natural Earth: https://www.naturalearthdata.com/downloads/
# (Put files in your shapefiles folder)
############

#### Slide 6 ####


### Setup data file paths

### Important - This is just an example. You must change folder locations for your system, depending on where you put the downloaded files, otherwise the code below will not work.
### Windows users should use "/" instead of "\" for file/folder locations

# Location of DEM/bathymetric data (change for your system)
demf <- "~/data/dems/"
# Location of vector data (change for your system)
shpf <- "~/data/shapefiles/"

############
### Large regional Map
############

### Load bathymetric data as a raster object
### Important - You may need to edit this for your system/file location. For example, if you put the downloaded data in a sub-folder within the "dems" direction, you might need to do something like this: bath <- rast(paste0(demf, "sub-folder-name/GEBCO_2022.nc"))
bath <- rast(paste0(demf, "GEBCO_2022.nc"))

### Crop to your area (otherwise its too big!)
# Set an "extent" - this is your area, defined using coordinates (min/max lon, min/max lat)
e_r <- ext(c(-14, 0, 28, 38))
# Crop
bath_r <- crop(bath, e_r)

### Hillshade
# Create a hill shade to give your base map some texture
# Need to create a slope and aspect layer - then use those to create the hillshade
bath_r_slope <- terrain(bath_r, "slope", unit="radians")
bath_r_aspect <- terrain(bath_r, "aspect", unit="radians")
bath_r_hill <- shade(bath_r_slope, bath_r_aspect)

#### Slide 7 ####

### Load vector data (shapefiles) using sf
# Then convert to terra spatial vector and crop to region
# This uses the built in R pipe |> (version 4.1+)
# dsn = source folder (i.e. where the shapefiles are located)
# layer = shapefile name (no need for the extension)

coast <- st_read(dsn=paste0(shpf, "ne_10m_coastline"), layer="ne_10m_coastline") |> vect() |> crop(e_r)
rivers <- st_read(dsn=paste0(shpf, "ne_10m_rivers_lake_centerlines"), layer="ne_10m_rivers_lake_centerlines") |> vect() |> crop(e_r)

#### Slide 8 ####

### Create some vector data (place names)
# Create as a data.frame, include longitude and latitude, and the name
# Convert to terra vector
# Must assign a CRS to the data (this should be the same as the basemap)
places <- data.frame(lon=c(-6.84, -7.6, -7.98), lat=c(33.96, 33.55, 31.61), name=c("Rabat", "Casablanca", "Marrakech")) |> vect(geom=c("lon", "lat"), crs="EPSG:4326")

#### Slide 9 ####

### Plot
# Check out ?terra::plot for all options

# Set up colours for the plot
# Since using bathymetric data, its useful to set up separate water and land colours
col_w <- colorRampPalette(c("navy", "blue", "skyblue", "lightblue2", "lightblue"))
# For this example, we will use the built in terrain.colors() function for land

# Set up "breakpoints" 
# Break points define the point where the colour should change, and also allow complete control over the colours of the map
# This is necessary since we are using bathymetric data, and we want to ensure water and land are clearly delineated
# For DEMs and other data, you don't have to use breakpoints
# - Inf and Inf ensure that the entire range of data will be assigned a colour
# Since we are using bathymetic data, we set two breakpoint ranges - one for water, and one for land
bp_w <- c(-Inf, seq(-5000, 0, 100))
bp_l <- c(seq(0, 4000, 50), Inf)

#### Slide 10 ####

# Plot base map
# Usually a good idea to plot axes separately
# The map is built up in layers
# Most of these steps are optional, and all can be customised.

# Hill shade layer first
plot(bath_r_hill, maxcell=500000, col=grey(0:100/100), legend=FALSE, axes=FALSE)
# Plot elevation layer on top (use alpha, and add=TRUE)
# For the example, we are supressing the legend, but typically, you would want to show some kind of legend
plot(bath_r, maxcell=500000, breaks=c(bp_w, bp_l[-1]), col=c(col_w(length(bp_w)-1), terrain.colors(length(bp_l))), alpha=0.6, add=TRUE, legend=FALSE, axes=FALSE)
# Add axes
axis(1)
axis(2)
axis(3)
axis(4)
# Add a grid
grid()

#### Slide 11 ####

# Add vector data
lines(coast)
lines(rivers, col="skyblue")
# For the placenames - first add the marker points, then add the text
points(places, cex=2)
text(places, label="name", pos=4)
# Add a bounding box (to indicate an area that will have a zoomed in map)
# Use the rect function, and use coordinates to define area
rect(xleft=-8.2, ybottom=30.9, xright=-7.4, ytop=31.4, border="red", lwd=2)
# Add a scalebar
sbar(200, xy="bottomright", type="bar", divs=4)
# Add a north arrow
arrows(x0=-13.5, x1=-13.5, y0=36.8, y1=37.8, lwd=3)
# Add a box
box()


#### Slide 13 ####

############
### Zoomed in map (code is not shown in presentation, but is available below)
############

# Since the bathymetric data is relatively low resolution, we'll use ASTER DEM (v3) for the zoomed in map

############
# Information and download links for the ASTER DEM: https://lpdaac.usgs.gov/products/astgtmv003/
# Download relevant files and place in your DEM directory
# You may need to download several files for full coverage of your region of interest
# You may need to register to access the files (but it is free)
############

### Load DEMs - ASTERv3
### Important - Check file location and change below to match your system
# For example, if you put the downloaded data in a sub-folder within the "dems" directory, you might need to do something like this: dem1 <- rast(paste0(demf, "sub-folder-name/ASTGTMV003_N30W009_dem.tif"))
dem1 <- rast(paste0(demf, "ASTGTMV003_N30W009_dem.tif"))
dem2 <- rast(paste0(demf, "ASTGTMV003_N31W009_dem.tif"))
dem3 <- rast(paste0(demf, "ASTGTMV003_N30W008_dem.tif"))
dem4 <- rast(paste0(demf, "ASTGTMV003_N31W008_dem.tif"))
# Merge DEMs into one
dem <- mosaic(dem1, dem2, dem3, dem4, fun="mean")


# Set up zoomed in extent and crop DEM
e_z <- ext(c(-8.2, -7.4, 30.9, 31.4))
# Crop
dem_z <- crop(dem, e_z)

# Create hillshade
dem_z_slope <- terrain(dem_z, "slope", unit="radians")
dem_z_aspect <- terrain(dem_z, "aspect", unit="radians")
dem_z_hill <- shade(dem_z_slope, dem_z_aspect)

# Add vector data (if required)
# For this example, random sample data
samples <- data.frame(lon=runif(20, -8.2, -7.4), lat=runif(20, 30.9, 31.4), name=1:20) |> vect(geom=c("lon", "lat"), crs="EPSG:4326")


# Plot
plot(dem_z_hill, maxcell=5000000, col=grey(0:100/100), legend=FALSE, axes=FALSE)
plot(dem_z, maxcell=5000000, col=terrain.colors(length(bp_l)), alpha=0.6, legend=FALSE, add=TRUE, axes=FALSE)
# Add axes
axis(1)
axis(2)
axis(3)
axis(4)
# Add a grid
grid()
# Vector data (samples)
points(samples, col=rainbow(20), cex=2)
text(samples, label="name", pos=4)
# Add a scalebar
sbar(10, xy="bottomleft", type="bar", divs=4)
# Add a north arrow
arrows(x0=-8.17, x1=-8.17, y0=31.32, y1=31.38, lwd=3)
# Add a box
box(col="red", lwd=2)


#### Slide 14 ####

############
### Climate map
############

############
# Information and download links for Worldclim data: https://www.worldclim.org/data/index.html
# Download relevant files and place in the relevant directory
############


# Load data into R
### Important - change the file location

# Temperature
wc_t <- rast(***WorldClim*** "wc2.1_10m_tavg_08.tif") |> crop(e_r)
# Precipitation
wc_p <- rast(***WorldClim*** "wc2.1_30s_prec_01.tif") |> crop(e_r)


# Plot base map
# Usually a good idea to plot axes separately
# The map is built up in layers
# Hillshade layer first
plot(bath_r_hill, maxcell=500000, col=grey(0:100/100), legend=FALSE, axes=FALSE)
# Plot elevation layer on top (use alpha, and add=TRUE)
# For the example, we are supressing the legend, but typically, you would want to show some kind of legend
plot(wc_t, maxcell=500000, col=rev(heat.colors(100)), alpha=0.6, add=TRUE, legend=FALSE, axes=FALSE)
# Add axes
axis(1)
axis(2)
axis(3)
axis(4)
# Add a grid
grid()
# Add vector data
lines(coast)
lines(rivers, col="skyblue")
# For the placenames - first add the marker points, then add the text
points(places, cex=2)
text(places, label="name", pos=4)
# Add a scalebar
sbar(200, xy="bottomright", type="bar", divs=4)
# Add a north arrow
arrows(x0=-13.5, x1=-13.5, y0=36.8, y1=37.8, lwd=3)
# Add a box
box()


### Colour ramp for precipitation
raincol_all <- colorRampPalette(c("#FF3F00", "#FF8500", "#FFC700", "#FAFB00", "#B8FB00", "#6DF900", "#00F801", "#01F943", "#00F991", "#00FCD5", "#00EBFF", "#00AFFF", "#006CFF", "#0932FF", "#5732FE", "#A337FF", "#E33CFE"))

# Plot base map
# Usually a good idea to plot axes separately
# The map is built up in layers
# Hillshade layer first
plot(bath_r_hill, maxcell=500000, col=grey(0:100/100), legend=FALSE, axes=FALSE)
# Plot elevation layer on top (use alpha, and add=TRUE)
# For the example, we are supressing the legend, but typically, you would want to show some kind of legend
plot(wc_p, maxcell=500000, col=raincol_all(100), alpha=0.6, add=TRUE, legend=FALSE, axes=FALSE)
# Add axes
axis(1)
axis(2)
axis(3)
axis(4)
# Add a grid
grid()
# Add vector data
lines(coast)
lines(rivers, col="skyblue")
# For the placenames - first add the marker points, then add the text
points(places, cex=2)
text(places, label="name", pos=4)
# Add a scalebar
sbar(200, xy="bottomright", type="bar", divs=4)
# Add a north arrow
arrows(x0=-13.5, x1=-13.5, y0=36.8, y1=37.8, lwd=3)
# Add a box
box()

Part 2: Creating your own raster data and spatial analysis - Full R code

####################################
### BES R Masterclass 2022
# Part 2: Creating your own raster data and spatial analysis - Full R code (with notes!)
# This code follows the presentation, available on my blog.
# You may find additional help/info for creating maps and working with spatial data on my blog.
# Please feel free to get in touch with any queries (contact details on my blog)
# Benjamin Bell: https://www.benjaminbell.co.uk
####################################

####################################
# Recommended 16 GB ram minimum to run this code!
# Tested with R v4.2.1 on Windows 11 
####################################

#### You may need to refer to part 1 of this code ####

### Install required packages
install.packages(c("sf", "terra"))

### Load packages (You must do this each time you start R!)
library(sf)
library(terra)

### Setup data file paths

### Important - This is just an example. You must change folder locations for your system, depending on where you put the downloaded files, otherwise the code below will not work.
### Windows users should use "/" instead of "\" for file/folder locations

# Location of DEM/bathymetric data (change for your system)
demf <- "~/data/dems/"
# Location of vector data (change for your system)
shpf <- "~/data/shapefiles/"


#### Slide 15 ####

############
### Zoomed in map 
############

# Since the bathymetric data is relatively low resolution, we'll use ASTER DEM (v3) for the zoomed in map

############
# Information and download links for the ASTER DEM: https://lpdaac.usgs.gov/products/astgtmv003/
# Download relevant files and place in your DEM directory
# You may need to download several files for full coverage of your region of interest
# You may need to register to access the files (but it is free)
############

### Important - This DEM will form the basis for the temperature map we are going to create

# If you have already done these steps (from part 1), you can skip this part of the code - skip to #### Slide 16 ####

### Load DEMs - ASTERv3
### Important - Check file location and change below to match your system
# For example, if you put the downloaded data in a sub-folder within the "dems" directory, you might need to do something like this: dem1 <- rast(paste0(demf, "sub-folder-name/ASTGTMV003_N30W009_dem.tif"))
dem1 <- rast(paste0(demf, "ASTGTMV003_N30W009_dem.tif"))
dem2 <- rast(paste0(demf, "ASTGTMV003_N31W009_dem.tif"))
dem3 <- rast(paste0(demf, "ASTGTMV003_N30W008_dem.tif"))
dem4 <- rast(paste0(demf, "ASTGTMV003_N31W008_dem.tif"))
# Merge DEMs into one
dem <- mosaic(dem1, dem2, dem3, dem4, fun="mean")

# Set up zoomed in extent and crop DEM
e_z <- ext(c(-8.2, -7.4, 30.9, 31.4))
# Crop
dem_z <- crop(dem, e_z)

# Create hillshade
dem_z_slope <- terrain(dem_z, "slope", unit="radians")
dem_z_aspect <- terrain(dem_z, "aspect", unit="radians")
dem_z_hill <- shade(dem_z_slope, dem_z_aspect)


#### Slide 16 ####

# Generate new temperature data based on lapse rate equation
# This equation is based on my paper: https://www.tandfonline.com/doi/full/10.1080/15230430.2022.2046897 
# Available open access
# Relevant for the High Atlas, Morocco

# New mean annual temperature raster
ha_mean_temp <- 20.83754 + -0.004391263 * dem_z

# Look at DEM
dem_z

# Take a look at the new temperature raster layer
ha_mean_temp


#### Slide 17 ####

# Plot temperature raster
plot(dem_z_hill, maxcell=1000000, col=grey(0:100/100), legend=FALSE, axes=FALSE)
plot(ha_mean_temp, maxcell=1000000, col=rev(heat.colors(30)), range=c(0, 30), alpha=0.6, legend=FALSE, add=TRUE, axes=FALSE)
# Add axes
axis(1)
axis(2)
axis(3)
axis(4)
# Add a grid
grid()
# Add a box
box()

#### Slide 18 ####

# Generate new temperature data based on lapse rate equations
# These equations are based on my paper: https://www.tandfonline.com/doi/full/10.1080/15230430.2022.2046897 
# Available open access
# Relevant for the High Atlas, Morocco

# Max temperature
ha_max_temp <- 29.702586 + -0.006256 * dem_z
# Min temperature
ha_min_temp <- 12.609510 + -0.002544 * dem_z


# Plot min temperature
plot(dem_z_hill, maxcell=1000000, col=grey(0:100/100), legend=FALSE, axes=FALSE)
plot(ha_min_temp, maxcell=1000000, col=rev(heat.colors(30)), range=c(0, 30), alpha=0.6, legend=FALSE, add=TRUE, axes=FALSE)
# Add axes
axis(1)
axis(2)
axis(3)
axis(4)
# Add a grid
grid()
# Add a box
box()


# Plot
plot(dem_z_hill, maxcell=1000000, col=grey(0:100/100), legend=FALSE, axes=FALSE)
plot(ha_max_temp, maxcell=1000000, col=rev(heat.colors(30)), range=c(0, 30), alpha=0.6, legend=FALSE, add=TRUE, axes=FALSE)
# Add axes
axis(1)
axis(2)
axis(3)
axis(4)
# Add a grid
grid()
# Add a box
box()

# Note that all temperature plots use the same colourscale: col=rev(heat.colors(30)), range=c(0, 30)


#### Slide 19 ####

############
### Simple spatial analysis
############

# Find area with mean temp of between 6 and 8 degrees
# First create a matrix
# Convert values you don't want to NA to remove them (e.g. those below 6, and those above 8)
# Convert other values (that you want to keep) to the same value (e.g. 7)
mat <- matrix(c(0, 6, NA, 6, 8, 7, 8, 35, NA), nrow=3, byrow=TRUE)

# It might be helpful to view the matrix to understand which values are changing:
mat

# Reclassify raster and convert to a single vector
t68 <- classify(ha_mean_temp, mat, include.lowest=TRUE) |> as.polygons()


#### Slide 20 ####

# Plot the basemap (hillshade) and the vector
plot(dem_z_hill, maxcell=1000000, col=grey(0:100/100), legend=FALSE, axes=FALSE)
plot(t68, col="red", add=TRUE)
# Add axes
axis(1)
axis(2)
axis(3)
axis(4)
# Add a box
box()

#### Slide 21 ####

# Use to mask DEM
dem_mask <- mask(dem_z, t68)

# Take a look at the masked DEM
dem_mask

#### Slide 22 ####

# Now, also mask only north-east facing areas 0 to 90 degrees

# First - mask the aspect DEM (created when making the hillshade) using the temperature data
aspect_mask <- mask(dem_z_aspect, t68)

# Create matrix to reclassify data (remove areas we don't want)
### Important: aspect layer uses radians as units
# 0 to 1.5708 radians = 0 to 90 degrees 
mat2 <- matrix(c(0, 1.5708, 1, 1.5708, 10, NA), nrow=2, byrow=TRUE)

# It might be helpful to view the matrix to understand which values are changing:
mat2

# Reclassify raster and convert to a single vector
aspect_ne <- classify(aspect_mask, mat2, include.lowest=TRUE) |> as.polygons()

#### Slide 23 ####

# Plot
plot(dem_z_hill, maxcell=1000000, col=grey(0:100/100), legend=FALSE, axes=FALSE)
plot(aspect_ne, col="red", add=TRUE)
# Add axes
axis(1)
axis(2)
axis(3)
axis(4)
# Add a box
box()

Thanks for reading, please leave any comments below.