[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