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

[Oscar Perpiñan]

Hello,

I couldn't run your code because I am not able to read the file.
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).

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

In order to show that this is not an issue strictly related to
rasterVis::levelplot, I have written some code with spplot. I think
that the concordance between the ggmap raster and the shapefile is
aceptable.

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)