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.

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:

library(raster) library(rgdal) library(ggplot2) library(reshape2) library(plyr) # Assuming you have a path 'Maps' that you store your spatial files in. This # is all downloaded from <a href="http://www.naturalearthdata.com/downloads/">http://www.naturalearthdata.com/downloads/</a> using the # 1:50m "Medium" scale data. nat.earth <- stack('./Maps/NE2_50M_SR_W/NE2_50M_SR_W.tif') ne_lakes <- readOGR('./Maps/NaturalEarth/ne_50m_lakes.shp', 'ne_50m_lakes') ne_rivers <- readOGR('./Maps/NaturalEarth/ne_50m_rivers_lake_centerlines.shp', 'ne_50m_rivers_lake_centerlines') ne_coast <- readOGR('./Maps/NaturalEarth/ne_50m_coastline.shp', 'ne_50m_coastline') # 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]) x.subset } 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(nat.earth, y=extent(domain)) rast.table <- data.frame(xyFromCell(nat.crop, 1:ncell(nat.crop)), getValues(nat.crop/255)) rast.table$rgb <- with(rast.table, rgb(NE2_50M_SR_W.1, NE2_50M_SR_W.2, NE2_50M_SR_W.3, 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.

Reblogged this on Geospheres and commented:

Interesting post from downwithtime showing viable alternatives to ArcGIS map visualizations.

Pingback: “Lies, damn lies, and statistics” (and math) | ScienceBorealis.ca Blog