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

[Agustin Lobo]

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.
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)
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo