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:

Now read the shape files, including those for global coastlines, adminstrative units (borders), large lakes, and a graticule. Set the filenames:

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.

## 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.

## 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.)

Plot everything (except ocean).

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.

## CRS arguments:
##  +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## CRS arguments: +proj=robin +lon_0=0w +ellps=WGS84

Now do the projections.

Plot the projected shapefiles.

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.

Plot the SpatialLines objects

2.3 Write out the projected shapefiles

Next, write out the projected shapefiles, first setting the output path.

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.

## 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]
## 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:

Write out the Caspian polygon.

## Warning in writeOGR(outshape, outshapefile, outfile, driver = "ESRI Shapefile", : Field names abbreviated
## for ESRI Shapefile driver

[Back to top]

3 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

## 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

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

Convert the points to polygons.

## 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.00

Plot the treecover polygons.

Write out the shapefiles.

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.00

## 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.00

Write out the projected shapefiles.

[Back to top]

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 tree-cover polygons (this can be slow–lots of polygons!)

## 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 fields

Set 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.

4.2 Read the GCDv3 data

Read the GCDv3 charcoal site-location data.

##         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        971
## 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.0

Set symbol size and colors for the charcoal points.

Project 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.

The resulting plot will look like this:

The resulting plot will look like this:

[Back to top]