How do you plot a grid over shapefile polygons using R? A straightforward procedure with a surprising lack of clear information on the internet. This post will cover how to create a grid using the Terra package, how to plot shapefiles with ggplot2 and the tidyterra package, and give some insights into spatial visualization for your projects! I’ve also included all the code as an annotated r-script you can download below. Let’s hop in!
Background
I spent the last week developing figures for my Master’s thesis on the World Turtles. However, I ran into numerous issues due to out-of-date tutorials, changes in API requirements, and the shift to the requirement of SF objects for spatial data. But I am happy to say that my figures are made, embedded, and ready for publishing!
When I posted about these issues on Social Media, a biologist friend of mine asked for help with their figures. Specifically, how do you get a grid over a SF object in R using the ggplot2 package? Luckily, the knowledge for how was fresh in my mind and I happened to need content for the weekly Learn Adventurously Newsletter (which you should sign up for using the form below).
Sign up for the Learn Adventurously Newsletter
What are we creating?
The desired figure takes a simple shapefile and plots a grid on top. Here is the final product that we are working towards.
A few things to note:
- I am using a shapefile I already have loaded. The sharp-eyed may notice that it is the same shapefile as the figure I showed at the top of the blog, however, we have cropped it (arbitrarily) to include everything east of 20 degrees.
- This tutorial should work with any Simple Feature object.
- You can achieve this visually using the panel.grid.major theme arguments as part of the ggplot2 package. However, this tutorial creates an actual rasterized grid based off the shapefile. The importance of this is explained later.
- We are using ggplot2 to create these visuals. This post does not cover the basics, but you can check out the free webinar on the fundamentals of ggplot2 or my self-paced online course covering only ggplot2. Luckily, this code is very straightforward if you know the basics!
Eventually, i’ll include this material for a “Spatial Data in R Course”. Stay tuned for that! If you want to download the full annotated code for this figure, all I need is your email:
Let’s break this down!
First, we load all the packages we need.
library(tidyverse)
library(sf)
library(terra)
library(tidyterra)tidyverse is a staple that should be loaded for almost all projects. It includes ggplot2.
sf stands for simple features and is the standard for working with spatial polygons in R
terra is used for all sorts of spatial analyses. It is also the replacement for the popular raster package.
Finally, tidyterra allows us to work with spatial data from Terra using a tidy framework. For example, plotting with ggplot2.
Let’s first plot our shapefile object sf so we can see what we are working with. Again, I already had this loaded in my r environment, but this tutorial should work with any sf shapefile.
ggplot(sf)+
geom_sf()
Our shapefile includes the entire world. You can crop the shapefile using a function such as st_crop, but we won’t need to here.
The goal was to create a grid that goes over this shapefile. There are ways to do this visually via ggplot2 themes, but creating an actual grid is a far more useful skill to learn. Often, I create grids based on a shapefile to create a presence-absence matrix of taxa occurrence, use it for randomization to create null models, or extract data such as climate, elevation, or taxa richness.
Luckily, its pretty easy to do!
I use the simple st_make_grid() function to make a grid from a shapefile. Let’s create that grid and then plot it using an additional geom_sf() layer. We’ll set the fill to NA to make them transparent and the line color to red.
grid <- st_make_grid(sf, cellsize = 10)
ggplot(sf)+
geom_sf()+
geom_sf(data = grid, fill = NA, col = "red")
st_make_grid() took 2 arguments, first, the shapefile and the other was cell size. The shapefile argument sets the extent (or the boundaries) of where a grid is created. Since we fed it an sf object, it uses the entire extent of that object! In this case, the entire world.
The cell size argument is the resolution, and it uses the units found in the SF object. Since here we are using decimal lat/lon, a resolution of 10 will create cells of 10 degrees by 10 degrees. We can see this clearly as between 0 degrees and 50 degrees E, there are 5 cells per row 🙂
1 degree is ~ 111km, so a 10×10 grid has each cell measure ~1,111km by ~1,111 km.
Note the exact distance calculated can change depending on projection and north/south distance from the equator. The joys of cartography and GIS!
We now have all the building blocks needed for our visualization. We just need to fine-tune our arguments to make this map good. First, we’re going to set an offset in the grid function. Offset determines where the bottom left corner of our grid will be. I often use this to get around cropping the underlying shapefile when creating visualizations en masse.
Let’s set our grid to start at (-20, -40) to capture all the continents of the Eastern Hemisphere. I also chose these coordinates, because they would line up the grid lines on major axes.
grid <- st_make_grid(sf, #plot
cellsize=10), #Resolution
offset = c(-20,-40)) #offset 
Our next step would be to crop the visualization to our grid. To do so we use the scale_x_continuous and scale_y_continuous functions, setting the limits arguments to our grid. We can manually set the minimum and maximum limits to whatever we want. For example let’s set the x limit between -20 and 50, while we set our y scale limits to between -40 and 40. This will crop our image to just Africa.
ggplot(sf)+ #base plot
geom_sf()+ #base plot
geom_sf(data = grid, fill = NA, col = "red")+ #grid plot
scale_x_continuous(limits = c(-20,50))+ #X limits
scale_y_continuous(limits = c(-40,40)) #y limits
But, for my work, I often need to produce dozens of visualizations automatically. Manually setting those boundaries sucks! Especially if you need to make a minor change and now need to update 30 different ggpltos. Instead, let’s set these limits using a bounding box.
In this use case, we’ll create a bounding box from the grid so that the visual is cropped to the scale of the grid. To do so we use the st_bbox() function from SF which takes the extent of an SF object and returns the minimum and maximum values for each of the x and y scales.
bb <- st_bbox(grid) #bounding box
#Code Output
> bb <- st_bbox(grid) #bounding box
> bb
xmin ymin xmax ymax
-20 -40 160 60 Well hey! Would you know it? The bounding box gives us the min and max limits for x and y. Those are the exact values we need for scale_x_continuous and scale_y_continuous! Let’s add that to our visualization using simple subsetting.
ggplot(sf)+ #base plot
geom_sf()+ #base plot
geom_sf(data = grid, fill = NA, col = "red")+ #grid plot
scale_x_continuous(limits = c(bb[1],bb[3]))+ #X limits
scale_y_continuous(limits = c(bb[2],bb[4]))+ #y limits
Common mistake alert! Pay attention to which values from the bounding box are min and max. My students regularly assume the output from st_bbox is organized as xmin, xmax, ymin, ymax; in reality it is xmin, ymin, xmax, ymax.
You may think this is done, but I want to extend an extra layer of customization by adjusting the theme elements of a ggplot. What I want to do is remove the background and major panel lines. The major panel lines in this visualization are not problematic, because they just so happen to align with the grid. However, let’s adjust our offset to visualize the issue you very well may have.
grid <- st_make_grid(sf, #plot
cellsize=10, #Resolution
offset = c(-23,-45))
bb <- st_bbox(grid) #bounding box
ggplot(sf)+ #base plot
geom_sf()+ #base plot
geom_sf(data = grid, fill = NA, col = "red")+ #grid plot
scale_x_continuous(limits = c(bb[1],bb[3]))+ #X limits
scale_y_continuous(limits = c(bb[2],bb[4]))
This visualization is very confusing to look at thanks to the major panel lines. You can adjust them to line up perfectly with our gridlines, but it is much easier to just hide them. I tend to follow the visualization practice of less is more. If it doesn’t add anything to your visualization, remove it!
To do so, we’ll set the color of panel.grid equal to NA. While we’re at it, let’s set the panel background to be clear as well. This plot below is still plotting using the adjusted offset from earlier. Themes can be confusing for newbies to ggplot2, but I do have it documented in my Free Fundamentals of R for Biologists e-textbook.
ggplot(sf)+ #base plot
geom_sf()+ #base plot
geom_sf(data = grid, fill = NA, col = "red")+ #grid plot
scale_x_continuous(limits = c(bb[1],bb[3]))+ #X limits
scale_y_continuous(limits = c(bb[2],bb[4]))+ #y limits
theme(panel.grid = element_line(color = NA),
panel.background = element_rect(fill = NA)) 
And just like that, we have an SF object plotted with a grid overlaid! From here you can adjust the resolution, the x/y limits, and all the colors/fill you can imagine. Of course, there are numerous ways to perform this operation, but this is the way I like to do it!
You can feel free to copy and paste the code from this blog, however, If you want to download the fully annotated r script, all you need is to fill out this form below. Donations are entirely optional but greatly appreciated!