stamen.R

library(sp)
library(ggmap)
library(RColorBrewer)

## Get data from OpenPV
## https://openpv.nrel.gov/search?&state=CA&minsize=100&maxsize=30600&pagenum=1&nppage=25
caPV <- read.csv('californiaOpenPV.csv')
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 function 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)