Making maps in R

I don’t have a lot of experience with geospatial data analysis but following a recent cruise I had the need to plot some geospatial data.  The go-to program for this seems to be Matlab (or even ArcGIS, for those who are serious about their map making), but I do almost all of my plotting in R.  This sent me on a bit of an exploration of the various map making options for R.  Although the maturity of the various R geospatial libraries is far below that of commercial products, it’s possible to make some pretty decent plots.  Here’s a synopsis of the method I ended up using.

To plot data on a map you first need a map to work from.  A map comes in the form of a shape file, which is nothing more than an collection of lines or shapes describing the outline of various map features.  Shape files come in various flavors, including lines and polygons.  For reasons that will become clear working with polygons is much preferable to working with lines, but I found it difficult (impossible) to consistently get the right polygon format into R.  R does have some built in polygons in the maps and mapdata packages.  Using these it is possible to access a global polygon map, however the resolution of the map is quite low.  The documentation doesn’t provide a maximum resolution, but it is low enough to be of little value on scales below hundreds of miles.  Here’s a quick plot of our study area (the Antarctic Peninsula) with no further information plotted, just to get a sense of the resolution.

library(maps)
library(mapdata)
map(database = 'world',
    regions = "antarctica",
    xlim = c(-77, -55),
    ylim = c(-70, -60),
    fill = T,
    col = 'grey',
    resolution = 0,
    bg = 'white',
    mar = c(1,1,2,1)
)

box()

 Not too bad, but the y dimension is over 1000 km.  If our map covered 10 km or even 100 km we’d have at best a large block to define any landmass.  To produce a more detailed map we need to use a shape file with a much higher resolution.  This is where things get a little tricky.  There are high resolution shape files available from a number of sources, primarily government agencies, but finding them can be difficult.  After a lot of searching I discovered the GSHHG High-resolution Geography Database available from NOAA.  The database is quite large, but it can be subset with the GEODAS Coastline Extractor which also makes use of the WDBII database (included with the GSHHG link) containing rivers and political boundaries.  It takes a little fussing around, but once you get the hang of the coastline extractor it isn’t too bad to work with.  On first use you have to point the extractor to the databases.  Then select “plot” from the drop down menu, and enter the desired lat and long range.  This returns a nice graphic of the selected region which you can export in the desired shape file format using “Save as”.

The export step is where I’m a bit stuck.  Ideally you would export the coastline as a polygon so that R could color it, if desired (as in the previous example).  I can only get R to recognize shape files containing lines however, exported in the ESRI format.  Here’s the region plotted earlier but now from the GSHHG database:

library(maps)
library(mapdata)
library(maptools)
plter_shape <- readShapeLines('plter_region_shore_shore',
                              repair = T)
map("world",
    "antarctica",
    xlim = xlim,
    ylim = ylim,
    fill = T,
    col = 'grey',
    resolution = 0,
    bg = 'white',
    mar = c(6,3,2,1),
    type = 'n'
)

plot(plter_shape, add = T)
box()

Which produces this plot:

That’s quite a bit more detail.  If we zoomed in the improvement would be even more pronounced.  You’ll notice that I plotted the ESRI shape file on a plot produced from maps, even though I had no intention of using maps to reproduce the coastline.  Using maps allows R to plot the data as a projection, which is particularly important at high latitudes.  The default maps projection isn’t great for 60 degrees south (see the maps documentation for a list of available projections), but it gets the point across.  A square projection would be much worse.  What’s particularly nice about all this is that R will use the projected coordinates for all future calls to this plot.  Points or lines added to the plot will thus show up in the right place.

Now that we have our basic map we can add some data.  One of our science tasks was to evaluate bacterial production along the Palmer Long Term Ecological Research Station grid.  Plotting depth-integrated values on our map lets us see where bacterial production was high.  I’m not going to bother showing how I arrived at depth integrated values because the code is too specific to my data format, but suffice to say that I used the interp1 function to interpolate between depth points at each station, then I summed the interpolated values.  This gives me a data frame (depth_int) of three variables: lon, lat, and bp (bacterial production).

First, let’s generate some colors to use for plotting:

z <- depth_int$bp
zcol <- colorRampPalette(c('white', 'blue', 'red'))(100)[as.numeric(cut(z, breaks = 100))]

Then we add the points to the plot.  This is no different than building any other R plot:

points(depth_int$lon,
       depth_int$lat,
       bg = zcol,
       pch = 21,
       cex = 2,
       col = 'grey'
)

Snazz it up a little by adding grid lines and axis labels:

axis(1, labels = T)
axis(2, labels = T)
grid()

Titles…

title <- expression(paste('Bacterial production gC m'^{-2},'d'^{-1}, sep = " "))
title(xlab = 'Longitude',
      ylab = 'Latitude',
      main = 'title')

And the coup de grâce, a scale bar.  I can’t find a good scale bar function, but it’s easy enough to build one from scratch:

## set the location and the colorbar gradation
xleft <- -74
xright <- -73
ybot <- -65
yint <- (65 - 62.5) / 100
ytop <- ybot + yint

## create the bar by stacking a bunch of colored rectangles
for(c in colorRampPalette(c('white', 'blue', 'red'))(100)){
  ybot = ybot + yint
  ytop = ytop + yint
  rect(xleft, ybot, xright, ytop, border = NA, col = c)
  print(c(xleft, xright, ybot, ytop, c))
}

## generate labels
labels <- round(seq(min(z), max(z), length.out = 5),2)

## add the labels to the plot
text(c(xright + 0.2),
     seq(-65, -62.5, length.out = 5),
     labels = as.character(labels),
     cex = 0.8,
     pos = 4)

Done correctly this should approximate the plot published in my earlier post.  Plenty of tweaks one could make (I’d start by shifting the location of the scale bar and drawing an outlined white box around it), but it’s not a bad start!