Maps in R
1 Introduction
R has the ability through the maps package and the base graphics to generate maps, but such “out-of-the-box” maps, like other base graphics-generated illustrations, these may not be suitable for immediate publication. Other options exit, for example, using packages like ggmaps, mapview, or mapmate. However, the ability of the sp and rgdal packages to handle the diverse kinds of spatial data, including shapefiles, provides a facility for creating publishable maps. The examples here use a set of shapefiles downloaded from the Natural Earth web page [http://www.naturalearthdata.com].
This first example recreates Figure 1 from Marlon, J.R., et al., 2016, Reconstructions of biomass burning from sediment-charcoal records to improve data–model comparisons. Biogeosciences 13:3225-3244. [http://www.biogeosciences.net/13/3225/2016/]. There are two general steps in the production of this map, the first involves 1) projecting the Natural Earth shapefiles into the Robinson projection (which are reusable for other applications), and the second the projection of the data and assembly into the finished map.
2 Set up outlines and polygons for global (Robinson projection) maps
2.1 Read Natural Earth shapefiles
The first step is to read a set of shapefiles downloaded from the Natural Earth web page [http://www.naturalearthdata.com], and project these into the Robinson projection. Two CRAN packages, rnaturalearth and rnaturalearthdata can be used to manage and download the data within R, but the Natural Earth web pages are worth looking at.
Begin by loading the appropriate packages:
library(sp)
library(rgdal)
library(rgeos)
library(raster)
library(rasterVis)
library(classInt)
library(RColorBrewer)Now read the shape files, including those for global coastlines, adminstrative units (borders), large lakes, and a graticule. Set the filenames:
# Natural Earth shape files -- global (Robinson) projections
# get shapefiles from http://www.naturalearthdata.com
shape_path <- "/Users/bartlein/Projects/ESSD/data/RMaps/source/"
coast_shapefile <- paste(shape_path, "ne_50m_coastline/ne_50m_coastline.shp", sep="")
ocean_shapefile <- paste(shape_path, "ne_50m_ocean/ne_50m_ocean.shp", sep="")
admin0_shapefile <- paste(shape_path, "ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp", sep="")
admin1_shapefile <- paste(shape_path, "ne_50m_admin_1_states_provinces_lakes/ne_50m_admin_1_states_provinces_lakes.shp", sep="")
lakes_shapefile <- paste(shape_path, "ne_50m_lakes/ne_50m_lakes.shp", sep="")
bb_shapefile <- paste(shape_path, "ne_50m_graticules_all/ne_50m_wgs84_bounding_box.shp", sep="")
grat30_shapefile <- paste(shape_path, "ne_50m_graticules_all/ne_50m_graticules_30.shp", sep="")Before reading each file, determine what kind it is (e.g. lines vs. polygons), then read and plot it. Using the ogrInfo() and readOGR functions from rgdal eliminates the necessity of selecting the appropriate function in maptools. In the ‘plot()’ functions below, use the col argument for lines, bor for polygon border colors.
# find out kind of shapefile (lines vs. polygons)
layer <- ogrListLayers(coast_shapefile)
ogrInfo(coast_shapefile, layer=layer)## Source: "/Users/bartlein/Projects/ESSD/data/RMaps/source/ne_50m_coastline/ne_50m_coastline.shp", layer: "ne_50m_coastline"
## Driver: ESRI Shapefile; number of rows: 1428
## Feature type: wkbLineString with 2 dimensions
## Extent: (-180 -85.19219) - (180 83.59961)
## CRS: +proj=longlat +datum=WGS84 +no_defs
## LDID: 87
## Number of fields: 2
## name type length typeName
## 1 scalerank 12 10 Integer64
## 2 featurecla 4 32 String## OGR data source with driver: ESRI Shapefile
## Source: "/Users/bartlein/Projects/ESSD/data/RMaps/source/ne_50m_coastline/ne_50m_coastline.shp", layer: "ne_50m_coastline"
## with 1428 features
## It has 2 fields
## Integer64 fields read as strings: scalerank## Object of class SpatialLinesDataFrame
## Coordinates:
## min max
## x -180.00000 180.00000
## y -85.19219 83.59961
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## scalerank featurecla
## 0:1428 Coastline:1428
Get the Proj4 string for later use>
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"Read the ocean polygons.
# find out kind of shapefile (lines vs. polygons)
layer <- ogrListLayers(ocean_shapefile)
ogrInfo(ocean_shapefile, layer=layer)## Source: "/Users/bartlein/Projects/ESSD/data/RMaps/source/ne_50m_ocean/ne_50m_ocean.shp", layer: "ne_50m_ocean"
## Driver: ESRI Shapefile; number of rows: 1
## Feature type: wkbPolygon with 2 dimensions
## Extent: (-180 -85.19219) - (180 90)
## CRS: +proj=longlat +datum=WGS84 +no_defs
## LDID: 87
## Number of fields: 2
## name type length typeName
## 1 scalerank 12 10 Integer64
## 2 featurecla 4 32 String## OGR data source with driver: ESRI Shapefile
## Source: "/Users/bartlein/Projects/ESSD/data/RMaps/source/ne_50m_ocean/ne_50m_ocean.shp", layer: "ne_50m_ocean"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings: scalerank## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -180.00000 180
## y -85.19219 90
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## scalerank featurecla
## 0:1 Ocean:1
Read the other shape files the same way. (Note that the summary output is suppressed here.)
layer <- ogrListLayers(admin0_shapefile)
ogrInfo(admin0_shapefile, layer=layer)
admin0_poly <- readOGR(admin0_shapefile, layer=layer)
summary(admin0_poly)
layer <- ogrListLayers(admin1_shapefile)
ogrInfo(admin1_shapefile, layer=layer)
admin1_poly <- readOGR(admin1_shapefile, layer=layer)
summary(admin1_poly)
layer <- ogrListLayers(lakes_shapefile)
ogrInfo(lakes_shapefile, layer=layer)
lakes_poly <- readOGR(lakes_shapefile, layer=layer)
summary(lakes_poly)
lrglakes_poly <- lakes_poly[as.numeric(lakes_poly$scalerank) <= 2 ,]
layer <- ogrListLayers(grat30_shapefile)
ogrInfo(grat30_shapefile, layer=layer)
grat30_lines <- readOGR(grat30_shapefile, layer=layer)
summary(grat30_lines)
layer <- ogrListLayers(bb_shapefile)
ogrInfo(bb_shapefile, layer=layer)
bb_poly <- readOGR(bb_shapefile, layer=layer)
summary(bb_poly)
bb_lines <- as(bb_poly, "SpatialLines")Plot everything (except ocean).
# plot everything
plot(coast_lines, col="black")
plot(admin0_poly, bor="gray", add=TRUE)
plot(admin1_poly, bor="pink", add=TRUE)
plot(lakes_poly, bor="lightblue", add=TRUE)
plot(lrglakes_poly, bor="blue", add=TRUE)
plot(grat30_lines, col="lightblue", add=TRUE)
plot(bb_lines, col="black", add=TRUE)
plot(coast_lines, col="purple", add=TRUE)
2.2 Project the shapefiles
Next, transform (project) the individual shapefiles into the Robinson projection using spTransform. Set the proj4 string for the projectdion using the CRS() function.
# set Robinson CRS
unproj_crs <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
unproj_crs## CRS arguments:
## +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0## CRS arguments: +proj=robin +lon_0=0w +ellps=WGS84Now do the projections.
# do projections
bb_poly_proj <- spTransform(bb_poly, robin_crs)
coast_lines_proj <- spTransform(coast_lines, robin_crs)
admin0_poly_proj <- spTransform(admin0_poly, robin_crs)
admin1_poly_proj <- spTransform(admin1_poly, robin_crs)
lakes_poly_proj <- spTransform(lakes_poly, robin_crs)
grat30_lines_proj <- spTransform(grat30_lines, robin_crs)
lrglakes_poly_proj <- spTransform(lrglakes_poly, robin_crs)Plot the projected shapefiles.
# plot projected files
plot(bb_poly_proj, col="gray95")
plot(admin0_poly_proj, col="white", bor="pink", add=TRUE)
plot(admin1_poly_proj, col="white", bor="pink", add=TRUE)
plot(lakes_poly_proj, bor="lightblue", add=TRUE)
plot(lrglakes_poly_proj, bor="purple", add=TRUE)
plot(coast_lines_proj, col="black", add=TRUE)
plot(grat30_lines_proj, col="lightblue", add=TRUE)
plot(bb_poly_proj, bor="black", add=TRUE)
Because the map we’re creating will use the borders and lakes only as outlines (as opposed to polygons that might be filled with a specific color), those polygons can be converted to SpatialLines objects.
# convert polygons to spatial lines
admin0_lines_proj <- as(admin0_poly_proj, "SpatialLines")
admin1_lines_proj <- as(admin1_poly_proj, "SpatialLines")
lakes_lines_proj <- as(lakes_poly_proj, "SpatialLines")
lrglakes_lines_proj <- as(lrglakes_poly_proj, "SpatialLines")
bb_lines_proj <- as(bb_poly_proj, "SpatialLines")Plot the SpatialLines objects
# test the SpatialLines shapefiles
plot(bb_poly_proj, col="gray95")
plot(coast_lines_proj, col="green", add=TRUE)
plot(admin0_lines_proj, col="lightblue", add=TRUE)
plot(admin1_lines_proj, col="lightblue", add=TRUE)
plot(lakes_lines_proj, col="blue", add=TRUE)
plot(lrglakes_lines_proj, col="purple", add=TRUE)
plot(grat30_lines_proj, col="gray", add=TRUE)
plot(coast_lines_proj, col="black", add=TRUE)
plot(bb_lines_proj, col="black", add=TRUE)
2.3 Write out the projected shapefiles
Next, write out the projected shapefiles, first setting the output path.
# write out the various shapefiles
outpath <- "/Users/bartlein/Projects/ESSD/data/RMaps/derived/glrob_50m/"
outshape <- coast_lines_proj
outfile <- "glRob_50m_coast_lines"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialLinesDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)It’s always good practice to test whether the shapefile has ideed been written out correctly. Read it back in and plot it.
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/bartlein/Projects/ESSD/data/RMaps/derived/glRob_50m/glRob_50m_coast_lines", layer: "glRob_50m_coast_lines"
## with 1428 features
## It has 1 fields
Write out the other shapefiles.
outshape <- bb_poly_proj
outfile <- "glRob_50m_bb_poly"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)
outshape <- bb_lines_proj
outfile <- "glRob_50m_bb_lines"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialLinesDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)
outshape <- admin0_poly_proj
outfile <- "glRob_50m_admin0_poly"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)
outshape <- admin0_lines_proj
outfile <- "glRob_50m_admin0_lines"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialLinesDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)
outshape <- admin1_poly_proj
outfile <- "glRob_50m_admin1_poly"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)
outshape <- admin1_lines_proj
outfile <- "glRob_50m_admin1_lines"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialLinesDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)
outshape <- lakes_poly_proj
outfile <- "glRob_50m_lakes_poly"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)
outshape <- lakes_lines_proj
outfile <- "glRob_50m_lakes_lines"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialLinesDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)
outshape <- lrglakes_poly_proj
outfile <- "glRob_50m_lrglakes_poly"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)
outshape <- lrglakes_lines_proj
outfile <- "glRob_50m_lrglakes_lines"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialLinesDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)
outshape <- grat30_lines_proj
outfile <- "glRob_50m_grat30_lines"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialLinesDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)2.4 Clip out a polygon for the Caspian
Another setup step, and again one that creates some resuable files, is to make a polygon shape file for the Caspian Sea, which can be plotted and filled with color.
The first step is to creat a “bounding box” that surrouds the Caspian.
caspian_bb <- as(extent(45, 56, 35, 50), "SpatialPolygons")
proj4string(caspian_bb) <- unproj_proj4string
summary(caspian_bb)## Object of class SpatialPolygons
## Coordinates:
## min max
## x 45 56
## y 35 50
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]# get the points that define the ouline of the Caspian
caspian_poly <- gIntersection(ocean_poly, caspian_bb)
proj4string(caspian_poly) <- unproj_proj4string
summary(caspian_poly)## Object of class SpatialPolygons
## Coordinates:
## min max
## x 46.70723 54.72324
## y 36.61450 47.11016
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Project and plot the Caspian polygon:
caspian_poly_proj <- spTransform(caspian_poly, robin_crs)
plot(bb_poly_proj, bor="black")
plot(coast_lines_proj, col="black", add=TRUE)
plot(caspian_poly_proj, col="red", add=TRUE)
Write out the Caspian polygon.
outshape <- caspian_poly_proj
outfile <- "glRob_50m_caspian_poly"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)## Warning in writeOGR(outshape, outshapefile, outfile, driver = "ESRI Shapefile", : Field names abbreviated
## for ESRI Shapefile driver3 Set up tree-cover data
One version of Fig. 1 in Marlon et al. (2016) includes a gray-shade “layer” of the “UMD” tree-cover data (Defries et al., 2000, A new global 1-km dataset of percentage tree cover derived from remote sensing, Global Change Biology 6:247-254) [http://glcf.umd.edu/data/treecover/data.shtml], used to provide a background context for the location of the GCDv3 sites. The data were converted to a netCDF file, and this is read, and in turn coverted to spatial points and spatial polygon dataframes.
3.1 Read the tree-cover data
# read a single-variable netCDF dataset using raster()
tree_path <- "/Users/bartlein/Projects/ESSD/data/nc_files/"
tree_name <- "treecov.nc"
tree_file <- paste(tree_path, tree_name, sep="")
tree <- raster(tree_file) # open treecov.nc
tree## class : RasterLayer
## dimensions : 360, 720, 259200 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : /Users/bartlein/Projects/ESSD/data/nc_files/treecov.nc
## names : treecov
## zvar : treecov# levelplot
mapTheme <- rasterTheme(region=brewer.pal(8,"Greens"))
plt <- levelplot(tree, margin=F, par.settings=mapTheme)
plt + latticeExtra::layer(sp.lines(coast_lines, col="black", lwd=0.5))
Create a SpatialPoints data frame, and plot it:
## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## x -179.75 179.75
## y -55.75 83.25
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0]
## Number of points: 60480
## Data attributes:
## treecov
## Min. :-2.00
## 1st Qu.:-1.00
## Median :-1.00
## Mean :20.22
## 3rd Qu.:42.00
## Max. :80.00# plot treecover as spatialpoints
plotclr <- brewer.pal(9,"Greens")
cutpts <- c(10,20,30,40,50,60,70,80,90)
color_class <- findInterval(tree_pts$treecov, cutpts)+1
plot(tree_pts$x, tree_pts$y, col=plotclr[color_class], pch=16, cex=0.25)
plot(coast_lines, add=TRUE)
Convert the points to polygons.
# spatial polygons (takes a while, especially to plot)
tree_poly <- as(tree, "SpatialPolygonsDataFrame")
summary(tree_poly)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -180 180.0
## y -56 83.5
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## treecov
## Min. :-2.00
## 1st Qu.:-1.00
## Median :-1.00
## Mean :20.22
## 3rd Qu.:42.00
## Max. :80.00Plot the treecover polygons.
# plot the polygons
#pdf(file="treecov_poly.pdf")
plotclr <- brewer.pal(9,"Greens")
cutpts <- c(10,20,30,40,50,60,70,80,90)
color_class <- findInterval(tree_poly$treecov, cutpts)+1
clr <- plotclr[color_class]
plot(bb_lines, col="black")
plot(tree_poly, col=clr, bor=clr, lwd=0.1, add=TRUE)
plot(coast_lines, lwd=0.25, add=TRUE)
Write out the shapefiles.
# write tree_ptws
outpath <- "/Users/bartlein/Projects/ESSD/data/RMaps/derived/treecov/"
outshape <- tree_pts
outfile <- "treecov_pts"
outshapefile <- paste(outpath,outfile,sep="")
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)
# write <- tree_poly
outshape <- tree_poly
outfile <- "treecov_poly"
outshapefile <- paste(outpath,outfile,sep="")
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)3.2 Project the treecover points and polygons
Project and plot the shapefiles
## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## x -14606972 16782196
## y -5913980 8319388
## Is projected: TRUE
## proj4string : [+proj=robin +lon_0=0w +ellps=WGS84]
## Number of points: 60480
## Data attributes:
## treecov
## Min. :-2.00
## 1st Qu.:-1.00
## Median :-1.00
## Mean :20.22
## 3rd Qu.:42.00
## Max. :80.00plotclr <- brewer.pal(9,"Greens")
cutpts <- c(10,20,30,40,50,60,70,80,90)
color_class <- findInterval(tree_pts_proj$treecov, cutpts)+1
plot(bb_lines_proj, col="black", axes=FALSE, xlab = "", ylab = "")
points(tree_pts_proj$x, tree_pts_proj$y, col=plotclr[color_class], pch=16, cex=0.25)
plot(coast_lines_proj, add=TRUE)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -14636489 16811803
## y -5939113 8334404
## Is projected: TRUE
## proj4string : [+proj=robin +lon_0=0w +ellps=WGS84]
## Data attributes:
## treecov
## Min. :-2.00
## 1st Qu.:-1.00
## Median :-1.00
## Mean :20.22
## 3rd Qu.:42.00
## Max. :80.00plotclr <- brewer.pal(9,"Greens")
cutpts <- c(10,20,30,40,50,60,70,80,90)
color_class <- findInterval(tree_poly_proj$treecov, cutpts)+1
clr <- plotclr[color_class]
plot(bb_lines_proj, col="black", axes=FALSE, xlab = "", ylab = "")
plot(tree_poly_proj, col=clr, bor=clr, lwd=0.1, add=TRUE)
plot(coast_lines_proj, lwd=0.25, add=TRUE)
Write out the projected shapefiles.
# write tree_pts_proj
outshape <- tree_pts_proj
outfile <- "treecov_pts_proj"
outshapefile <- paste(outpath,outfile,sep="")
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)
# write <- tree_poly_proj
outshape <- tree_poly_proj
outfile <- "treecov_poly_proj"
outshapefile <- paste(outpath,outfile,sep="")
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)4 Map the GCDv3 charcoal records
Figure 1 of Marlon et al. (2016) shows the distribution of charcoal records, and additionally shows by means of symbol color the number of samples in each record. Two versions of the figure will be produced: 1) as published, and 2) with the tree-cover data as a gray-shaded background.
4.1 Read the projected shape files
First read the projected (Robinson) shape files produced above.
# read the projected Robinson shapefiles
shapepath <- "/Users/bartlein/Projects/ESSD/data/RMaps/derived/glRob_50m/"
coast_lines_proj <- readOGR(paste(shapepath, "glRob_50m_coast_lines", sep=""))
admin0_lines_proj <- readOGR(paste(shapepath, "glRob_50m_admin0_lines", sep=""))
admin0_poly_proj <- readOGR(paste(shapepath, "glRob_50m_admin0_poly", sep=""))
bb_lines_proj <- readOGR(paste(shapepath, "glRob_50m_bb_lines", sep=""))
bb_poly_proj <- readOGR(paste(shapepath, "glRob_50m_bb_poly", sep=""))
grat30_lines_proj <- readOGR(paste(shapepath, "glRob_50m_grat30_lines", sep=""))
lrglakes_poly_proj <- readOGR(paste(shapepath, "glRob_50m_lrglakes_poly", sep=""))
caspian_poly_proj <- readOGR(paste(shapepath, "glRob_50m_caspian_poly", sep=""))Read the projected tree-cover polygons (this can be slow–lots of polygons!)
# tree-cover shape file
treecover_path <- "/Users/bartlein/Projects/ESSD/data/RMaps/derived/treecov/"
tree_poly_proj <- readOGR(paste(treecover_path, "treecov_poly_proj.shp", sep="")) ## OGR data source with driver: ESRI Shapefile
## Source: "/Users/bartlein/Projects/ESSD/data/RMaps/derived/treecov/treecov_poly_proj.shp", layer: "treecov_poly_proj"
## with 60480 features
## It has 3 fieldsSet cutpoints and gray-scale color numbers for plotting the tree-cover data. These were calculated to follow a “Munsell sequence” of equally perceived class intervals.
# colors for tree cover
tree_clr_upper <- c( 1, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 999)
tree_clr_gray <- c(1.000, 0.979, 0.954, 0.926, 0.894, 0.859, 0.820, 0.778, 0.733, 0.684, 0.632, 0.576, 0.517, 0.455, 0.389)
colnum <- findInterval(tree_poly_proj$treecov, tree_clr_upper)+1
clr <- gray(tree_clr_gray[colnum], alpha=NULL)4.2 Read the GCDv3 data
Read the GCDv3 charcoal site-location data.
# read the data
csvpath <- "/Users/bartlein/Projects/ESSD/data/csv_files/"
csvname <- "GCDv3_MapData_Fig1.csv"
gcdv3 <- read.csv(paste(csvpath, csvname, sep=""))
gcdv3_pts <- data.frame(cbind(gcdv3$Long,gcdv3$Lat,gcdv3$samples22k))
names(gcdv3_pts) <- c("lon","lat","samples22k")
head(gcdv3_pts)## lon lat samples22k
## 1 -110.6161 44.6628 668
## 2 -110.3467 44.9183 567
## 3 -114.9867 45.7044 367
## 4 -114.2619 45.8919 387
## 5 -114.6503 46.3206 435
## 6 -113.4403 45.8406 971coordinates(gcdv3_pts) <- ~lon+lat
proj4string(gcdv3_pts) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
summary(gcdv3_pts)## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## lon -179.56000 179.5000
## lat -54.88333 69.3833
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0]
## Number of points: 736
## Data attributes:
## samples22k
## Min. : 5.0
## 1st Qu.: 29.0
## Median : 58.0
## Mean : 162.0
## 3rd Qu.: 118.5
## Max. :8166.0Set symbol size and colors for the charcoal points.
# symbol size and colors
nsamp_cutpts <- c(10,100,1000,10000)
nsamp_colors <- c("green3","deepskyblue2","slateblue4","purple")
nsamp_cex <- c(0.5,0.5,0.5,0.5)
nsamp_num <- findInterval(gcdv3_pts$samples22k, nsamp_cutpts)+1Project the charcoal point data.
4.3 Plot the maps
First, plot the sites with no background. The code below creates a .pdf file in the working directory.
# version with no background
pdffile <- "gcdv3_nsamp.pdf"
pdf(paste(pdffile,sep=""), paper="letter", width=8, height=8)
plot(bb_poly_proj, col="gray90", bor="black", lwd=0.1)
plot(grat30_lines_proj, col="black", lwd=0.3, add=TRUE)
plot(admin0_poly_proj, col="white", bor="gray50", lwd=0.4, add=TRUE)
plot(lrglakes_poly_proj, col="gray90", bor="black", lwd=0.2, add=TRUE)
plot(caspian_poly_proj, col="gray90", bor="black", lwd=0.2, add=TRUE)
plot(coast_lines_proj, col="black", lwd=0.5, add=TRUE)
plot(bb_lines_proj, col="black", lwd=1.0, add=TRUE)
plot(gcdv3_pts.proj, pch=2, col=nsamp_colors[nsamp_num], cex=nsamp_cex[nsamp_num], lwd=0.6, add=TRUE)
text(-17000000, 9100000, pos=4, cex=0.8, "GCDv3 -- Number of Samples (Since 22 ka)")
legend(-17000000, -5000000, legend=c("< 10","10 - 100","100 - 1000","> 1000"), bg="white",
title="Number of Samples", pch=2, pt.lwd=0.6, col=nsamp_colors, cex=0.5)
dev.off()The resulting plot will look like this:
# version with treecover background
pdffile <- "gcdv3_nsamp_treecov.pdf"
pdf(pdffile, paper="letter", width=8, height=8)
plot(bb_poly_proj, col="aliceblue", bor="black", lwd=0.1)
plot(grat30_lines_proj, col="black", lwd=0.2, add=TRUE)
plot(tree_poly_proj, col=clr, bor=clr, lwd=0.01, ljoin="bevel", add=TRUE)
plot(admin0_lines_proj, col="gray50", lwd=0.2, add=TRUE)
plot(lrglakes_poly_proj, col="aliceblue", bor="black", lwd=0.2, add=TRUE)
plot(caspian_poly_proj, col="aliceblue", bor="black", lwd=0.2, add=TRUE)
plot(coast_lines_proj, col="black", lwd=0.2, add=TRUE)
plot(bb_lines_proj, col="black", lwd=1.0, add=TRUE)
# print the symbols in white to knockout underlying tree cover
plot(gcdv3_pts.proj, pch=2, col="white", cex=nsamp_cex[nsamp_num], lwd=1.5, add=TRUE)
plot(gcdv3_pts.proj, pch=2, col=nsamp_colors[nsamp_num], cex=nsamp_cex[nsamp_num], lwd=0.6, add=TRUE)
text(-17000000, 9100000, pos=4, cex=0.8, "GCDv3 -- Number of Samples (Since 22 ka)")
legend(-17000000, -5000000, legend=c("< 10","10 - 100","100 - 1000","> 1000"), bg="white",
title="Number of Samples", pch=2, pt.lwd=0.6, col=nsamp_colors, cex=0.5)
dev.off()The resulting plot will look like this: 