The CRS of the object returned by get_map() and its consequences for plotting overlaid maps in R

Agustin.Lobo@ictja.csic.es
20140506
Just R code

The CRS of objects returned by get_map() is not documented. Here we use an empirical approach to confirm that, as it is usual in tools interacting with GM, the map is in th so called “Google PseudoMercator” while coordinates are provided in “Longitude, Latitude” (datum WGS84 in both cases), and show how to produce correct plots of overlaid maps

NOTE 20140508: This is kind of obsolete information as it is more simple using dismo:gmap() instead of ggmap:get_map() to download the GM layer. This is because dismo:gmap(), unlike ggmap:get_map(), returns a raster layer from package raster including complete CRS information

Requirements

require(utils)
require(rgdal)

## Loading required package: rgdal
## Loading required package: sp
## rgdal: version: 0.8-16, (SVN revision 498)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.10.0, released 2013/04/24
## Path to GDAL shared files: /usr/share/gdal/1.10
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)

require(raster)

## Loading required package: raster

require(ggmap)

## Loading required package: ggmap
## Loading required package: ggplot2

Test object

gmap <- get_map(location = c(5, 51), maptype = "hybrid", source = "google", 
    crop = FALSE, zoom = 6)

1. The CRS of get_map() output

1.1 Convert ggmap output to raster brick layer (package raster)

mgmap <- as.matrix(gmap)
vgmap <- as.vector(mgmap)
vgmaprgb <- col2rgb(vgmap)
gmapr <- matrix(vgmaprgb[1, ], ncol = ncol(mgmap), nrow = nrow(mgmap))
gmapg <- matrix(vgmaprgb[2, ], ncol = ncol(mgmap), nrow = nrow(mgmap))
gmapb <- matrix(vgmaprgb[3, ], ncol = ncol(mgmap), nrow = nrow(mgmap))
rgmaprgb <- brick(raster(gmapr), raster(gmapg), raster(gmapb))
rm(gmapr, gmapg, gmapb)

1.2 Assume output of get_map() is in “Lon,Lat WGS84” (epsg:4326)

projection(rgmaprgb) <- CRS("+init=epsg:4326")
extent(rgmaprgb) <- unlist(attr(gmap, which = "bb"))[c(2, 4, 1, 3)]
rgmaprgb

## class       : RasterBrick 
## dimensions  : 1280, 1280, 1638400, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 0.01099, 0.006919  (x, y)
## extent      : -2.02, 12.04, 46.35, 55.21  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer.1, layer.2, layer.3 
## min values  :       0,       0,       0 
## max values  :     254,     254,     254

plotRGB(rgmaprgb)

Save as GTiff:

writeRaster(rgmaprgb, file = "rgmaprgb", format = "GTiff", overwrite = TRUE, 
    datatype = "INT1U")

## class       : RasterBrick 
## dimensions  : 1280, 1280, 1638400, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 0.01099, 0.006919  (x, y)
## extent      : -2.02, 12.04, 46.35, 55.21  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /media/alobo/Iomega_HDD/RrasterGuide/rasterVisTutorial/rgmaprgb.tif 
## names       : rgmaprgb.1, rgmaprgb.2, rgmaprgb.3 
## min values  :          0,          0,          0 
## max values  :        254,        254,        254

1.3 Assume output of get_map() is in “Google PseudoMercator” (epsg:3857)

rgmaprgbGM <- rgmaprgb
projection(rgmaprgbGM) <- CRS("+init=epsg:3857")

In order to set the extent (required to define a correct raster layer), we have to find out the extent in the units of the new projection
We create an SpPolDF and then use spTransform

unlist(attr(gmap, which = "bb"))[c(2, 4, 1, 3)]

## ll.lon ur.lon ll.lat ur.lat 
##  -2.02  12.04  46.35  55.21

rprobextSpDF <- as(extent(unlist(attr(gmap, which = "bb"))[c(2, 4, 1, 3)]), 
    "SpatialPolygons")
projection(rprobextSpDF) <- CRS("+init=epsg:4326")
rprobextGM <- spTransform(rprobextSpDF, CRS("+init=epsg:3857"))
rprobextGM@bbox

##       min     max
## x -224895 1340536
## y 5837356 7402786

extent(rgmaprgbGM) <- c(rprobextGM@bbox[1, ], rprobextGM@bbox[2, ])
rgmaprgbGM

