projection of ggmap:get_map() output

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 04/24/2014 03:52 PM, Barry Rowlingson wrote:

So do you mean (and I've not tried your code yet, sorry) that the
corners may well be at the given lat-long points but half way along any
edge might not be at the halfway-point in lat-long?

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."

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

if you convert the bbox coords to epsg:3857 then maybe you have a
correctly projected raster in those coordinates?

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.

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)

I've avoided ggmap et al because I've been unconvinced about their
geographic rigour vs their prettiness....

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.

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo