Etiquetas
ggmap, ggplot2, google map, maps, openstreetmap, photovoltaics, proportional symbol, R, spplot, stamen
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 (automap
and 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
ggmap
but a similar approach can be followed with theOpenStreetMap
package). - 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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(sp) | |
library(ggmap) | |
library(RColorBrewer) | |
## Get data from OpenPV. Open this URL: | |
## https://openpv.nrel.gov/search?&state=CA&minsize=100&maxsize=30600&pagenum=1&nppage=25 | |
## 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)[3] <- '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) | |
## http://spatialreference.org/ref/sr-org/6864/ | |
## 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, | |
x=lonCenter, y=latCenter, | |
width=width, height=height, | |
default.units='native') | |
## 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 | |
spplot(caPV["Pdc.kW"], | |
cuts = breaks, | |
col.regions=brewer.pal(n=5, 'Greens'), | |
cex=sizes, | |
edge.col='black', alpha=0.7, | |
scales=list(draw=TRUE), key.space='right', | |
sp.layout=sp.raster) | |
Pingback: Mapping Flows in R … with data.table and lattice | Omnia sunt Communia!
Pingback: Google Maps and ggmap « Software for Exploratory Data Analysis and Statistical Modelling - Statistical Modelling with R
Jannes dijo:
Hola Oscar,
thank you for this contribution! Do you have any idea, why the same script does not work with spatial polygons (instead of spatial points as you have used)? It seems the polygons are plotted behind the sp.raster.
Thank you!
Cheers,
Jannes
Oscar Perpiñán Lamigueiro dijo:
Hola!
The behavior of
sp.layout
is different if you use aSpatialPolygonsDataFrame
. From the help page ofspplot
:There is a
first
argument to change this default behaviour, but I can get it work. Instead, you can use the+.trellis
andlayer
mechanisms defined in thelatticeExtra
package.I have created a new gist with an example: https://gist.github.com/oscarperpinan/7482848.
Jannes dijo:
Hola Oscar,
perfecto, muchas gracias por tu ayuda, la aprecio mucho!
Que estés bien,
Jannes
Roland dijo:
Hola Oscar,
very nice code (https://gist.github.com/oscarperpinan/7482848), especially the latticeExtra section was very helpful to me.
But may I ask a question…
If I put that stuff inside a function, even if I use print(spplot…), it does not work. It says that the object gmap is not found. Do you have any idea what is the problem and how to solve it?
Thank you and greetings!
(Sorry, I Think here is the better place to share this question than to post on gist, or what do ypu think?)
Oscar Perpiñán Lamigueiro dijo:
Hola Roland,
Thanks for your comments! Glad to know this code helps you.
I suspect your problem is related to the
latticeExtra::layer
mechanism. From its help page:I will post a new gist to show a suitable solution.
PD. Both the blog and gist are good places for comments.
Roland dijo:
Hola Oscar,
thank you very very much for your help and the gist post. And sorry for bothering you again….
I’m able to reproduce your example.
But I’m still having some problem and I’m not sure if I understand it in the right way…
a) The following does not work out, even if I pass the arguments (e.g. gmap) to the ‘data’ argument of ‘layer’:
f1 <- function(){
SIDS <- …
gmap <- get_map(…)
bbMap gmap is missing again.
b) Moreover it seems not to be possible to create the data which should be plotted inside the function (withLayerInFunction()) from arguments given to that function. E.g. something like that:
withLayerInFunction <- function(a,b){
sids <- "create data from a and b"
…
}
My conclusions: If I pass arguments to the withLayerInFunction() function, these arguments must be known globaly ("you can usually refer to objects which exist in the global environment") ?!?
Or everything has to be defined inside the print(spplot…) resp. the layer() function ("you typically refer to variables as if inside the panel function") ?!?
Is that right?
And what can I do if I need something like example a) or b) (besides reorganizing all my code) !?
Would be the "<<-" Operator helpful to create global variables, that might be found by the function layer() ?
Thank you very much again and greetings!
Roland, little bit confused…
Oscar Perpiñán Lamigueiro dijo:
Hi!
I have added more versions of the function in the gist. The last one uses external arguments. I think it is what you need.
Cheers,
Oscar.
Roland dijo:
Hola Oscar!
Muchas gracias! And thank you for your patience!
I had an (stupid) error in the «data» argument in «layer».
Thanks again for your help!
Cheers!
Roland
John Sampson dijo:
Thank you for this script. I am new to R and would like to create a historical map of where ancestors in my family have lived. I tried to re-create your map but ran into a problem with the bounding box. I think I may just not know what to do with the url you provided in the comments. this line of code:
> latCenter <- with(bbMap, 11.lat + ur.lat)/2
produces this error:
Error: unexpected symbol in "latCenter <- with(bbMap, 11.lat"
Would you be able to provide a suggestion?
Thank you so much!
John
Oscar Perpiñán Lamigueiro dijo:
Hi John,
The error is because you are using 11 (number) instead of ll (two l’s). I have checked the code I posted at github and it is correct.
On the other hand, you have to open the url and export the result as a CSV file. I will add some comments to clarify this point.
Thanks for your comments.
Oscar.
Pingback: Stamen maps with spplot | R for Journalists | S...
Pingback: Stamen maps with spplot | Mapas con R | Scoop.it