## class       : RasterBrick 
## dimensions  : 1280, 1280, 1638400, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 1223, 1223  (x, y)
## extent      : -224895, 1340536, 5837356, 7402786  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:3857 +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 
## data source : in memory
## names       : layer.1, layer.2, layer.3 
## min values  :       0,       0,       0 
## max values  :     254,     254,     254

plotRGB(rgmaprgbGM)

Save as GTiff:

writeRaster(rgmaprgbGM, file = "rgmaprgbGM", format = "GTiff", overwrite = TRUE, 
    datatype = "INT1U")

## class       : RasterBrick 
## dimensions  : 1280, 1280, 1638400, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 1223, 1223  (x, y)
## extent      : -224895, 1340536, 5837356, 7402786  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs 
## data source : /media/alobo/Iomega_HDD/RrasterGuide/rasterVisTutorial/rgmaprgbGM.tif 
## names       : rgmaprgbGM.1, rgmaprgbGM.2, rgmaprgbGM.3 
## min values  :            0,            0,            0 
## max values  :          254,          254,          254

1.4 Overlay of a reference vector layer on the GTiff images using QGIS

In order to avoid any confusion introduced by overlaying within R, we overlay a reference vector layer on the raster layers exported as GTiff in QGIS
“Google PseudoMercator” (epsg:3857)

The overlay of a reference vector layer on the raster layer from get_map() assuming “Google PseudoMercator” is correct
“Lon,Lat WGS84” (epsg:4326)
The overlay of a reference vector layer on the raster layer from get_map() assuming “Lon,Lat WGS84” is shifted

According to these findings, the object returned by get_map()

str(gmap)

##  chr [1:1280, 1:1280] "#0B130B" "#1B221B" "#C5C9C5" "#F2F2F2" ...
##  - attr(*, "class")= chr [1:2] "ggmap" "raster"
##  - attr(*, "bb")='data.frame':   1 obs. of  4 variables:
##   ..$ ll.lat: num 46.4
##   ..$ ll.lon: num -2.02
##   ..$ ur.lat: num 55.2
##   ..$ ur.lon: num 12

consists of:

• a matrix which grid is projected in *“Google PseudoMercator” (epsg:3857)*
• the coordinates of the bounding box provided in *epsg:4326*

2. Consequences for the Display of overlaid maps within R

Care must be taken with the different CRS of the grid and the coordinates for overlaying layers in a plot.
In order to test the overlayig within R, we import the reference vector layer into R and create a reprojected version:

download.file("https://dl.dropboxusercontent.com/u/3180464/NWEur.zip", "NWEur.zip", 
    method = "curl")
unzip("NWEur.zip", exdir = "NWEur")
NWEur <- readOGR(dsn = "./NWEur", layer = "NWEur", stringsAsFactors = FALSE)

## OGR data source with driver: ESRI Shapefile 
## Source: "./NWEur", layer: "NWEur"
## with 9 features and 18 fields
## Feature type: wkbPolygon with 2 dimensions

NWEurGM <- spTransform(NWEur, CRS("+init=epsg:3857"))

2.1 Standard plots

Overlay using standard plots is correct if we assume the output of get_map() is in “Google PseudoMercator” projection and use the vector layer with the same CRS

plotRGB(rgmaprgbGM)
plot(NWEurGM, add = TRUE, lwd = 2, ext = extent(rgmaprgbGM))

# 

2.2. Using levelplot() to plot the output of get_map() and layer() to overlay the reference polygons layer

(or eventually another raster)

require(latticeExtra)

## Loading required package: latticeExtra
## Loading required package: RColorBrewer
## Loading required package: lattice
## 
## Attaching package: 'latticeExtra'
## 
## The following object is masked from 'package:ggplot2':
## 
##     layer

require(rasterVis)

## Loading required package: rasterVis
## Loading required package: hexbin
## Loading required package: grid

# NOTE: the order of loading these packages matters!

levelplot() requires a raster layer, thus the output of get_map() cannot be used directly: We use the raster version that we calculated above.
The plot is shifted if we assume get_map() output is in “Longitude, Latitude”:

levelplot(rgmaprgb[[3]], margin = FALSE, par.settings = GrTheme) + layer(sp.polygons(NWEur, 
    col = "yellow"))

and correct if use layers in “Google PseudoMercator”:

