[R-sig-Geo] plot raster over google map

Barry Rowlingson b.rowlingson at lancaster.ac.uk
Thu Feb 28 18:52:06 CET 2013


On Thu, Feb 28, 2013 at 5:03 PM, Ross Ahmed <rossahmed at googlemail.com> wrote:
> I've downloaded the following google map:
>
> library(raster)
> myExtent <- matrix(c(-1.87710, -1.772096, 55.61399, 55.682171), ncol=2,
> byrow=T)
> g <- gmap(myExtent, type='hybrid')
> plot(g, interpolate=T)
>
> I'm trying to plot the raster over the google map:
>
> rownames(myExtent) <- c('long', 'lat')
> colnames(myExtent) <- c('min', 'max')
> rasterExtent <- extent(myExtent)
> myRaster   <- raster(ext=rasterExtent,crs="+proj=merc +a=6378137 +b=6378137
> +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null
> +no_defs")
> values(myRaster) <- rnorm(ncell(myRaster))

 You've created a raster with an extent in lat-long coordinates, but
assigned it a very non lat-long projection CRS.

 gmap has done a couple of things, it has got the google map for the
lat-long extent you requested, but returned a raster in google's CRS.
It's not a lat-long grid.

 > extent(g)
class       : Extent
xmin        : -213661.6
xmax        : -192526.8
ymin        : 7476426
ymax        : 7500886

 > extent(myExtent)
class       : Extent
xmin        : -1.8771
xmax        : -1.772096
ymin        : 55.61399
ymax        : 55.68217


> plot(myRaster, add=T, alpha=0.5)
>
> However, there's no sign of the raster in the plot. I presume theres a
> projection issue here. How can I get the raster to plot over the top of the
> google map?

 You need to stop lying to your raster, and tell it it really is a
raster in lat-long, then project it to the google map CRS:

 > projection(myRaster)="+init=epsg:4326"  # lat-long WGS84
 > pRaster =projectRaster(myRaster, crs=projection(g)) # reproject to
the google CRS (no need to type it all in again)
 > plot(pRaster,add=TRUE) # shazam

Barry



More information about the R-sig-Geo mailing list