I heart ggplot2

Since discovering the R package ggplot2, I’ve never looked back. Though I will admit that much of the documentation is a headache, millions have shared their code online. I cannot recommend enough spending a little time getting to know how to use the package, as it will enable you in a few simple steps, to produce very cool graphics, as can be seen below:

ggplot1

Here I will walk through some code for drawing a pacific-centred globe in ggplot, and adding either points or lines. Most default maps are centered around Europe, making Pacific-centered maps difficult to plot as many adjustments need to be made. Here’s how I have gone about dealing with these issues in ggplot.

First and foremost, set your working directory setwd() and upload your packages:

#load packages
library(maps)
library(geosphere)
library(plyr)
library(ggplot2)
library(sp)

Now the package sp has a bug with the fortify function, which you can fix by running this lump of code, awesomely shared from here.

#Deals with problem with fortify

#' Fortify method for classes from the sp package.
#'
#' To figure out the correct variable name for region, inspect
#' /code{as.data.frame(model)}.
#'
#' @param model \code{SpatialPolygonsDataFrame} to convert into a dataframe.
#' @param data not used by this method
#' @param region name of variable used to split up regions
#' @param ... not used by this method
#' @name fortify.sp
#' @examples
#' if (require("maptools")) {
#' sids <- system.file("shapes/sids.shp", package="maptools")
#' nc1 <- readShapePoly(sids,
#' proj4string = CRS("+proj=longlat +datum=NAD27"))
#' nc1_df <- fortify(nc1)
#' }
NULL
#' @rdname fortify.sp
#' @export
#' @method fortify SpatialPolygonsDataFrame
fortify.SpatialPolygonsDataFrame <- function(model, data, region = NULL, ...) {
 attr <- as.data.frame(model)
 # If not specified, split into regions based on polygons
 if (is.null(region)) {
 coords <- ldply(model@polygons,fortify)
 message("Regions defined for each Polygons")
 } else {
 cp <- sp::polygons(model)
 try_require("maptools")
 # Union together all polygons that make up a region
 unioned <- unionSpatialPolygons(cp, attr[, region])
 coords <- fortify(unioned)
 coords$order <- 1:nrow(coords)
 }
 coords
}
#' @rdname fortify.sp
#' @export
#' @method fortify SpatialPolygons
fortify.SpatialPolygons <- function(model, data, ...) {
 ldply(model@polygons, fortify)
}
#' @rdname fortify.sp
#' @export
#' @method fortify Polygons
fortify.Polygons <- function(model, data, ...) {
 subpolys <- model@Polygons
 pieces <- ldply(seq_along(subpolys), function(i) {
 df <- fortify(subpolys[[model@plotOrder[i]]])
 df$piece <- i
 df
 })
 within(pieces,{
 order <- 1:nrow(pieces)
 id <- model@ID
 piece <- factor(piece)
 group <- interaction(id, piece)
 })
}
#' @rdname fortify.sp
#' @export
#' @method fortify Polygon
fortify.Polygon <- function(model, data, ...) {
 df <- as.data.frame(model@coords)
 names(df) <- c("long", "lat")
 df$order <- 1:nrow(df)
 df$hole <- model@hole
 df
}
#' @rdname fortify.sp
#' @export
#' @method fortify SpatialLinesDataFrame
fortify.SpatialLinesDataFrame <- function(model, data, ...) {
 ldply(model@lines, fortify)
}
#' @rdname fortify.sp
#' @export
#' @method fortify Lines
fortify.Lines <- function(model, data, ...) {
 lines <- model@Lines
 pieces <- ldply(seq_along(lines), function(i) {
 df <- fortify(lines[[i]])
 df$piece <- i
 df
 })
 within(pieces,{
 order <- 1:nrow(pieces)
 id <- model@ID
 piece <- factor(piece)
 group <- interaction(id, piece)
 })
}
#' @rdname fortify.sp
#' @export
#' @method fortify Line
fortify.Line <- function(model, data, ...) {
 df <- as.data.frame(model@coords)
 names(df) <- c("long", "lat")
 df$order <- 1:nrow(df)
 df
}

Once it has been run, you no longer need to run it again. Now one problem with re-centring a map around the Pacific Ocean is that it mucks up your longitudes, which you also need to re-centre. You also have the problem using the default worldmap in the “maps” package, that it cuts your lovely polygons and does weird things as can be seen here. The bit of code below deals with this.

# Function to regroup split lines and polygons
# takes dataframe, column with long and unique group variable,
#returns df with #added column named group.regroup
RegroupElements <- function(df, longcol, idcol){
 g <- rep(1, length(df[,longcol])) if (diff(range(df[,longcol])) > 300) {
 d <- df[,longcol] > mean(range(df[,longcol]))
 g[!d] <- 1
 g[d] <- 2
 }
 g <- paste(df[, idcol], g, sep=".")
 df$group.regroup <- g
 df
}

# Function to close regrouped polygons
# takes dataframe, checks if 1st and last longitude value are the same,
#if not, inserts first as last and reassigns order variable
ClosePolygons <- function(df, longcol, ordercol){
 if (df[1,longcol] != df[nrow(df),longcol]) {
 tmp <- df[1,]
 df <- rbind(df,tmp)
 }
 o <- c(1: nrow(df)) # rassign the order variable
 df[,ordercol] <- o
 df
}

For this example, I use tracking data from this paper by Driscoll and Ueta The migration route and behaviour of Eastern Curlews Numenius madagascariensis for my maps. The data is of a tracked bird and in the paper is the individual marked ID16091.

#Load data
speciestrax<-matrix(0,11,5)
colnames(speciestrax)<-c("slong","slat","elong","elat","number")
speciestrax[,1]<-c(153.1590,134.4637,120.7381,119.6376,120.2324,130.1927,130.2522,124.3630,121.7456,120.6749,141.5546)
speciestrax[,2]<-c(-27.178338,7.485653,18.468052,26.189403,36.111755,46.200918,49.591650,46.914750,40.906617,33.232855,-2.815979)
speciestrax[,3]<-c(134.4637,120.7381,119.6376,120.2324,130.1927,130.2522,124.3630,121.7456,120.6749,141.5546,143.9936)
speciestrax[,4]<-c(7.485653,18.468052,26.189403,36.111755,46.200918,49.591650,46.914750,40.906617,33.232855,-2.815979,-13.939959)
speciestrax[,5]<-1:11
rownames(speciestrax)<-c("stage1","stage2","stage3","stage4","stage5","stage6","stage7","stage8","stage9","stage10","stage11")

Now that we have everything that we need, we can start plotting!

Recentering points

In my first example, I map points (stopover sites for birds) which therefore need to be recentered around the pacific. This is achieved in the following section of code which determines what longitude we want to be the centre, and then recalculates the longitudes of our points.

#Plotting stuff
center=160 #describe where you want the centre of the map to be

#data<-speciestrax
data<-data.frame(as.matrix(speciestrax))

# shift coordinates to recenter great circles
data$slong.recenter <- ifelse(data$slong < center - 180 , data$slong + 360, data$slong)
data$elong.recenter <- ifelse(data$elong < center - 180 , data$elong + 360, data$elong)

# shift coordinates to recenter worldmap
worldmap2 <- map_data("world")
worldmap2$long.recenter <- ifelse(worldmap2$long < center - 180 , worldmap2$long + 360, worldmap2$long)

# now regroup
worldmap2.rg <- ddply(worldmap2, .(group), RegroupElements, "long.recenter", "group")

# close polys
worldmap2.cp <- ddply(worldmap2.rg, .(group.regroup), ClosePolygons, "long.recenter", "order") # use the new grouping var

Now for some plotting! Like previously mentioned, I love doing these sorts of maps in ggplot. It enables you to draw lovely globes and add points to them quite simply.

print(ggplot() +
        geom_polygon(aes(long.recenter,lat,group=group.regroup),
                      data=worldmap2.cp) +
        geom_point(aes(slong.recenter,slat,
                       color=number),data=data, size=5,alpha=0.8) +
        theme(panel.background = element_rect(fill = "white"),
              panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
              axis.ticks = element_blank(), axis.title.x = element_blank(),
              axis.title.y = element_blank(), axis.text.x = element_blank(),
              legend.title=element_blank(),axis.text.y = element_blank()   ) +
        coord_map("ortho", orientation=c(10,140,0))
)

globe2

As you can see, the ggplot2 default is blue. The function geom_polygon() draws the countries, geom_point() draws the dots, theme() describes the titles of the axes and legend (you will notice I have set most of these to element_blank() so that there is no title or x/y axis lables) and finally coord_map() is used to determine the x,y and z rotation of the globe. So you can play around with these to rotate the Earth… Pretty neat.

Now we can play around with some of the colours and settings.

print(ggplot() +
        geom_polygon(aes(long.recenter,lat,group=group.regroup),
                     size = 0.2, fill="white", colour = "black", data=worldmap2.cp) +
        geom_point(aes(slong.recenter,slat,
                       color=number,size=number),data=data, alpha=0.8) +
        scale_color_gradient(high="gold" ,low="brown")+
	  scale_size(range=c(7,20))+
        theme(panel.background = element_rect(fill = "black"),
              panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
              axis.ticks = element_blank(), axis.title.x = element_blank(),
              axis.title.y = element_blank(), axis.text.x = element_blank(),
              legend.title=element_blank(),axis.text.y = element_blank(), legend.position = "none") +
        coord_map("ortho", orientation=c(10,140,0))
)

globe3

As you will notice, it’s also possible to change the colours of the globe using fill and colour in geom_polygon(). theme() can be used to change the background colour. You may also notice that the points are coloured and sized according to “number” a random column that I created. This is done via the scale_color_gradient and scale_size functions.

