Several R packages provide an interface to query map services (Google Maps, Stamen Maps or OpenStreetMap) to obtain raster images from them. As far as I know, there are three packages devoted to this task: RgoogleMaps, OpenStreetMap and ggmap. The latter two are increasingly popular with a wide collection of providers.
Both of them use the background image to configure the graphical window defining functions (
ggmap, respectively) to display the images with ggplot2. Additional information can be layered upon this background image with ggplot2 functions. In my opinion, this approach is somewhat strange. I feel more comfortable with the approach implemented in the
spplot function of the
sp package, which relies on the main data to be displayed to configure the graphical window instead of using an auxiliary image as the starting point.
Although none of these two packages include functions to work with
spplot, it is not difficult to build a mixed solution. The key point is that the result of their queries is mainly a raster image which can be easily displayed with the
grid.raster function after some corrections.
The next code displays the location of the photovoltaic systems installed in California grouped by nominal power using data available from the OpenPV project. The procedure is:
- Download the data and define the coordinates and projection.
- Download an image from the Stamen Maps service according to the spatial extent of the data (this example uses
ggmapbut a similar approach can be followed with the
- Use grid.raster to display the image using the coordinates of its corners to define the width, height and center location.
- Display the variable with circles of different sizes and colours over the background image.
|## Get data from OpenPV. Open this URL:|
|## and export the results as a CSV file named "californiaOpenPV.csv"|
|caPV <- read.csv('californiaOpenPV.csv')|
|## Longitude and Latitude are the names of columns where the spatial information is stored.|
|## With the coordinates<- method caPV is now an SpatialPointsDataFrame|
|coordinates(caPV) <- ~ Longitude + Latitude|
|proj4string(caPV) <- CRS("+proj=longlat +datum=WGS84")|
|names(caPV) <- 'Pdc.kW'|
|## Download stamen tiles using the bounding box of the SpatialPointsDataFrame object|
|bbPoints <- bbox(caPV)|
|gmap <- get_map(c(bbPoints), maptype='watercolor', source='stamen', crop=FALSE)|
|## Bounding box of the map to resize and position the image with grid.raster|
|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 sp.layout of spplot: a list with the name of the function|
|## ('grid.raster') and its arguments|
|sp.raster <- list('grid.raster', gmap,|
|## Define classes and sizes of the circle for each class|
|breaks <- c(100, 200, 500, 1e3, 25e3)|
|classes <- cut(caPV$Pdc.kW, breaks)|
|meds <- tapply(caPV$Pdc.kW, classes, FUN=median)|
|sizes <- (meds/max(meds))^0.57 * 1.8|
|## Finally, the spplot function|
|cuts = breaks,|