projection of ggmap:get_map() output

On 05/06/2014 05:33 PM, Agustin Lobo wrote:

I think I've clarified this question here:
http://rpubs.com/alobo/getmapCRS

In short, we can get exact positioning of both layers with
* plot(r); plot(v, add=TRUE)
* spplot() + layer() and
* levelplot() + layer()
but we must be careful of CRS details.

Please send comments and corrections.

Useful!

These rpubs read nicely, but are very discouraging to try a rerun. Did
you upload the R-markdown file somewhere?

Anyway, I cannot see how the script can run when levelplot is called
before lattice is loaded, and also the "+" and layer() come, I believe
from latticeExtra which has not been loaded up to that point (2.2).

Does the "+" operator in

levelplot() + layer()

or

spplot() + layer()

check correspondence of the proj4string of the operands on both sides of
the + , can it be made to do so, or do we have remember your conclusions?

Agus

On Tue, Apr 29, 2014 at 11:40 AM, Oscar Perpi?an <oscar.perpinan at upm.es> wrote:

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)

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