More Graphics with Google earth

Dr. Paul, R graphics guru, blessed us with his rendition of transparent contour maps drawn on a google earth image: Cool stuff. I’ll be taking his code and turning it into a function and sharing it back: here is the picture his code creates:

That is just plain slick.  While I was looking over his code a couple more examples hit my inbox. One using ggplot which I havent got to ( see the comments on the previous post) and the other example came from the raster guru, Robert. First, lets take a look at his code:

library(dismo)
r <- raster('world_avg.tif')
e <- extent(-122.6, -122.3, 37.75, 37.9)
rc <- crop(r, e)
g <- gmap(e, type='satellite')
mrc <- projectRaster(rc, g)
plot(g, interpolate=TRUE)
contour(mrc, add=TRUE, lwd=2, levels=1:7 * 20, col=rev(heat.colors(7)))
pt <- cbind(-122.45, 37.8)
mpt <- Mercator(pt)
points(mpt, col='light blue', pch='+', cex=2)

First we load the package dismo. Now ordinarily you would not go looking for resources to do graphics in a package that dealt with species distribution. In addition to working on raster Robert also works on dismo. It has a couple of tools we need. The first is g <- gmap(e, type=’satellite’).  The call to gmap replaces my calls to GetMap.bbox() from RgoogleMaps, in fact it is based on that code. I happen to like the interface of gmap better. Also, there is no annoying *rda file written for every map. And the pesky GDAL warnings are gone ( you could just suppress the warnings dummy) Anyway, gmap issues a call to the static map server and the pixels are returned to the object “g”. next comes a projection: mrc <- projectRaster(rc, g). This function projects the “rc” raster to the mrc raster using the coordinate system of “g” ( err, I think) basically, the g raster has a mercator projection, and rc is unprojected. so we need to get a mercator projection of rc. Then, we just plot g. On top of that we add a contour of mrc ( which is just rc projected in a mercator projection). Then we set the levels and select heat.colors. I have to study this more to understand how to select the number of levels. We’ll see how this selection plays out. ( I prolly need to vary it based on the range of the data)  Ideally, I’d like the levels to correspond to the Imhoff urbanity levels. Then we draw the cross hair. We pick a point ( where the station is at) then then we take that point and transform that point using a mercator projection. Then we draw the point and we are done.  Below find the three charts we did as samples in the previous post. The key to the first two was seeing that the nightlights data was high frequency so that position accuracy of the station matters. The third chart shows how minor errors in positioning of the station can lead to mis identification. The station location given by the inventory is not coincident with the location of the place name. If we look at the place name, where the station is likely to be positioned, we see different nightlights. So positional innacuries, even small ones, can change how we categorize sites Rural/periurban/urban. hence, it is important for  folks to either improve the accuracy of the data or employ exclusionary rules to preclude errors that may arise from mislocations.