Making Maps in GGPlot2

by Mark R Payne
DTU-Aqua, Charlottenlund, Denmark
http://www.staff.dtu.dk/mpay

Mon Dec 14 12:43:15 2015

This tutorial gives a really quick overview of some of the basics of making good looking and useful maps in GGPlot2. It assumes that the reader is familar with the basics of GGPlot2. Note that it also assumes that you are using at least ggplot2 v1.1.0 - at the time of writing, this version was still in development, but can be install from the development repository, following the instructions given at the bottom of the page, here: https://github.com/hadley/ggplot2

The mapping funcctions in ggplot2 are pretty good, and are getting better all the time. Like most things in ggplot2, they are extremely productive, and you can make really good plots for exploring and visualising data really quickly. In this regard it is much better than either lattice, sp or base graphics, which are the main contendors. However, if you really want superstar maps for your publications, I suggest that you may end up going to GMT or base graphics, which still give you the finest level of control.

Note that here I am focusing on using ggplot2 for oceanographic data, which are typically either point data or raster data. If you are dealing with area data you might want to look elsewhere, particularly in the chloropleth community which has lots of good examples, and in the help files of ggplot2.

Setup

Firstly, lets setup up some fake data and make a map of it using the normal (base graphics) plotting routines. For this, you’ll need:

library(maps); 
library(mapdata)

Generate some synthetic data. Note that I’m not going to get too stressed about the validity of this data - just some random points, some of which may be on land and some in the ocean.

n <- 50
dat <- data.frame(lon=runif(n,7,10),
                  lat=runif(n,55,60))

Plot it and overlay the coastlines

plot(dat,pch=16,col="red")
map("worldHires",add=TRUE,fill=TRUE)

which is a bit ugly really, even if its technically correct. Lets make the map first, and then overlay the points

xlims <- range(pretty(dat$lon))
ylims <- range(pretty(dat$lat))
map("worldHires",xlim=xlims,ylim=ylims,fill=TRUE,mar=c(5,2,0,0))
points(dat,pch=16,col="red")
map.axes()

Note that map() takes care of the aspect ratio of the plot so that it actually looks like a map. You can achieve a similar result by specifying the asp= argument like so

plot(dat,pch=16,col="red",asp=2)
map("worldHires",add=TRUE,fill=TRUE)

There are numerous ways to arrive at the “right” asp - guessing is a good one. Atleratnive, you can plot a blank map first, then your points, and then add the map after that. You can also use the mapasp() function in the sp package.

GGPlot2

And now to GGPlot2. Here’s a starter:

library(ggplot2)
ggplot(dat,aes(x=lon,y=lat))+geom_point(col="red",pch=16)

..which is not so useful. We can add a coastline layer through the borders() functionality

ggplot(dat,aes(lon,lat))+borders() + geom_point(col="red",pch=16)

Whether you should include political borders is a case of ‘smag og behag’, ie its a matter of taste. Personally, I hate them, and usually just make all land black (but that’s because I’m interested in the ocean)

ggplot(dat,aes(lon,lat))+borders(fill="black",colour="black") +
  geom_point(col="red",pch=16)

This map is technically correct, but visually useless. We’d like it a little more focused. What GGPlot2 is doing here is taking the polygons from the maps and mapdata packages and plotting them as a layer on your plot. We can zoom in to the region of interest using the coord_fixed() option

ggplot(dat,aes(lon,lat))+borders(fill="black",colour="black") +
  geom_point(col="red",pch=16) +
  coord_fixed(xlim = xlims,ylim=ylims)

which looks better. Note that if you us xlim() and ylim() to set the coordinates, you will break your polygons like so

ggplot(dat,aes(lon,lat))+borders(fill="black",colour="black") +
  geom_point(col="red",pch=16) +
  xlim(xlims) + ylim(ylims)

## Warning: Removed 101812 rows containing missing values (geom_polygon).

… because xlim() and ylim() throw out all points that are outside of the bounds. Generally it’s better ot use coord_fixed() instead. However, this plotting arrangement is actually plotting the entire global database from the maps package, which can get pretty slow (especially if you use the “worldHires” map). You can confine the size of the polygon database that is extracted using the xlim and ylim arguments to borders e.g.

ggplot(dat,aes(lon,lat))+
  borders(fill="black",colour="black",xlim=xlims,ylim=ylims) +
  geom_point(col="red",pch=16) 

Note that this has only taken polygons that intersect our region of interest (ROI) rather than the entire global database. However, its still usually necessary to restrict the coordinates down again with coord_fixed() (this is a significant sourcec of irration, but from my understanding of how ggplot2 works at the lowest level, is essentially unavoidable)

ggplot(dat,aes(lon,lat))+
  borders(fill="black",colour="black",xlim=xlims,ylim=ylims) +
  geom_point(col="red",pch=16) +
  coord_fixed(xlim = xlims,ylim=ylims)

The default “projection” here sets a 1:1 ratio, which is not appropriate at non-tropical latitudes. The quickest way to fix this is using the quick_map() function to specifiy the coordinate limits

ggplot(dat,aes(lon,lat))+
  borders(fill="black",colour="black",xlim=xlims,ylim=ylims) +
  geom_point(col="red",pch=16) +
  coord_quickmap(xlim = xlims,ylim=ylims)

which is better. You can do a whole lot of stuff with projections here with the coord_map() function, which I’m not going to explore too much here. But just to prove the point

ggplot(dat,aes(lon,lat))+
  borders(fill="black",colour="black",xlim=xlims,ylim=ylims) +
  geom_point(col="red",pch=16) +
  coord_map("mercator")

I’m not sure how to specify both the xlim and ylim coordinates and the projection.

You can use this approach in conjunction the full plotting power ggplot2, which is really useful. e.g.

dat$facet <- sample(letters[1:3],n,replace=TRUE)
g.base <- ggplot(dat,aes(lon,lat))+
        borders(fill="black",colour="black",xlim=xlims,ylim=ylims) +
        coord_quickmap(xlim = xlims,ylim=ylims)
print(g.base)

g1 <- g.base + geom_point(aes(col=facet,shape=facet),size=2)
print(g1)

g2 <- g1 + facet_wrap(~facet)
print(g2)

Finally, if you want better quality coastlines or other maps, use the worldhires map database in borders.e.g.

ggplot()+
  borders("worldHires",fill="black",colour="black",xlim=xlims,ylim=ylims) +
  coord_quickmap(xlim = xlims,ylim=ylims)

ggplot() +borders("nz",fill="black",colour="black")+coord_quickmap()

See the help files in the maps and mapdata packages for more details about available databases.

Working with rasters

Finally, a brief example of working with raster data in GGPlot2. To illustrate the point, we’ll use data from HadISST.

library(raster)
b.all <- brick("data/HadISST_sst.nc")
b <- b.all[[c(1,6)]]
b[b < -300] <- NA

So, firstly with the default functions in raster

plot(b)

And now with ggplot2. There is a not a particularly straightforward way to do this - first we have to extract the data from the data frame, then add the coordinates in manually.

b.dat <- data.frame(coordinates(b),as.data.frame(b))
head(b.dat)

##        x    y X1870.01.16 X1870.06.16
## 1 -179.5 89.5          NA          NA
## 2 -178.5 89.5          NA          NA
## 3 -177.5 89.5          NA          NA
## 4 -176.5 89.5          NA          NA
## 5 -175.5 89.5          NA          NA
## 6 -174.5 89.5          NA          NA

Next we need to reshape it into the long-format used by ggplot2 using the “melt” function from the reshape package

library(reshape)
library(maptools)
plt.dat <- melt(b.dat,id.vars = c("x","y"))
head(plt.dat)

##        x    y    variable value
## 1 -179.5 89.5 X1870.01.16    NA
## 2 -178.5 89.5 X1870.01.16    NA
## 3 -177.5 89.5 X1870.01.16    NA
## 4 -176.5 89.5 X1870.01.16    NA
## 5 -175.5 89.5 X1870.01.16    NA
## 6 -174.5 89.5 X1870.01.16    NA

And lets roll

ggplot(plt.dat,aes(x=x,y=y,fill=value))+
  geom_raster()+borders(fill="black",colour="black") + facet_wrap(~variable)+
  scale_fill_distiller(palette="RdPu")+ coord_quickmap()

Note that you should be careful about the size of the data frame that you end up creating here - it can be very easy to make massive objects when you start playing with rasters. One strategy is to use the support for ggplot in the rasterVis package which achieves a lot of the same stuff directly e.g.

library(rasterVis)
gplot(b)+geom_raster(aes(fill=value))+
  borders(fill="black",colour="black") + facet_wrap(~variable)+
  scale_fill_distiller(palette="RdPu")+ coord_quickmap()

This work by Mark R Payne is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License. For details, see http://creativecommons.org/licenses/by-nc-sa/3.0/deed.en_US Basically, this means that you are free to “share” and “remix” for non-commerical purposes as you see fit, so long as you “attribute” me for my contribution. Derivatives can be distributed under the same or similar license.

This work comes with ABSOLUTELY NO WARRANTY or support.

This work should also be considered as BEER-WARE. For details, see http://en.wikipedia.org/wiki/Beerware