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

Fri Apr 25 20:15:42 CEST 2014

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:
> 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 ﬁ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