Goal

In this document we will use tmap to visualize the landscape around a geographic point and create several circular buffers around the point. Mapping with tmap is a relatively easy way to make maps, and it’s becoming increasingly popular. Although this exercise is meant to be a quick starting point and relatively accessible for beginners, tmap seems efficient at working with large datasets, very customizable and aesthetic, and excellent at producing online interactive maps.

Loading required packages, importing data, and data cleaning

For mapping, we will be using sf to manage our spatial objects, terra for raster data manipulation, tmap for making the maps, and readr and tidyverse for data management.

library(readr)
library(tidyverse)
library(sf)
library(terra)
library(landscapemetrics)
# rinstall.packages("tmap", repos = c("https://r-tmap.r-universe.dev", "https://cloud.r-project.org"))
library(tmap)

Read in site coordinates from a data frame

For visualization purposes, we want the example site we map to show all 6 of our categorical landcover types within a 2000 m circular buffer. So, here we load in a data set containing the lat/lon coordinates of the site we want to map.

da <- read.csv("./data/site_coordinates.csv")


We now have the lat and lon data for our one point.
Now we need to make R recognize it as a geographic point with the sf package. This point is in the province of Alberta, Canada in lat/lon format, so the first thing we need to do is define the coordinate reference system of the map and then assign that reference system to the specified columns for the lon/lat coordinates (x,y).

#set the appropriate coordinate reference system (crs)
crs_wgs84 <- st_crs(4326) # WGS84 has EPSG code 4326

#make the sf
Site <- st_as_sf(da, coords = c("lon", "lat"), crs=crs_wgs84)

  Let’s double check it worked

sf::st_crs(Site)$input 
## [1] "EPSG:4326"


And we need to transfrom the lat/lon coordinates to UTM coordinates to match the landcover data we will work with later.

# transform to utm coords so it matches the projection of the rasters
crs<- st_crs(3857) #define the crs for UTM 
Site <- sf::st_transform(Site, crs = crs)

#check projection
sf::st_crs(Site)$input
## [1] "EPSG:3857"


Now if we want to plot our site, it looks like this . . . which is not very informative at all. We need to put this point on a map.

plot(Site, col="black", pch=19)


But first let’s create some circular buffers around our site ranging from 150 - 2000 m.

#create buffers around sites
buffer150 <- st_buffer(Site, dist = 150)
buffer510 <- st_buffer(Site, dist = 510)
buffer990 <- st_buffer(Site, dist = 990)
buffer2000 <- st_buffer(Site, dist = 2000)
buffer5000 <- st_buffer(Site, dist=5000)

Working with raster maps

The map were working with is a raster file in .tif format, which we can read in with the terra::rast function.

#import raster
map15 <- terra::rast("./data/AlbertaLandcoverMap2015_5km_around_site.tif")

#give the map an appropriate name
names(map15)<- "2015 Alberta Landcover Map"

#take a look at it's description
map15
## class       : SpatRaster 
## dimensions  : 334, 333, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : -12546030, -12536040, 7075508, 7085528  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
## source      : AlbertaLandcoverMap2015_5km_around_site.tif 
## name        : 2015 Alberta Landcover Map 
## min value   :                          1 
## max value   :                          9


This map is pretty small, so we can plot the map and our site directly in base R if we want. Although, this doesn’t tell us much at the moment.

plot(map15)
plot(Site, add=TRUE,col="black", pch=19)

First, let’s crop the map to the largest buffer we made earlier. If we were using a larger raster, we would likely want to crop the map first before plotting it because large rasters can sometimes be difficult or slow to visualize in R.

#create a bounding box around the largest buffer
mybbox <- terra::ext(buffer2000)*1.35
map15 <- crop(map15, mybbox, extend=TRUE)
map15
## class       : SpatRaster 
## dimensions  : 180, 180, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : -12543720, -12538320, 7077818, 7083218  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
## source(s)   : memory
## varname     : AlbertaLandcoverMap2015_5km_around_site 
## name        : 2015 Alberta Landcover Map 
## min value   :                          1 
## max value   :                          9


Looks good, and we see that the coordinate references matches the one we used to transform the site coordinates earlier.

Creating maps with tmap

Now we can start creating our map with tmap.

Let’s take a look.