Recentering lines

Now instead of points I will use lines to show the migratory route used by a bird. It’s slightly more complex in terms of the recentering, but with this bit of code it is easy enough.

rts <- gcIntermediate(speciestrax[,1:2],speciestrax[,3:4], 100, breakAtDateLine=FALSE, addStartEnd=TRUE, sp=TRUE)
rts.ff <- fortify.SpatialLinesDataFrame(rts)
trax.ll<-speciestrax
trax.ll$id <-as.character(c(1:nrow(speciestrax))) # that rts.ff$id is a char
gcircles <- merge(rts.ff, trax.ll, all.x=T, by="id") # join attributes, we keep them all, just in case

# shift coordinates to recenter great circles
gcircles$long.recenter <- ifelse(gcircles$long < center - 180 , gcircles$long + 360, gcircles$long)

# shift coordinates to recenter worldmap
worldmap <- map_data ("world")
worldmap$long.recenter <- ifelse(worldmap$long < center - 180 , worldmap$long + 360, worldmap$long)

# now regroup
gcircles.rg <- ddply(gcircles, .(id), RegroupElements, "long.recenter", "id")
worldmap.rg <- ddply(worldmap, .(group), RegroupElements, "long.recenter", "group")

# close polys
worldmap.cp <- ddply(worldmap.rg, .(group.regroup), ClosePolygons, "long.recenter", "order") # use the new grouping var

Now that this has been done we can plot it out

print(ggplot() +
        geom_polygon(aes(long.recenter,lat,group=group.regroup),
                     size = 0.2, fill="grey", colour = "grey", data=worldmap2.cp) +
        geom_line(aes(long.recenter,lat,group=group.regroup,color=id),
	size=1.3, data= gcircles.rg) + # set transparency here
	  scale_size(range=c(7,20))+
        theme(panel.background = element_rect(fill = "white"),
              panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
              axis.ticks = element_blank(), axis.title.x = element_blank(),
              axis.title.y = element_blank(), axis.text.x = element_blank(),
              legend.title=element_blank(),axis.text.y = element_blank(), legend.position = "none") +
        coord_map("ortho", orientation=c(10,140,0))
)

globe4

As you can see the default line colours are quite colourful, and personally I don’t care much for them. You will also note that I’ve changed the colour of the land to grey (just because I can). To make the lines a different colour using a colourpalette and to vary the thickness according to number (the fake column I created in the data) you can use the following bit of code:

print(ggplot() +
        geom_polygon(aes(long.recenter,lat,group=group.regroup),
                     size = 0.2, fill="pink", colour = "black", data=worldmap2.cp) +
        geom_line(aes(long.recenter,lat,group=group.regroup,color=id,size=id),
	data= gcircles.rg,alpha=0.5) + # set transparency here
	scale_colour_brewer(palette="PuOr") +
	scale_size_manual(values=c(3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8)) + #scale_size(to = c(1, 2)) +
	  theme(panel.background = element_rect(fill = "white"),
              panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
              axis.ticks = element_blank(), axis.title.x = element_blank(),
              axis.title.y = element_blank(), axis.text.x = element_blank(),
              legend.title=element_blank(),axis.text.y = element_blank(), legend.position = "none") +
        coord_map("ortho", orientation=c(10,140,0))
)

globe5
I have again changed the colour of the earth for fun and made the lines transparent using alpha in geom_line().

There is of course plenty more you can do, this is a simple starting point. Enjoy!

Advertisements

2 thoughts on “I heart ggplot2

    1. Not that I know of I’m afraid. But there seem to be some more elegant work arounds for the polygon splicing at http://stackoverflow.com/questions/5353184/fixing-maps-library-data-for-pacific-centred-0%C2%B0-360%C2%B0-longitude-display.
      Specifically this function:

      plot.map<- function(database,center,…){
      Obj <- map(database,…,plot=F)
      coord <- cbind(Obj[[1]],Obj[[2]])

      # split up the coordinates
      id <- rle(!is.na(coord[,1]))
      id <- matrix(c(1,cumsum(id$lengths)),ncol=2,byrow=T)
      polygons <- apply(id,1,function(i){coord[i[1]:i[2],]})

      # split up polygons that differ too much
      polygons <- lapply(polygons,function(x){
      x[,1] <- x[,1] + center
      x[,1] 180,x[,1]-360,x[,1])
      if(sum(diff(x[,1])>300,na.rm=T) >0){
      id <- x[,1] < 0
      x <- rbind(x[id,],c(NA,NA),x[!id,])
      }
      x
      })
      # reconstruct the object
      polygons <- do.call(rbind,polygons)
      Obj[[1]] <- polygons[,1]
      Obj[[2]] <- polygons[,2]

      map(Obj,…)
      }

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s