Skip to content

projection of ggmap:get_map() output

11 messages · Barry Rowlingson, Jon Olav Skoien, Edzer Pebesma +2 more

Original
Edzer Pebesma
#
On 04/24/2014 03:52 PM, Barry Rowlingson wrote:
Indeed:

library(ggmap)
ggmap(get_map(matrix(c(-10,20,40,80),2,2)))

and look at the tics along the y-axis: they're not linear, neither are
they drawn by google (at least, they're not on the ggmapTemp.png file
written).

http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf

mention on p. 146 that "... the coordinate system is ?xed to the
Mercator projection."
The google API wants long/lat in WGS84, and what gets returned has a
different bbox anyway, but I believe that a long/lat aligned "box" in
WGS84 remains a box in Mercator.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Oscar Perpiñan
    

    

      
      #
    

  


  

          

      
      
      
    

        
Just as additional confirmation: inspecting the code, ggmap uses the
mercator projection for displaying the "raster":
https://github.com/dkahle/ggmap/blob/master/R/ggmap.R#L562

            

      
      
      
    

        
I have been playing with code to show the differences. If my code is
correct, the differences are very small, at least for the bounding box
I have chosen:

library(geosphere) ## provides a mercator function

## Bounding box: lower-left and upper-right
ll <- c(-10, 30)
ur <- c(10, 50)
## Bounding box in mercator projection
llxy <- mercator(ll)
urxy <- mercator(ur)

## Matrix of coordinates in mercator projection
xy <- expand.grid(x=seq(llxy[1], urxy[1], length=100),
                  y=seq(llxy[2], urxy[2], length=100))

## Matrix of lon-lat values using bounding box
lonlat <- expand.grid(x=seq(ll[1], ur[1], length=100),
                  y=seq(ll[2], ur[2], length=100))

## Mercator coordinates using matrix of lon-lat
xyproj <- mercator(lonlat)

## Differences between xy and xyproj
dif <- 1 - xyproj/xy
names(dif) <- paste0('dif', names(dif))
summary(dif)

library(lattice)
difXY <- cbind(xy, dif)
levelplot(difx + dify ~ x*y, data=difXY)

            

      
      
      
    

        
IMHO ggmap and friends are useful only when you need an easy way to
add an attractive context or "decorative" material to a graphic.

Best,

Oscar.

  
  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
According to the results
https://dl.dropboxusercontent.com/u/3180464/ggmap_GM_compare.odp
of the code below, the matrix downloaded by get_map() is on epsg:3857,
the so-called
Google Pseudo Mercator, while the bounding box is given in epsg:4326
coordinates.
The odp file shows the raster as created by the code below out of a
get_map() command,
and the same area displayed by plugin OpenLayers.
Note there is still a bit of shift/distortion between the 2 when you
compare alternatively,
but this is much less than what I had assuming get_map() was
downloading in epsg:4326

mibb <- matrix(c(-69.88634, -48.78450, -34.05533, -18.63424),byrow=TRUE,nrow=2)
gmap <- get_map(location = mibb, maptype = "hybrid", source= "google",
crop = FALSE, zoom = 5)

# Make a raster brick object from gmap
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))

#version in "Google PseudoMercator" epsg:3857
rgmaprgbGM <- rgmaprgb
projection(rgmaprgbGM) <- CRS("+init=epsg:3857")
#We have to find out the extent in the units of the new projection
unlist(attr(gmap,which="bb"))[c(2,4,1,3)]
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 at bbox
extent(rgmaprgbGM) <- c(rprobextGM at bbox[1,],rprobextGM at bbox[2,])
rgmaprgbGM
plotRGB(rgmaprgbGM)
writeRaster(rgmaprgbGM,file="rgmaprgbGM",format="GTiff",overwrite=TRUE,datatype="INT1U")

I therefore think that overlaying layers using levelplot() + layer()
might have to be done with the user's layers
in epsg:4326, but I'll confirm soon. It could be as well that ggmap
reprojects from 4326 to 3857.
The query for get_map() must be done in epsg:4326, though, that's for sure.

ggmap should specify the projection of the object retrieved by
get_map(). The mentions you find to "mercator" are not unambiguous, as
"mercator" is not one projection but a projection family.

Geo-graphics must be pretty **and** exact.

Agus

        
On Fri, Apr 25, 2014 at 12:13 PM, Oscar Perpi?an <oscar.perpinan at upm.es> wrote:

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
I've tested the effect of the projections on overlaying graphs using
levelplot() + layer()
In short, the user's object must be in epsg:4326, probably because
layer() carries out a reprojection.
The bad news is that it carries out a wrong reprojection, perhaps
because the incorrect
(or insufficient) definition of the projection as just "mercator" in layer()

Example:
(shape file in
https://dl.dropboxusercontent.com/u/3180464/NWEur.zip)

gmap2 <- get_map(location = c(5,51), maptype = "hybrid", source=
"google", crop = FALSE, zoom = 6)
ggmap(gmap2)
#convert gmap2 to rgb raster brick
mgmap2 <- as.matrix(gmap2)
vgmap2 <- as.vector(mgmap2)
vgmap2rgb <- col2rgb(vgmap2)
gmap2r <- matrix(vgmap2rgb[1,],ncol=ncol(mgmap2),nrow=nrow(mgmap2))
gmap2g <- matrix(vgmap2rgb[2,],ncol=ncol(mgmap2),nrow=nrow(mgmap2))
gmap2b <- matrix(vgmap2rgb[3,],ncol=ncol(mgmap2),nrow=nrow(mgmap2))
rgmap2 <- brick(raster(gmap2r),raster(gmap2g),raster(gmap2b))
extent(rgmap2) <- unlist(attr(gmap2,which="bb"))[c(2,4,1,3)]
projection(rgmap2) <- CRS("+init=epsg:3857")
plotRGB(rgmap2)
plot(rgmap2[[3]],col=grey.colors(32,start=0,end=1,gamma=1))
#shapefiles to be overlaid: coastlines easy to follow in the resulting graphic
NWEur <- readOGR(dsn="countries_shp",layer="NWEur",stringsAsFactors=FALSE)
projection(NWEur)
NWEurGM <- spTransform(NWEur,CRS("+init=epsg:3857"))
#levelplots
levelplot(rgmap2[[3]],margin=FALSE,par.settings=GrTheme) +
  layer(sp.polygons(NWEur))
levelplot(rgmap2[[3]],margin=FALSE,par.settings=GrTheme) +
  layer(sp.polygons(NWEurGM))

overlay with the layer in epsg3857 is totaly out of place,
but the one with the layer in epsg4326 is clearly shifted.

Thus, we cannot use levelplot() + layer() with get_map() output at the moment.

Agus

        
On Fri, Apr 25, 2014 at 3:55 PM, Agustin Lobo <alobolistas at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    
2 days later

    

      
      

        


  

        

      Oscar Perpiñan
    

    

      
      #
    

  


  

      
Hello,

I couldn't run your code because I am not able to read the file.
Anyway, some comments about your last email:

- "layer" is a function of the latticeExtra package. This function
modifies the panel function of a "trellis" graphic produced with a
function of the lattice package. This function, together with
"+.trellis" (of the same package) allows for overlaying layers. It is
important to note that it does not carry out any projection. It only
displays the graphical result of another function (for example,
sp.polygons or whatever).

- You create the canvas with "levelplot" displaying the ggmap raster,
and then overlay the data with layer. This is the same approach that
ggmap uses. Instead, I prefer to rely on the main data to be displayed
to configure the graphical window instead of using an auxiliary image
as the starting point.

- I'm not a fan of the ggmap approach. I agree with you that the
resulting object should include the projection information and even
adhere to the sp classes. However, I think that their results are
useful if managed with care.

In order to show that this is not an issue strictly related to
rasterVis::levelplot, I have written some code with spplot. I think
that the concordance between the ggmap raster and the shapefile is
aceptable.

Best,

Oscar.

##########################################################
library(sp)
library(ggmap)
## latticeExtra must be loaded after ggmap because both ggplot2 and
## latticeExtra define a 'layer' function. We need the definition from
## latticeExtra.
library(latticeExtra)

## Shape example
library(maptools)
owd <- getwd()
setwd(system.file("shapes", package = "maptools"))
sc90 <- readShapeSpatial("co45_d90")
proj4string(sc90) <- CRS("+proj=longlat +datum=NAD27")
setwd(owd)

## Get the stamen map
bbPoly <- bbox(sc90)
gmap <- get_map(c(bbPoly), maptype='toner', source='stamen', crop=FALSE)
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)

## Use latticeExtra::+.trellis and latticeExtra::layer
spplot(sc90['AREA'], fill='transparent', col='indianred1') +
                  layer(grid.raster(gmap,
                                    x=lonCenter, y=latCenter,
                                    width=width, height=height,
                                    default.units='native'),
                        under=TRUE)

  
  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
Nice, thank you!

Markus Loecher, author of RgoogleMaps, explained that ggmap uses

https://github.com/dkahle/ggmap/blob/master/R/LonLat2XY.R

to go from long/lat to XY (pixels in google map tiles), which was taken
from the RgoogleMaps one:

https://github.com/markusloecher/RgoogleMaps/blob/master/R/LatLon2XY.R

which points to code by someone called Neil Young in a google groups:

https://groups.google.com/forum/?hl=en#!topic/Google-Maps-API/0hA6wp6VaW8

I would expect that your cylindrical, equidistant projection (linear in
long/lat) must collide more with the Mercator used by google when the
plotted area becomes larger.

        
On 04/28/2014 11:58 AM, Oscar Perpi?an wrote:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
Oscar,

        
On Mon, Apr 28, 2014 at 11:58 AM, Oscar Perpi?an <oscar.perpinan at upm.es> wrote:

            

      
      
      
    

        
Note dsn is the directory where you have the 4 files and layer is the
name of the layer, which is the name of the files without extension.

            

      
      
      
    

        
There has to be a reprojection on the fly because the the two layers of
the plot are on different coordinate systems. I think it is clear now that
ggplot uses the Google "PseudoMercator" for the matrix while expects
all user input on epsg4326

            

      
      
      
    

        
Users cannot handle with care undocumented idiosyncratic features.

            

      
      
      
    

        
I agree this is not an specific problem of reasterVis:levleplot.

            

      
      
      
    

        
I think "acceptable" depends on the scale only. At coarser scale you
just do not see the error
because it is relative to the extent and the result converges to the
correct  reprojection.

I do not think we should settle for this kind of errors.

In the code you provide, how do you know NAD27 is the appropriate datum?
There is no prj file (one of the big shortcomings of shapefiles is
that the presence
of crs information in the form of a prj file is not compulsory). Note
that you define
a datum != WGS84 and you do not define the towgs argument.

Agus

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
This is great information, thanks.
It is clear that without an infrastructure such as the proj4 used by
spTransform() the overlays will never be
without error.

Agus

On Mon, Apr 28, 2014 at 4:57 PM, Edzer Pebesma

        
<edzer.pebesma at uni-muenster.de> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Oscar Perpiñan
    

    

      
      #
    

  


  

          

      
      
      
    

        
My fault, sorry.

            

      
      
      
    

        
I am not sure what do you mean with "reprojection on the fly" but
"layer" does not make any reprojection.
In fact, it is only a wrapper around
"eval(substitute(expression(...)))" to modify easily the panel
function of a trellis graphic.
However, in my code I am using "grid.raster" with its arguments width,
height, x and y, to resize the ggmap raster and to position it in the
graphic window created with spplot.

            

      
      
      
    

        
I think it depends on the extent. I have tested the code again with
the administrative boundaries of Spain and some of its regions
(http://biogeo.ucdavis.edu/data/gadm2/shp/ESP_adm.zip).
It fails for the country extent (as Edzer suggested and you have found
with your file) but it works correctly for a smaller extent (code
copied below).

            

      
      
      
    

        
It is an example extracted from the p.121 of the 2008 edition of
"Applied Spatial Data Analysis with R". That chapter is now available
here: http://cran.r-project.org/web/packages/maptools/vignettes/combine_maptools.pdf

Best,

Oscar.

#######################################################
library(ggmap)
library(latticeExtra)
library(sp)
library(maptools)

## http://biogeo.ucdavis.edu/data/gadm2/shp/ESP_adm.zip
setwd('/tmp')
spain <- readShapePoly('ESP_adm/ESP_adm1')
proj4string(spain) <- CRS("+proj=longlat +datum=WGS84")

andalucia <- spain[spain$ID_1 == 1,]
catalonia <- spain[spain$ID_1 == 6,]

## myPoly <- spain
myPoly <- andalucia
## myPoly <- catalonia

bbPoly <- bbox(myPoly)
gmap <- get_map(c(bbPoly), maptype='watercolor', source='stamen', crop=FALSE)
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)

spplot(myPoly["ID_1"], fill='transparent') +
                  layer(grid.raster(gmap,
                                    x=lonCenter, y=latCenter,
                                    width=width, height=height,
                                    default.units='native'),
                        under=TRUE)