[R-sig-Geo] projection of ggmap:get_map() output

[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:
>> 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 fixed 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