Home > Uncategorized > More Graphics with Google earth

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:

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.

About these ads
Categories: Uncategorized
  1. November 11, 2010 at 5:10 PM | #1

    Can we clean up the spam please? (sign of popularity, by the way!)

  2. Steven Mosher
    November 11, 2010 at 10:17 PM | #2

    thanks. done

  3. Gláuber
    June 26, 2012 at 9:26 AM | #3

    How can I add lat and lon axes in the image like the first image?

    • Steven Mosher
      June 26, 2012 at 9:35 AM | #4

      Arrg. That was something I never figured out. I think if you add

      ylab = “Latitude”, xlab = “Longitude” that it might work. either to the plot
      command or the contour command

  4. July 26, 2013 at 1:10 AM | #5

    Hello Mr Steven
    When I was running your lines in R.
    First Send me error in -world_avg.tif-

    do you have world avg. tif file.

    and second send me Error in crop

    thanks a lot

    Hernando Hernandez.

    look at these

    > r e rc g mrc plot(g, interpolate=TRUE)
    > contour(mrc, add=TRUE, lwd=2, levels=1:7 * 20, col=rev(heat.colors(7)))
    Error en contour(mrc, add = TRUE, lwd = 2, levels = 1:7 * 20, col = rev(heat.colors(7))) :
    error in evaluating the argument ‘x’ in selecting a method for function ‘contour’: Error: objeto ‘mrc’ no encontrado
    > pt mpt points(mpt, col=’light blue’, pch=’+’, cex=2)

  1. No trackbacks yet.

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


Get every new post delivered to your Inbox.

Join 31 other followers

%d bloggers like this: