Custom geographic heatmap in R with ggmap

Question

The goal is to build something like http://rentheatmap.com/sanfrancisco.html

I got a map with ggmap and was able to draw points on top of it.

library('ggmap') map <- get_map(location=c(lon=20.46667, lat=44.81667), zoom=12, maptype='roadmap', color='bw') positions <- data.frame(lon=rnorm(100, mean=20.46667, sd=0.05), lat=rnorm(100, mean=44.81667, sd=0.05), price=rnorm(10, mean=1000, sd=300)) ggmap(map) + geom_point(data=positions, mapping=aes(lon, lat)) + stat_density2d(data=positions, mapping=aes(x=lon, y=lat, fill=..level..), geom="polygon", alpha=0.3) 

This is a good image based on density. Does anyone know how to do something similar, but uses the position $ property to draw contours and scale?

I looked carefully through stackoverflow.com and did not find a solution.

EDIT 1

 positions$price_cuts <- cut(positions$price, breaks=5) ggmap(map) + stat_density2d(data=positions, mapping=aes(x=lon, y=lat, fill=price_cuts), alpha=0.3, geom="polygon") 

Results of five independent stat_density statistics graphs:

EDIT 2 (from hrbrmstr )

 positions <- data.frame(lon=rnorm(10000, mean=20.46667, sd=0.05), lat=rnorm(10000, mean=44.81667, sd=0.05), price=rnorm(10, mean=1000, sd=300)) positions$price <- ((20.46667 - positions$lon) ^ 2 + (44.81667 - positions$lat) ^ 2) ^ 0.5 * 10000 positions <- data.frame(lon=rnorm(10000, mean=20.46667, sd=0.05), lat=rnorm(10000, mean=44.81667, sd=0.05)) positions$price <- ((20.46667 - positions$lon) ^ 2 + (44.81667 - positions$lat) ^ 2) ^ 0.5 * 10000 positions <- subset(positions, price < 1000) positions$price_cuts <- cut(positions$price, breaks=5) ggmap(map) + geom_hex(data=positions, aes(fill=price_cuts), alpha=0.3) 

Results in:

It also creates a decent picture on real data. This is the best result. Other suggestions are welcome.

EDIT 3: Here are the test data and the results of the method above:

https://raw.githubusercontent.com/artem-fedosov/share/master/kernel_smoothing_ggplot.csv

 test<-read.csv('test.csv') ggplot(data=test, aes(lon, lat, fill=price_cuts)) + stat_bin2d(, alpha=0.7) + geom_point() + scale_fill_brewer(palette="Blues") 

I believe that some method that uses a kernel other than a density kernel should be used to calculate regular polygons. It seems that the function should be in ggplot out of the box, but I cannot find it.

EDIT 4: I appreciate your time and efforts to figure out the right solution to this seemingly not too complicated question. I voted for your answers as a good approximation to the goal.

I found one problem: data with circles is too artificial, and approaches do not work so well for reading world data.

Paul's approach gave me the plot:

It seems to be capturing data patterns that are cool.

jazzurro approage gave me this plot:

He also received samples. However, both graphs do not seem as pretty as the stat_density2d graph by default. I still wait a couple of days to see if some other solution comes. If not, I will reward the generosity of jazzurro, as this will be the result that I will use.

There is an open version of python + google_maps code. Maybe someone will find inspiration here: https://github.com/jeffkaufman/apartment_prices

Answer 1

It seems to me that the map in the link you provided was created using interpolation. With that in mind, I wondered if I could achieve the same ascetic by superimposing an interpolated raster on ggmap.

 library(ggmap) library(akima) library(raster) ## data set-up from question map <- get_map(location=c(lon=20.46667, lat=44.81667), zoom=12, maptype='roadmap', color='bw') positions <- data.frame(lon=rnorm(10000, mean=20.46667, sd=0.05), lat=rnorm(10000, mean=44.81667, sd=0.05), price=rnorm(10, mean=1000, sd=300)) positions$price <- ((20.46667 - positions$lon) ^ 2 + (44.81667 - positions$lat) ^ 2) ^ 0.5 * 10000 positions <- data.frame(lon=rnorm(10000, mean=20.46667, sd=0.05), lat=rnorm(10000, mean=44.81667, sd=0.05)) positions$price <- ((20.46667 - positions$lon) ^ 2 + (44.81667 - positions$lat) ^ 2) ^ 0.5 * 10000 positions <- subset(positions, price < 1000) ## interpolate values using akima package and convert to raster r <- interp(positions$lon, positions$lat, positions$price, xo=seq(min(positions$lon), max(positions$lon), length=100), yo=seq(min(positions$lat), max(positions$lat), length=100)) r <- cut(raster(r), breaks=5) ## plot ggmap(map) + inset_raster(r, extent(r)@xmin, extent(r)@xmax, extent(r)@ymin, extent(r)@ymax) + geom_point(data=positions, mapping=aes(lon, lat), alpha=0.2) 

http://i.stack.imgur.com/qzqfu.png

Unfortunately, I could not figure out how to change color or alpha using inset_raster ... perhaps due to my ignorance with ggmap.

EDIT 1

This is a very interesting problem that makes me scratch my head. The interpolation didn’t look exactly like I thought it applied to real data; the polygon fits itself and jazzurro, of course, looks much better!

It's amazing why the raster approach looked so jagged, I looked at the map you attached a second time and noticed a visible buffer around the data points ... I wondered if I could use some rgeos tools to try to reproduce the effect:

 library(ggmap) library(raster) library(rgeos) library(gplots) ## data set-up from question dat <- read.csv("clipboard") # load real world data from your link dat$price_cuts <- NULL map <- get_map(location=c(lon=median(dat$lon), lat=median(dat$lat)), zoom=12, maptype='roadmap', color='bw') ## use rgeos to add buffer around points coordinates(dat) <- c("lon","lat") polys <- gBuffer(dat, byid=TRUE, width=0.005) ## calculate mean price in each circle polys <- aggregate(dat, polys, FUN=mean) ## rasterize polygons r <- raster(extent(polys), ncol=200, nrow=200) # define grid r <- rasterize(polys, r, polys$price, fun=mean) ## convert raster object to matrix, assign colors and plot mat <- as.matrix(r) colmat <- matrix(rich.colors(10, alpha=0.3)[cut(mat, 10)], nrow=nrow(mat), ncol=ncol(mat)) ggmap(map) + inset_raster(colmat, extent(r)@xmin, extent(r)@xmax, extent(r)@ymin, extent(r)@ymax) + geom_point(data=data.frame(dat), mapping=aes(lon, lat), alpha=0.1, cex=0.1) 

PS I found out that the color matrix needs to be sent to inset_raster to configure the overlay

Answer 2

Here is my approach. The geom_hex approach geom_hex good. When it came out, I really liked it. I'm still doing it. Since you asked for something else, I tried the following. I think my result is similar to the one with stat_density2d . But I could have avoided the problems you had. Basically, I created a shapefile and drew polygons. I multiplied price zone data (price_cuts) and drew polygons from the edge to the center of the zone. This approach is in line EDIT 1 and 2. I think there is still some distance to reach your final goal if you want to draw a map with a large area. But hopefully this will allow you to move forward. Finally, I'd like to say thanks to a few SO users who asked great questions related to polygons. I could not come up with this answer without them.

 library(dplyr) library(data.table) library(ggmap) library(sp) library(rgdal) library(ggplot2) library(RColorBrewer) ### Data set by the OP positions <- data.frame(lon=rnorm(10000, mean=20.46667, sd=0.05), lat=rnorm(10000, mean=44.81667, sd=0.05)) positions$price <- ((20.46667 - positions$lon) ^ 2 + (44.81667 - positions$lat) ^ 2) ^ 0.5 * 10000 positions <- subset(positions, price < 1000) ### Data arrangement positions$price_cuts <- cut(positions$price, breaks=5) positions$price_cuts <- as.character(as.integer(positions$price_cuts)) ### Create a copy for now ana <- positions ### Step 1: Get a map map <- get_map(location=c(lon=20.46667, lat=44.81667), zoom=11, maptype='roadmap', color='bw') ### Step 2: I need to create SpatialPolygonDataFrame using the original data. ### http://stackoverflow.com/questions/25606512/create-polygon-from-points-and-save-as-shapefile ### For each price zone, create a polygon, SpatialPolygonDataFrame, and convert it ### it data.frame for ggplot. cats <- list() for(i in unique(ana$price_cuts)){ foo <- ana %>% filter(price_cuts == i) %>% select(lon, lat) ch <- chull(foo) coords <- foo[c(ch, ch[1]), ] sp_poly <- SpatialPolygons(list(Polygons(list(Polygon(coords)), ID=1))) bob <- fortify(sp_poly) bob$area <- i cats[[i]] <- bob } cathy <- as.data.frame(rbindlist(cats)) ### Step 3: Draw a map ### The key thing may be that you subet data for each price_cuts and draw ### polygons from outer side given the following link. ### This link was great. This is exactly what I was thinking. ### http://stackoverflow.com/questions/21748852/choropleth-map-in-ggplot-with-polygons-that-have-holes ggmap(map) + geom_polygon(aes(x = long, y = lat, group = group, fill = as.numeric(area)), alpha = .3, data = subset(cathy, area == 5))+ geom_polygon(aes(x = long, y = lat, group = group, fill = as.numeric(area)), alpha = .3, data =subset(cathy, area == 4))+ geom_polygon(aes(x = long, y = lat, group = group, fill = as.numeric(area)), alpha = .3, data = subset(cathy, area == 3))+ geom_polygon(aes(x = long, y = lat, group = group, fill = as.numeric(area)), alpha = .3, data = subset(cathy, area == 2))+ geom_polygon(aes(x = long, y = lat, group = group, fill = as.numeric(area)), alpha= .3, data = subset(cathy, area == 1))+ geom_point(data = ana, aes(x = lon, y = lat), size = 0.3) + scale_fill_gradientn(colours = brewer.pal(5,"Spectral")) + scale_x_continuous(limits = c(20.35, 20.58), expand = c(0, 0)) + scale_y_continuous(limits = c(44.71, 44.93), expand = c(0, 0)) + guides(fill = guide_legend(title = "Property price zone"))