Scott Rinnan

Quantitative Ecology and Resource Management, University of Washington

Troubleshooting the gmap function for making maps in R

2 Comments

There are several wonderful tools and packages available for incorporating map features from into plots in R. Some of the more well-known examples are the ggmap function from the aptly-named ggmap package (which follows the conventions of Hadley Wickham’s Grammar of Graphics)1, and the GetMap function from the RgoogleMaps package. These both can query and pull data from Google Maps:

library(ggmap)
rainier<-get_map(location=c(-121.760374,46.852886),maptype="satellite",zoom=11)
ggmap(rainier)

unnamed-chunk-1-1

You can specify the location by coordinates, or you can type in a location just like you would in a Google Maps searchbar. You can even liven things up a little with some more playful styles:

seattle<-get_map(location="Seattle, WA",maptype="watercolor")
ggmap(seattle)

unnamed-chunk-2-1

After several months of working with the dismo package for other research purposes, I found a new (to me) function in the package that I really like. The gmap function in the dismo package is unique in that it returns the map query to you in raster format, which can be extremely useful for plotting purposes when you’re working with spatial data.

library(dismo)
wash<-gmap("Washington state",scale=2)
plot(wash)

unnamed-chunk-3-1unnamed-chunk-3-2

Since this map is in raster format, all the familiar raster and other spatial analysis functions can be used, such as crop, mask, etc. You can also make a map by giving the gmap function an extent object. Seriously, this is sooooo useful!

library(maptools)
sagegrouse<-readShapePoly("~/Google Drive/JDSM/Cent_urop_pl.shp")
grouse_habitat<-gmap(extent(sagegrouse),lonlat=T,type="satellite")
plot(grouse_habitat)
plot(sagegrouse,add=T,col="#00FF0040")

unnamed-chunk-4-1unnamed-chunk-4-2

Here I’ve read in a polygon shapefile that represents the historical distribution of greater sage-grouse, and used the extent of the shapefile to get a background map. And since the map is just a raster object, we can use the regular plot function for the map, and it will automatically use the plot method for plotting rasters. Pretty slick, huh?

But wait, what’s up with all that white space that precedes all those plots?

Well, turns out there’s something wonky about the format of the raster that gets returned from the Google database. I don’t know what it is, precisely, but as near as I can tell, when you do the command plot(map) where map is a raster object that’s been created with the gmap function, this gets plotted with the rasterImage method from the graphics package, rather than the the plot method for rasters in the raster package. Why does this matter? Because the rasterImage function has all sorts of stability issues, and is not even supported on some operating systems. I don’t know why the maps use that method rather the plot.raster method, but that’s what they’re doing, and that big white space is an artifact of that. This is a particularly frustrating issue if you are trying to write a figure to a file using a function like pdf; if you do this, it will create a file with two pages, the first one blank, the second one your plot.

Even worse, the rasterImage plotting method isn’t compatible with really useful tools, like a legend. What happens if we try to add a legend to one of these plots?

sagegrouse_breadth<-raster("~/Google Drive/Breadth_research/Historical_breadth_maps/Birds/Hist_cent_urop.tif")
plot(grouse_habitat)
cols<-colorRampPalette(c("black","darkgreen","olivedrab3","yellow","red3"))
plot(sagegrouse_breadth,add=T,legend=F,col=cols(200))
plot(sagegrouse_breadth,legend.only=T,col=cols(200))

unnamed-chunk-5-1

Yeah, that’s not helpful.

In certain cases, you can force which type of method you want to use by making a call to that method directly. I illustrate that here by creating some fake time series data, and plotting that with both the generic plot method and the plot.ts method for time series.

x<-1990:2010
y<-sample(20,21,replace=T)
plot(x,y)
plot.ts(x,y)

unnamed-chunk-6-1unnamed-chunk-6-2Unfortunately, we can’t do the same with the method for plotting rasters from the raster package. Or maybe we can, I don’t know; this is running up against the limits of my R knowledge.

Instead, I wrote a function that takes the gmap object and forces it to get plotted with the raster method, rather than the rasterImage method. It’s not an ideal fix, but it gets the job done.

plot_gmap<-function(ras,...){
cols<-ras@legend@colortable
x<-unique(ras)
par(mar=c(0,0,0,0))
plot(ras,col=cols[x+1],legend=F,box=F,axes=F,...)
}

Using this function to plot the gmap object, the results should now be in a format that more closely resembles what we want:

plot_gmap(grouse_habitat)
cols<-colorRampPalette(c("black","darkgreen","olivedrab3","yellow","red3"))
plot(sagegrouse_breadth,add=T,legend=F,col=cols(200))
plot(sagegrouse_breadth,legend.only=T,col=cols(200))

unnamed-chunk-8-1

There we go, that’s more like it! The plot_gmap function can accept any argument that you would use when plotting a raster object with the raster package, such as alpha (controls transparency), ext (controls extent), as well as any other graphical parameters such as main (makes a title).

There are still some plotting issues, to be sure. If you try and crop a gmap object, you’ll notice that the color scheme gets messed up:

seattle<-gmap("Seattle, WA",type="satellite",scale=2)
plot_gmap(seattle)
seattle_crop<-crop(seattle,extent(seattle)/2)
plot_gmap(seattle_crop)

unnamed-chunk-9-1unnamed-chunk-9-2

As near as I can tell, this is not a fault of the function I wrote, and more a consequence of cropping rasters that use a colortable to specify colors. Colortables are a relatively new feature in the raster package, and several of the package’s functions have not yet been rewritten to preserve the structure of colortables when objects gets manipulated.

Long story short: gmap is a great function, but there are frustrating plotting issues that come along with it. Hopefully the plot_gmap function I wrote above helps you get around those issues. Good luck, and happy plotting!

References:

  1. Wickham, H. (2010). A layered grammar of graphics. Journal of Computational and Graphical Statistics, 19(1), 3-28.
Advertisements

2 thoughts on “Troubleshooting the gmap function for making maps in R

  1. Hi Scott,

    Thanks for the post, very useful! It inspired me to look into the plotting code of raster package, and it turns out that to avoid the blank first plot we just need to comment a line calling (incorrectly) plot.new.

    In case future readers may find it useful, I have uploaded the modified function here: https://github.com/Pakillo/rSDM/blob/master/R/plot_gmap.R.

    Cheers!

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