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

Thu Apr 24 20:00:05 CEST 2014

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

> 
> 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've avoided ggmap et al because I've been unconvinced about their
> geographic rigour vs their prettiness....
> 
> 
> 
> 
> 
> On Thu, Apr 24, 2014 at 12:43 PM, Agustin Lobo <alobolistas at gmail.com
> <mailto:alobolistas at gmail.com>> wrote:
> 
>     But the fact that the output of get_map() includes the bounding box in
>     geographic
>     coordinates (epsg:4326) does not imply that the matrix is on the same
>     projection.
>     This should be specified.
>     I can provide you a utm raster along with bounding box coordinates in
>     geographic coordinates (epsg:4326).
>     Actually, according to the code I just showed, it is not clear at all
>     that the matrix provided
>     in the ggmap "raster" object is in epsg:4326.
> 
>     Agus
> 
> 
> 
>     On Thu, Apr 24, 2014 at 10:21 AM, Edzer Pebesma
>     <edzer.pebesma at uni-muenster.de
>     <mailto:edzer.pebesma at uni-muenster.de>> wrote:
>     > EPSG:4326, or WGS84, is only one way to refer to places on the
>     earth by
>     > longitudes and latitudes; it seems that google maps come in a Mercator
>     > projection, EPSG:3857, read e.g. this:
>     >
>     > http://docs.openlayers.org/library/spherical_mercator.html
>     >
>     > or
>     >
>     >
>     https://developers.google.com/maps/documentation/javascript/maptypes#MapCoordinates
>     >
>     > Iliffe and Lott: Datums and Map Projections: For Remote Sensing,
>     GIS and
>     > Surveying
>     >
>     > provide an excellent explanation on what goes on here, including what
>     > google maps does.
>     >
>     > On 04/24/2014 09:56 AM, Agustin Lobo wrote:
>     >> (follow up of previous message, I pressed enter as if I were in the R
>     >> session and sent the message...)
>     >>
>     >> On Wed, Apr 23, 2014 at 4:37 PM, Barry Rowlingson
>     >> <b.rowlingson at lancaster.ac.uk
>     <mailto:b.rowlingson at lancaster.ac.uk>> wrote:
>     >>>  - I imagine the web services expect EPSG:4326.
>     >>
>     >> What they expect, yes. But what about what we retrieve? The bbox
>     >> coordinates might be in geographic coordinates (epsg:4326) but
>     >> the actual projection of the matrix could be any other. Actually, I
>     >> think that no geographic object without CRS information should be
>     >> admitted in R.
>     >>
>     >> In order to investigate this further, I've done the following:
>     >>
>     >> mibb <- matrix(c(-69.88634, -48.78450, -34.05533,
>     -18.63424),byrow=TRUE,nrow=2)
>     >> delme <- get_map(location = mibb, maptype = "hybrid", source=
>     >> "google", crop = FALSE, zoom = 5)
>     >>
>     >> 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))
>     >> projection(rgmaprgb) <- CRS("+init=epsg:4326")
>     >> extent(rgmaprgb) <- unlist(attr(gmap,which="bb"))[c(2,4,1,3)]
>     >> rgmaprgb
>     >> plotRGB(rgmaprgb)
>     >>
>     >> This looks good, but when I save as GTiff
>     >>
>     writeRaster(rgmaprgb,file="rgmaprgb",format="GTiff",overwrite=TRUE,datatype="INT1U")
>     >>
>     >> and check with QGIS, there is a very significant shift versus the GM
>     >> layer that is downladed by plugin OpenLayers.
>     >>
>     >> Thus I worry that my overlay using rasterVis:levelplot and ggmap:
>     >> layer could actually be wrong as well.
>     >>
>     >> I'm going to test on a more local area and let you know.
>     >>
>     >> Agus
>     >>
>     >> _______________________________________________
>     >> R-sig-Geo mailing list
>     >> R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>
>     >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>     >>
>     >
>     > --
>     > Edzer Pebesma
>     > Institute for Geoinformatics (ifgi), University of Münster
>     > Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251
>     > 83 33081 http://ifgi.uni-muenster.de GPG key ID 0xAC227795
>     >
>     >
>     > _______________________________________________
>     > R-sig-Geo mailing list
>     > R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>
>     > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>     >
> 
>     _______________________________________________
>     R-sig-Geo mailing list
>     R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>
>     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 
> 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251
83 33081 http://ifgi.uni-muenster.de GPG key ID 0xAC227795

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20140424/abae40c4/attachment.bin>