ggplotting shapefiles

By Carly Prior

Load these libraries

 1library(ggplot2)
 2library(maptools)
 3library(maps)
 4library(rgdal)
 5library(sp)
 6library(raster)
 7library(plyr)
 8library(dplyr)
 9library(knitr)
10library(tidyverse)

What is a shapefile?

A shapefile is a specific format for geospatial vector data describing points, lines or polygons. Shapefiles are most often used in GIS software like ArcGIS, but today we'll be using them in R as polygon layers on a ggplot map.

Load the shapefile with raster::shapefile

This shapefile is for the Laurentide Ice Sheet, which covered large portions of North America 18 thousand years ago during the Last Glacial Maximum.

 1# Note that while we are only loading the .shp file the .dbf, .shx, .sbn, .sbx, and .prj files all must be in the same folder as the .shp file to use the shapefile (more info in link below)
 2shapeLGM <- shapefile("LGM/ice18k_ICE.shp")
 3
 4#Exploring the class of our shapefile (link below for further info)
 5class(shapeLGM)
 6
 7## [1] "SpatialPolygonsDataFrame"
 8## attr(,"package")
 9## [1] "sp"
10
11#Exploring the shapefile's coordinate reference system (crs)
12crs(shapeLGM)
13
14## CRS arguments:
15##  +proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +x_0=0 +y_0=0
16## +ellps=clrk66 +units=m +no_defs

Note that this initial shape file has the Lambert conformal conic (lcc) projection

Load some basic maps with ggplot2::map_data

1world <- map_data("world")
2kable(head(world))
long lat group order region subregion
-69.89912 12.45200 1 1 Aruba NA
-69.89571 12.42300 1 2 Aruba NA
-69.94219 12.43853 1 3 Aruba NA
-70.00415 12.50049 1 4 Aruba NA
-70.06612 12.54697 1 5 Aruba NA
-70.05088 12.59707 1 6 Aruba NA
1states <- map_data("state") # We'll use this one later
2
3# Let's make a quick world map:
4ggplot() +
5  geom_polygon(data = world, aes(x = long, y = lat, group = group), fill = NA, color = "black", size = 0.25) +
6  coord_fixed(1.3) +
7  NULL

Now let's prepare our shapefile and plot it on the map

 1#Just for the example's sake we'll change the name of the file
 2testLGM <- shapeLGM
 3
 4#explicitly set the id so that later we can collapse the SpatialPolygonsDataFrame into a regular DataFrame
 5testLGM@data$id <- rownames(testLGM@data)
 6
 7#collapse to a regular DataFrame
 8testLGM.points <- fortify(testLGM, region = "id")
 9
10#Matches ids and joins the two dataframes in those spaces (explicit join)
11testLGM.df <- join(testLGM.points, testLGM@data, by="id")
12
13# Let's map it!
14ggplot() +
15   geom_polygon(data = testLGM.df, aes(x = long, y = lat, group = group), fill = "deepskyblue2", alpha = 0.2) +
16  coord_fixed(1.3) +
17  NULL

Well, that's not quite what we're looking for. This is because the shapefile has the lcc projection, and most of us are used to a rectangular projection, so let's change that and then use it as a layer over a map.

Note, you could also make a map with this specific projection!

 1#explicitly set the id so that later we can collapse the SpatialPolygonsDataFrame into a regular DataFrame
 2shapeLGM@data$id <- rownames(shapeLGM@data)
 3
 4#change the coordinate reference system to longlat projection using WGS84
 5shapeLGM.proj <- spTransform(shapeLGM, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
 6
 7#collapse to a regular DataFrame
 8lgmPoints <- fortify(shapeLGM.proj, region = "id")
 9
10#Matches ids and joins the two dataframes in those spaces (explicit join)
11lgm.df <- merge(lgmPoints, shapeLGM@data, by = "id")
12
13# Let's make a map!
14ggplot() +
15  geom_polygon(data = world, aes(x = long, y = lat, group = group), fill = NA, color = "black", size = 0.25) +
16   geom_polygon(data = lgm.df, aes(x = long, y = lat, group = group), fill = "deepskyblue2", alpha = 0.2) +
17  coord_fixed(1.3) +
18  NULL

Wow! What a map. But maybe we really care about the area specifically affected by this particular ice sheet.

1ggplot() +
2  geom_polygon(data = subset(world, long > -180 & long  30), aes(x = long, y = lat, group = group), fill = NA, color = "black", size = 0.25) +
3  geom_polygon(data = lgm.df, aes(x = long, y = lat, group = group), fill = "deepskyblue2", alpha = 0.2) +
4  coord_fixed(1.3) +
5  NULL

Super cool! Now, I specifically learned these methods to produce maps that show how the ice sheet encroached on an eastern North American species' extant geographic range, implying a historical range contraction as the ice sheets expanded and subsequently a historical range expansion as the ice sheets shrank.

1#Let's look at the ice sheet over a very specific area in the United States
2ggplot() +
3  geom_polygon(data = subset(states, long > -100 & long < -75), aes(x = long, y = lat, group = group), fill = NA, color = "black", size = 0.25) +
4   geom_polygon(data = lgm.df, aes(x = long, y = lat, group = group), fill = "deepskyblue2", alpha = 0.2) +
5  coord_fixed(1.3) +
6  NULL

Yikes.

1# Let's try again and reduce the size of the shape file like we did with the states object
2ggplot() +
3  geom_polygon(data = subset(states, long > -100 & long  -100 & long < -75), aes(x = long, y = lat, group = group), fill = "deepskyblue2", alpha = 0.2) +
4  coord_fixed(1.3) +
5  NULL

OK, still not great but let's subset on the latitude as well.

1ggplot() +
2  geom_polygon(data = subset(states, long > -100 & long  -100 & long < -75 & lat < 50), aes(x = long, y = lat, group = group), fill = "deepskyblue2", alpha = 0.2) +
3  coord_fixed(1.3) +
4  NULL

Much worse. Why?

Well our polygon has no defined edge at 50N, so there is no way to properly display what we asked for.

Let's fix that!

 1# We can reduce the dataframe to just longitude and latitude
 2lgmpoly <- lgm.df[,2:3]
 3
 4# Then subset that polygon to the desired area
 5lgmpoly.1  -100 & long < -75 & lat < 50)
 6
 7# Finally, serially add points along the latitude = 50 line
 8for (i in -100:-75){
 9  lgmpoly.1[nrow(lgmpoly.1) + 1,]  -100 & long < -75), aes(x = long, y = lat, group = group), fill = NA, color = "black", size = 0.25) +
10    geom_polygon(data = lgmpoly.1, aes(x = long, y = lat), fill = "deepskyblue2", alpha = 0.2) +
11  coord_fixed(1.3) +
12  NULL

A thing of beauty! On top of this map you can easily apply other points, lines and polygons!