The CRS of the object returned by get_map() and its consequences for plotting overlaid maps in R

The CRS of objects returned by get_map() is not documented. Here we use an empirical approach to confirm that, as it is usual in tools interacting with GM, the map is in th so called “Google PseudoMercator” while coordinates are provided in “Longitude, Latitude” (datum WGS84 in both cases), and show how to produce correct plots of overlaid maps NOTE 20140508: This is kind of obsolete information as it is more simple using dismo:gmap() instead of ggmap:get_map() to download the GM layer. This is because dismo:gmap(), unlike ggmap:get_map(), returns a raster layer from package raster including complete CRS information In order to set the extent (required to define a correct raster layer), we have to find out the extent in the units of the new projection We create an SpPolDF and then use spTransform In order to avoid any confusion introduced by overlaying within R, we overlay a reference vector layer on the raster layers exported as GTiff in QGIS “Google PseudoMercator” (epsg:3857) The overlay of a reference vector layer on the raster layer from get_map() assuming “Google PseudoMercator” is correct “Lon,Lat WGS84” (epsg:4326) The overlay of a reference vector layer on the raster layer from get_map() assuming “Lon,Lat WGS84” is shifted Care must be taken with the different CRS of the grid and the coordinates for overlaying layers in a plot. In order to test the overlayig within R, we import the reference vector layer into R and create a reprojected version: Overlay using standard plots is correct if we assume the output of get_map() is in “Google PseudoMercator” projection and use the vector layer with the same CRS levelplot() requires a raster layer, thus the output of get_map() cannot be used directly: We use the raster version that we calculated above. The plot is shifted if we assume get_map() output is in “Longitude, Latitude”: The drawback of this solution is that we are missing the color in the raster layer, although making a pseudocolor version of the get_map() output is possible. TBD layer() does not work with raster layers, therefore cannot use the rgmaprgbGM object. layer() does not work with the output of get_map() directly either, but can use grid.raster(). grid.raster() requires some explicit geographic information: As above, both layers are shifted if we assume raster and vector layers are in “Longitude, Latitude”: The vector layer must be in “Google PseudoMercator” to match the output of get_map(), but note that using the geographic attributes directly provides a wrong georeferencing and, as a consequence, no overlay: The geographic attributes must be reprojected to “Google PseudoMercator” (epsg:3857) to be consistent with the grid: Here we take advantage of the extent of the raster layer calculated in section 1.3 Source.


Яндекс.Метрика Рейтинг@Mail.ru Free Web Counter
page counter
Last Modified: April 22, 2016 @ 6:08 pm