We call the map with tmap_shape and tell it what type data the map has with tmap_raster.

We can then plot our site on top of that using another tmap_shape and specify what kind of display withtm_dots, and it looks like we’ve got some work to do. We have 9 landcover types, although we’re only concerned with 6, and landcover type number 2 does not appear in this map. Also we should identify the landcover types and choose some more appropriate colors.

m1 <- 
  tm_shape(map15) +
  tm_raster(style="cat")+ #style = "cat" means a categorical raster
  tm_shape(buffer2000) + tm_borders("black", lwd=2)+
  tm_shape(Site) + tm_dots(size=.25, col="black")+
  tm_shape(buffer150) + tm_borders("black", lwd=2)+
  tm_shape(buffer990) + tm_borders("black", lwd=2)+
  tm_shape(buffer510) + tm_borders("black", lwd=2)
m1

First, let’s identify these 1-9 landcover types.

Landcover Key
“1” = “water”
“2” = “barren”
“3” = “developed”
“4” = “grassland”
“5” = “wetland”
“6” = “non_flowering_crop”
“7” = “flowering_crop”
“8” = “canola”
“9” = “forest”

First let’s plot our site, landcover types, and buffers more aestheitcally, and then we can get rid of those landcovers we’re not interested in.

m1 <- tm_shape(map15)+
  tm_raster(style="cat",
            labels = c("Water", "Developed", "Grassland",  "Wetland" ,"Non-flowering Crop","Flowering Crop", "Canola", "Forest"), 
            palette = c("darkblue", "grey55", "darkolivegreen2", "lightskyblue3", "tan4", "lightpink2", "gold3", "olivedrab"))+ 
tm_layout(legend.outside = TRUE, legend.show = TRUE, frame = FALSE)+
tm_shape(Site) + tm_dots(size=.25, col="black")+
tm_shape(buffer150) + tm_borders("black", lwd=2)+
tm_shape(buffer990) + tm_borders("black", lwd=2)+
tm_shape(buffer510) + tm_borders("black", lwd=2)+
tm_shape(buffer2000) + tm_borders("black", lwd=2)+
tm_layout(legend.frame = FALSE)
  #tm_compass(position = c("left", "bottom"))+
  #tm_scale_bar(position = c("left", "bottom"))+
m1

Looks pretty good!

Now let’s remove the landcover types we’re not interested in: non-flowering crop and water.

To do that we need to reclassify (or “substitute” using subst) landcover types 1 and 6, and we’ll convert them to a new a new class of 10. We can “hide” this class in the legend if we choose by setting an empty label and choosing to display this new class in white, so it appears as an empty section of map.

Now lets compare our maps, and add a scale bar and north arrow

#Reclassify landcover types 1 and 6 to 10
map15_v2 <- subst(map15, 1, 10, raw =TRUE)
map15_v2 <- subst(map15_v2, 6, 10, raw =TRUE)

m2 <- tm_shape(map15_v2)+
  tm_raster(style="cat",labels = c("Developed", "Grassland",  "Wetland" ,"Flowering Crop", "Canola", "Forest", ""), palette = c( "grey55", "darkolivegreen2", "lightskyblue3",  "lightpink2", "gold3", "olivedrab","white"))+ 
  tm_layout(legend.outside = TRUE)+
tm_shape(Site) + tm_dots(size=.25, col="black")+
tm_shape(buffer990) + tm_borders("black", lwd=2)+
tm_shape(buffer510) + tm_borders("black", lwd=2)+
tm_shape(buffer150) + tm_borders("black", lwd=2)+
tm_shape(buffer2000) + tm_borders("black", lwd=2)+
  # tm_compass(position = c(0.125, 0.1))+
  tm_scale_bar(position = c(0.07, 0.08))+
  tm_compass(position =c("right", "bottom"))+
  #tm_scale_bar(position =c("left", "bottom"))+
  tm_layout(legend.show = TRUE, frame = FALSE, legend.frame = FALSE,
            outer.margins = c(0,0,0,0))

fig <- tmap_arrange(m1, m2, ncol = 1)
fig

Looks pretty good! And if we want, we can save this map as a pdf with the code below.

#tmap_save(fig, "Map_Landcover_Comparison_Buffers.pdf", height = 8, width = 7 )