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))
| -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!