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

[Agustin Lobo]

Edzer,

Thanks for the corrections. I had moved chunks and left loading the
packages in the wrong place.
It is correct now. Also, I've included a link to a file with just R
code and the link to the vector shapefile.

As far as I can tell, "levelplot() + layer()"
or "spplot() + layer()"
do not check correspondence of the proj4string of the operands on both sides of
the +. Actually, both commands totally ignore any CRS information.
This is why we've got into this mess
and must be aware of providing consistent layers.

Agus

On Tue, May 6, 2014 at 6:02 PM, Edzer Pebesma
<edzer.pebesma at uni-muenster.de> wrote:
>
>
> 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)
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> 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
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>