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

[Edzer Pebesma]

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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20140506/7ede7b12/attachment.bin>