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

[Agustin Lobo]

Oscar,

On Mon, Apr 28, 2014 at 11:58 AM, Oscar Perpiñan <oscar.perpinan at upm.es> wrote:
> Hello,
>
> 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.

> Anyway, some comments about your last email:
>
> - "layer" is a function of the latticeExtra package. This function
> modifies the panel function of a "trellis" graphic produced with a
> function of the lattice package. This function, together with
> "+.trellis" (of the same package) allows for overlaying layers. It is
> important to note that it does not carry out any projection. It only
> displays the graphical result of another function (for example,
> sp.polygons or whatever).

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

> - You create the canvas with "levelplot" displaying the ggmap raster,
> and then overlay the data with layer. This is the same approach that
> ggmap uses. Instead, I prefer to rely on the main data to be displayed
> to configure the graphical window instead of using an auxiliary image
> as the starting point.
>
> - I'm not a fan of the ggmap approach. I agree with you that the
> resulting object should include the projection information and even
> adhere to the sp classes. However, I think that their results are
> useful if managed with care.

Users cannot handle with care undocumented idiosyncratic features.

> In order to show that this is not an issue strictly related to
> rasterVis::levelplot,
I agree this is not an specific problem of reasterVis:levleplot.

> I have written some code with spplot. I think
> that the concordance between the ggmap raster and the shapefile is
> aceptable.

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 do not think we should settle for this kind of errors.

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.

Agus

> Best,
>
> Oscar.
>
> ##########################################################
> library(sp)
> library(ggmap)
> ## latticeExtra must be loaded after ggmap because both ggplot2 and
> ## latticeExtra define a 'layer' function. We need the definition from
> ## latticeExtra.
> library(latticeExtra)
>
> ## Shape example
> library(maptools)
> owd <- getwd()
> setwd(system.file("shapes", package = "maptools"))
> sc90 <- readShapeSpatial("co45_d90")
> proj4string(sc90) <- CRS("+proj=longlat +datum=NAD27")
> setwd(owd)
> ## Get the stamen map
> bbPoly <- bbox(sc90)
> gmap <- get_map(c(bbPoly), maptype='toner', 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)
>
> ## Use latticeExtra::+.trellis and latticeExtra::layer
> spplot(sc90['AREA'], fill='transparent', col='indianred1') +
>                   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