Monthly Archives: April 2014

Finally, an easy way to fix the horizontal lines in ggplot2 maps

ggplot2 tries to make it super easy to add country or state borders to your map, and for the most part it works great as long as you include the entire map region in your plot (all the states or the entire world map for example). A long standing issue with the ggplot2 borders() function is when you try and zoom in a way that removes data outside of the plotting region, countries and states get cut off and not reconnected correctly, resulting in spurious horizontal lines in the plot. A forum post and ggplot2 issue both allude to this problem.

Here is a simple example showing this issue. First, a simple unprojected map that looks great because we are not trying to zoom in at all.

ggplot() +
    borders("world") + borders("state") +
    theme_bw()

Alt text

An obvious extension would be to zoom in on a region and project the map:

ggplot() +
    borders("world") + borders("state") +
    theme_bw() +
    coord_map(xlim=xlim,ylim=ylim) +
    labs(y="",x="")

Alt text

Which pretty clearly illustrates the issue. Things look even worse if we try and add color:

ggplot() +
    borders("world",fill="darkseagreen") +
    borders("state",fill="darkseagreen") +
    theme_bw() +
    coord_map(xlim=xlim,ylim=ylim) +
    labs(y="",x="")

Alt text

The underlying issue is that when we zoom with coord_map() the polygons get cut up, with some parts appearing on either side of the map, hence the horizontal lines when ggplot tries to connect a country or state that appears on two sides of the map. The solution to this issue is to get the raw map data then use one of the various polygon clipping libraries in R. I used PBSmapping. Here is a complete example

library("PBSmapping")
library("ggplot2")
library("maps")
library("data.table")

# plot limits
xlim = c(-165,-90)
ylim = c(10,60)

worldmap = map_data("world")
setnames(worldmap, c("X","Y","PID","POS","region","subregion"))
worldmap = clipPolys(worldmap, xlim=xlim,ylim=ylim, keepExtra=TRUE)

statemap = map_data("state")
setnames(statemap, c("X","Y","PID","POS","region","subregion"))
statemap = clipPolys(statemap, xlim=xlim,ylim=ylim, keepExtra=TRUE)

p = ggplot() +
    coord_map(xlim=xlim,ylim=ylim) +
    geom_polygon(data=worldmap,aes(X,Y,group=PID),
        fill = "darkseagreen",color="grey50") +
    geom_polygon(data=statemap,aes(X,Y,group=PID),
        fill = "darkseagreen",color="grey50") +
    labs(y="",x="") +
    theme_bw()
 print(p)

Alt text

This could also be used to fix the fill issues when doing ggplot filled contour plots.