Natural Earth Data and R in ggplot2

I really like the data layers available from naturalearthdata.  They take a simple map and make it look beautiful, which I suppose is the point.

Figure 1.  Here's an amazing background for your maps using Natural Earth Data and ggplot2
Figure 1. Here’s an amazing background for your maps using Natural Earth Data and ggplot2

One thing that I’ve found frustrating though is the need to constantly go back and forth between R and ArcGIS/QGIS/GRASS to make figures with spatial data if I want them to look really pretty.  In preparation for AGU2013 this year (come see our poster) I decided I’d finally solve this annoying problem by figuring out how to get a pretty base raster layer into ggplot, and then superimpose the lakes and streams that serve to make a nice map look pretty.

Here’s the solution to make this fine figure:


#  Assuming you have a path 'Maps' that you store your spatial files in.  This
#  is all downloaded from <a href=""></a> using the
#  1:50m "Medium" scale data. <- stack('./Maps/NE2_50M_SR_W/NE2_50M_SR_W.tif')

ne_lakes <- readOGR('./Maps/NaturalEarth/ne_50m_lakes.shp',

ne_rivers <- readOGR('./Maps/NaturalEarth/ne_50m_rivers_lake_centerlines.shp',

ne_coast <- readOGR('./Maps/NaturalEarth/ne_50m_coastline.shp',

#  I have a domain I'm interested in, but there's no reason you can't define something else:
quick.subset <- function(x, longlat){

# longlat should be a vector of four values: c(xmin, xmax, ymin, ymax)
  x@data$id <- rownames(x@data)

  x.f = fortify(x, region="id")
  x.join = join(x.f, x@data, by="id")

  x.subset <- subset(x.join, x.join$long > longlat[1] & x.join$long < longlat[2] &
                           x.join$lat > longlat[3] & x.join$lat < longlat[4])


domain <- c(-98.6, -66.1, 36.5, 49.7)
lakes.subset <- quick.subset(ne_lakes, domain)
river.subset <- quick.subset(ne_rivers, domain)
coast.subset <- quick.subset(ne_coast, domain)

nat.crop <- crop(, y=extent(domain))

rast.table <- data.frame(xyFromCell(nat.crop, 1:ncell(nat.crop)),

rast.table$rgb <- with(rast.table, rgb(NE2_50M_SR_W.1,
# et voila!

ggplot(data = rast.table, aes(x = x, y = y)) +
  geom_tile(fill = rast.table$rgb) +
  geom_polygon(data=lakes.subset, aes(x = long, y = lat, group = group), fill = '#ADD8E6') +
  scale_alpha_discrete(range=c(1,0)) +
  geom_path(data=lakes.subset, aes(x = long, y = lat, group = group), color = 'blue') +
  geom_path(data=river.subset, aes(x = long, y = lat, group = group), color = 'blue') +
  geom_path(data=coast.subset, aes(x = long, y = lat, group = group), color = 'blue') +
  scale_x_continuous(expand=c(0,0)) +
  scale_y_continuous(expand=c(0,0)) +
  xlab('') + ylab('')

Some issues still remain for me. For one, there’s a problem with holes in this plotting. You’ll notice that Manitoulin Island and half of the Keweenaw Peninsula are basically lake (the light blue background) outlined in blue, and not the textured raster. I’m not sure exactly how to fix that, there’s plenty of discussion about how to fix holes in polygons using ggplot, but I haven’t worked out the solution yet.


If you have a solution I’m all ears.

Published by


Assistant scientist in the Department of Geography at the University of Wisconsin, Madison. Studying paleoecology and the challenges of large data synthesis.

9 thoughts on “Natural Earth Data and R in ggplot2”

  1. I am really excited about this post! But it’s not working for me. I also got my data from Natural Earth: no .tif. I went back & looked by following your link & there isn’t a .tif file! Just your typical dbf, prj, shp, shx. I just need a simple map of the US cut off at the 102 meridian, so I figured I could just read in the ne_50m_land.shp file. I got as far as”
    rast.table <- data.frame …
    And it gives me error:
    Error in (function (classes, fdef, mtable) :
    unable to find an inherited method for function ‘xmin’ for signature ‘"SpatialPolygonsDataFrame"’
    Even with your exact code (except .shp , no .tif, but I don't think that matters at this point).
    I know this isn't StackOverload, but I wanted to let you know. There's a small chance your code needs to be updated, though I, R newbie, probably unwittingly messed it up somewhere.

    1. If you’re not going to use a raster base layer then you don’t need the stack command (where I reference the ‘tif’ file). That’s why the ‘rast.table’ command isn’t working, because it’s looking for a raster stack. When I look here:

      I see links to the raster files, maybe you’re at too high a scale?

      My email address is linked in the ‘About’ section, feel free to send me the code and we can work it out together.

  2. Hi, I am also very excited about making such nice maps in R using the natural earth data.

    Unfortunately I am also an R newbie and am having quite some difficulty working around the .tif problem previously mentioned. I am working at the highest resolution, which I need for the maps I want to create, could you perhaps explain in some more detail how to get around this?

    Thanks for your patience

      1. If you could show what code would be needed for using the large scale data (10m) that would be really helpful. I have to use that resolution and I’ve tried tweaking your code (using the population places, the urban, and the roads downloads), but am getting stuck.

  3. Nice post. This helped me make a nice map. Instead of subsetting the shape files I think it’d be easier to use coord_cartesian(xlim=…,ylim=..) and just zoom into the area you want!
    Also, one word of caution for others, I tried to use the maxColorValue argument from rgb() but it doesn’t give you the same hex values as in this example (where you divide by 255 to normalize to 1). Pretty easy to tell which one is correct!

    1. Thanks! I’ve actually been using coord_cartesian in my own code for about a month now, I thought about editing the post but have never gotten around to it :) I appreciate you mentioning it here.

      The one limitation is that some of the files are very big, so if you’re doing something that’s otherwise fairly data intensive it makes sense to clip everything anyway, but maybe use a bigger buffer so that you don’t get weird clipping artifacts.

Leave a Reply

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

You are commenting using your 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