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

Oscar Perpiñan oscar.perpinan at upm.es
Tue Apr 29 11:40:03 CEST 2014


>> I couldn't run your code because I am not able to read the file.
> Note dsn is the directory where you have the 4 files and layer is the
> name of the layer, which is the name of the files without extension.

My fault, sorry.

> There has to be a reprojection on the fly because the the two layers of
> the plot are on different coordinate systems. I think it is clear now that
> ggplot uses the Google "PseudoMercator" for the matrix while expects
> all user input on epsg4326

I am not sure what do you mean with "reprojection on the fly" but
"layer" does not make any reprojection.
In fact, it is only a wrapper around
"eval(substitute(expression(...)))" to modify easily the panel
function of a trellis graphic.
However, in my code I am using "grid.raster" with its arguments width,
height, x and y, to resize the ggmap raster and to position it in the
graphic window created with spplot.

> I think "acceptable" depends on the scale only. At coarser scale you
> just do not see the error
> because it is relative to the extent and the result converges to the
> correct  reprojection.

I think it depends on the extent. I have tested the code again with
the administrative boundaries of Spain and some of its regions
(http://biogeo.ucdavis.edu/data/gadm2/shp/ESP_adm.zip).
It fails for the country extent (as Edzer suggested and you have found
with your file) but it works correctly for a smaller extent (code
copied below).

> In the code you provide, how do you know NAD27 is the appropriate datum?
> There is no prj file (one of the big shortcomings of shapefiles is
> that the presence
> of crs information in the form of a prj file is not compulsory). Note
> that you define
> a datum != WGS84 and you do not define the towgs argument.

It is an example extracted from the p.121 of the 2008 edition of
"Applied Spatial Data Analysis with R". That chapter is now available
here: http://cran.r-project.org/web/packages/maptools/vignettes/combine_maptools.pdf

Best,

Oscar.

#######################################################
library(ggmap)
library(latticeExtra)
library(sp)
library(maptools)

## http://biogeo.ucdavis.edu/data/gadm2/shp/ESP_adm.zip
setwd('/tmp')
spain <- readShapePoly('ESP_adm/ESP_adm1')
proj4string(spain) <- CRS("+proj=longlat +datum=WGS84")

andalucia <- spain[spain$ID_1 == 1,]
catalonia <- spain[spain$ID_1 == 6,]

## myPoly <- spain
myPoly <- andalucia
## myPoly <- catalonia

bbPoly <- bbox(myPoly)
gmap <- get_map(c(bbPoly), maptype='watercolor', source='stamen', crop=FALSE)
bbMap <- attr(gmap, 'bb')
latCenter <- with(bbMap, ll.lat + ur.lat)/2
lonCenter <- with(bbMap, ll.lon + ur.lon)/2
height <- with(bbMap, ur.lat - ll.lat)
width <- with(bbMap, ur.lon - ll.lon)

spplot(myPoly["ID_1"], fill='transparent') +
                  layer(grid.raster(gmap,
                                    x=lonCenter, y=latCenter,
                                    width=width, height=height,
                                    default.units='native'),
                        under=TRUE)



More information about the R-sig-Geo mailing list