levelplot(rgmaprgbGM[[3]], margin = FALSE, par.settings = GrTheme) + layer(sp.polygons(NWEurGM, 
    col = "yellow"))

The drawback of this solution is that we are missing the color in the raster layer, although making a pseudocolor version of the get_map() output is possible. TBD

2.3 Using levelplot() to plot the output of get_map() and spplot() to overlay the reference polygons layer

We can also use spplot() instead of layer() for plotting the polygons:

levelplot(rgmaprgbGM[[3]], margin = FALSE, par.settings = GrTheme) + spplot(NWEurGM, 
    zcol = "NAME", fill = "transparent", col = "yellow", xlim = c(extent(rgmaprgbGM)@xmin, 
        extent(rgmaprgbGM)@xmax), ylim = c(extent(rgmaprgbGM)@ymin, extent(rgmaprgbGM)@ymax), 
    colorkey = FALSE)

2.4 Using spplot() to plot the reference polygons layer and layer() to plot the output of get_map() under it

layer() does not work with raster layers, therefore cannot use the rgmaprgbGM object.
layer() does not work with the output of get_map() directly either, but can use grid.raster().
grid.raster() requires some explicit geographic information:

bbMap <- attr(gmap, "bb")
latCenter <- with(bbMap, ll.lat + ur.lat)/2
lonCenter <- with(bbMap, ll.lon + ur.lon)/2
height <- with(bbMap, ur.lat - ll.lat)
width <- with(bbMap, ur.lon - ll.lon)

As above, both layers are shifted if we assume raster and vector layers are in “Longitude, Latitude”:

spplot(NWEur, zcol = "NAME", xlim = unlist(bbMap[c(2, 4)]), ylim = unlist(bbMap[c(1, 
    3)]), fill = "transparent", col = "yellow", colorkey = FALSE) + layer(grid.raster(gmap, 
    x = lonCenter, y = latCenter, width = width, height = height, default.units = "native"), 
    under = TRUE)

The vector layer must be in “Google PseudoMercator” to match the output of get_map(), but note that using the geographic attributes directly provides a wrong georeferencing and, as a consequence, no overlay:

spplot(NWEurGM, zcol = "NAME", fill = "transparent", col = "yellow", colorkey = FALSE) + 
    layer(grid.raster(gmap, x = lonCenter, y = latCenter, width = width, height = height, 
        default.units = "native"), under = TRUE)

The geographic attributes must be reprojected to “Google PseudoMercator” (epsg:3857) to be consistent with the grid: Here we take advantage of the extent of the raster layer calculated in section 1.3

xcenter <- (extent(rgmaprgbGM)@xmax + extent(rgmaprgbGM)@xmin)/2
ycenter <- (extent(rgmaprgbGM)@ymax + extent(rgmaprgbGM)@ymin)/2
height <- extent(rgmaprgbGM)@xmax - extent(rgmaprgbGM)@xmin
width <- extent(rgmaprgbGM)@ymax - extent(rgmaprgbGM)@ymin

The correct display is thus:

spplot(NWEurGM, zcol = "NAME", fill = "transparent", col = "yellow", xlim = c(extent(rgmaprgbGM)@xmin, 
    extent(rgmaprgbGM)@xmax), ylim = c(extent(rgmaprgbGM)@ymin, extent(rgmaprgbGM)@ymax), 
    colorkey = FALSE) + layer(grid.raster(gmap, x = xcenter, y = ycenter, width = width, 
    height = height, default.units = "native"), under = TRUE)

3. Conclusions

3.1 The object returned by ggmap:get_map() consists of

• a matrix which grid is projected in “Google PseudoMercator” (epsg:3857)
• the coordinates of the bounding box provided in epsg:4326 ### 3.2 The object returned by dismo:gmap() is a raster layer (in pseudocolor) (from package raster) including CRS information (“Google PseudoMercator” (epsg:3857)) ### 3.3 GTiff output must include the CRS and extent details in epsg:3857 ### 3.4 Plotting overlaid maps must be done with both layers on the same CRS:
• In standard R plots, displaying the Raster layer and then overlaying the vector
  • plot( r ); plot(v, add=TRUE)
• Displaying the Raster layer first and then overlaying the vector (or another raster) layer on top:
  • levelplot() + layer()
  • levelplot() + spplot()
• Displaying the Vector layer first and then puting the raster layer under it:
  • spplot() + layer()

a