The other day, I wanted to tweet out a map showing the distribution of some wood frog tissue samples compared to the entire range of the species. I’m not much for GIS and I didn’t need anything complicated so I wanted to plot in R. If you find yourself in a similar situation needing to plot a simple range map with GGplot, then this post is for you.

This post WILL NOT explain how to estimate a species range. There’s a lot more involved in that process and there are dedicated packages like rangemappr to help.

In this post, we will take a pre-made polygon of a species range and slap it onto a basemap.

The first task is to find a polygon of the range of your favorite species. Unless you already have a shapefile in your possession, one of the easiest places to acquire one is through the IUCN database. You can follow the links to go down the taxonomic rabbit hole for specific species. For this tutorial I’ll be using the red salamander (Pseudotriton ruber), just because I think that they are gorgeous.

Here’s the  IUCN page for the red salamander. In the top right of the page, you’ll see a drop down button labelled, “Download”. You will want to select and download the “Range data – Polygons (SHP)” option. You will need to agree to some conditions in order to download the file, but it is quick and painless. Once you download the zipped folder, extract it into your project directory.

Here are the packages we will be using.

library(tidyverse)
library(rgdal)
library(ggspatial)

First, we need to do a bit of work to get the shapefiles into R. First, we will use rgdal package to read in the shapefiles. Then, we will coerce that object into a dataframe that can play nicely with GGplot.

# Read in the shapefile
PSRU_shape <- readOGR(dsn = "./PSRU_range", layer = "data_0")
# Coerce into a dataframe to play nicely with GGplot
PSRU_shape_df <- fortify(PSRU_shape)

Now we can start mapping. We will use the ggmap package in GGplot that all comes bundled in the tidyverse. First, we define regions for our basemap. In this case, red salamanders range is entirely within the lower 48 United States. Then we plot the basemap and layer on the range map on top. It’s really that simple.

Low48_map <- map_data("state") # Create basemap form GGplot

ggplot(Low48_map, aes(x = long, y = lat, group = group)) +
geom_polygon() +
geom_polygon(data = PSRU_shape_df, fill = "orangered")

But we can make this look a lot nicer. First off, we can ditch the ugly GGplot default theme and add our own color scheme. We can also force GGplot to use a polyconic projection (you’ll need to make sure that you also have the mapproj package installed).

ggplot(Low48_map, aes(x = long, y = lat, group = group)) +
geom_polygon(fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = PSRU_shape_df, fill = "darkorange3", alpha = 0.8) +
theme_void() +
theme(panel.background = element_rect(fill = "cornsilk")) +
coord_map(projection = "polyconic")

Beyond the basics:

For species with larger ranges, you might want to incorporate more of the continent. We can easily do that by defining larger regions for our basemap. Here, I am still plotting the Lower 48 states basemap to get the state borders.

NAm_map <- map_data("world", region = c("Mexico", "Canada")) # Create basemap form GGplot

ggplot(NAm_map, aes(x = long, y = lat, group = group)) +
geom_polygon(fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = Low48_map, fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = PSRU_shape_df, fill = "darkorange3", alpha = 0.8) +
theme_void() +
theme(panel.background = element_rect(fill = "cornsilk")) +
coord_map(projection = "gilbert")

One issue is that it is extremely difficult to crop in if, for instance, we don’t need all of Canada’s norther islands. Using xlim and ylim truncates the polygons in weird ways. And you should NEVER crop with xlim/ylim, anyway.

ggplot(NAm_map, aes(x = long, y = lat, group = group)) +
geom_polygon(fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = Low48_map, fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = PSRU_shape_df, fill = "darkorange3", alpha = 0.8) +
theme_void() +
theme(panel.background = element_rect(fill = "cornsilk")) +
coord_map(projection = "gilbert") +
xlim(-125, -60) +
ylim(25, 64)

The better way to crop is to set the xlim and ylim of the clipping mask with coord_cartesian(). But, then we lose the projection.

ggplot(NAm_map, aes(x = long, y = lat, group = group)) +
geom_polygon(fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = Low48_map, fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = PSRU_shape_df, fill = "darkorange3", alpha = 0.8) +
theme_void() +
theme(panel.background = element_rect(fill = "cornsilk")) +
coord_map(projection = "mercator") +
coord_cartesian(xlim = c(-125, -60), ylim = c(25, 52))

There are a handful of work-arounds to the problem.

First, there is a code workaround here, but it is a bit complicated.

Second, you can simply export the image at an appropriate size to get a resulting ratio that works. One potential issue with this is that GGplot exports ALL of the information in the PDF, even the polygons outside of the image. In this case, that means a lot of data points to outlines all those tiny Arctic Islands. This can make for some really large PDF files. A very slick solution to this is to rasterize the maps with ggraster before exporting the images. You can still export the image at twice or three times the size to retain high resolution, but the file size will be much smaller.

Third, you could export the uncropped version as a pdf (which is a vector graphic) and open it in a vector graphics program like Adobe Illustrator or Inkscape. From there, you can crop in where ever you want.