Geographische Heatmap einer benutzerdefinierten Eigenschaft in R mit ggmap

Question

Das Ziel ist, so etwas wie zu bauen http://rentheatmap.com/sanfrancisco.html

Ich habe Karte mit ggmap und in der Lage Punkte oben drauf plotten.

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) 

Dies ist ein schönes Bild auf Dichte basiert. Weiß jemand, wie man etwas macht, das gleich aussieht, aber die Position $ property verwendet, um Konturen und Skalierungen zu erstellen?

Ich schaute gründlich durch stackoverflow.com und fand keine Lösung.

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") 

Ergebnisse in fünf unabhängigen stat_density Stellplätze:

EDIT 2 (von 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) 

Ergebnisse in:

Es schafft ein anständiges Bild auf echten Daten. Dies ist das bisher beste Ergebnis. Weitere Vorschläge sind willkommen.

EDIT 3: Hier Testdaten und Ergebnisse eines Verfahrens oben:

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") 

Ich glaube, dass es einige Verfahren, die andere verwendet als Dichte Kernel, um richtige Polygone zu berechnen. Es scheint, dass das Feature in ggplot out of the box sein sollte, aber ich kann es nicht finden.

EDIT 4: Ich schätze Sie Zeit und Mühe, um die richtige Lösung für diese scheinbar nicht zu komplizierte Frage zu finden. Ich habe beide Antworten als gute Annäherung an das Ziel gewählt.

Ich habe ein Problem aufgedeckt: Die Daten mit Kreisen sind zu künstlich und die Ansätze funktionieren auf Lesendaten nicht so gut.

Paul Ansatz gab mir die Handlung:

Es scheint, dass es Muster der Daten erfasst, die kühl ist.

jazzurro des approage hat mir dieses Grundstück:

es auch die Muster bekam. Beide Plots scheinen jedoch nicht so schön zu sein wie das standardmäßige plot stat_density2d. Ich werde noch ein paar Tage warten, um zu sehen, ob eine andere Lösung auftauchen wird.Wenn nicht, werde ich das Bounty an Jazzurro vergeben, da dies das Ergebnis sein wird, das ich verwenden werde.

Es gibt eine offene Python + google_maps Version des erforderlichen Codes. Vielleicht wird jemand hier Inspiration finden: https://github.com/jeffkaufman/apartment_prices

Answer 1

Es sieht für mich Sie, wie die Karte in der Verbindung gebunden durch Interpolation erzeugt wurde. In diesem Sinne fragte ich mich, ob ich durch Überlagerung eines interpolierten Rasters auf einer ggmap einen ähnlichen Asketen erreichen könnte.

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

Leider konnte ich nicht herausfinden, wie die Farbe oder die Alpha ändern inset_raster mit ... wahrscheinlich wegen meiner mangelnder Vertrautheit mit ggmap.

EDIT 1

Dies ist ein sehr interessantes Problem, das ich meinen Kopf kratzen. Die Interpolation hatte nicht ganz das Aussehen, das ich dachte, wenn sie auf reale Daten angewandt würde; das polygon nähert sich von selbst und jazzurro sieht sicherlich viel besser aus!

Ich frage mich, warum der Raster-Ansatz so gezackt aussieht. Ich habe einen zweiten Blick auf die beigefügte Karte geworfen und einen scheinbaren Puffer um die Datenpunkte bemerkt ... Ich habe mich gefragt, ob ich einige RGEOS-Tools verwenden könnte, um den Effekt zu replizieren :

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 Ich fand heraus, dass eine Matrix von Farben an inset_raster geschickt werden muss, um das Overlay anzupassen

Answer 2

Hier ist mein Ansatz. Der Ansatz ist nett. Als das herauskam, habe ich es wirklich gemocht. Ich mache immernoch. Da du etwas mehr gefragt hast, habe ich folgendes versucht. Ich denke, mein Ergebnis ist ähnlich einem mit stat_density2d. Aber ich konnte die Probleme vermeiden, die Sie hatten. Ich habe im Grunde selbst ein Shapefile erstellt und Polygone gezeichnet. Ich untergliederte Daten nach Preiszonen (price_cuts) und zeichnete Polygone von der Mitte zur Mitte der Zone. Dieser Ansatz liegt in der Linie von EDIT 1 und 2. Ich denke, es ist noch ein wenig Abstand, um Ihr endgültiges Ziel zu erreichen, wenn Sie eine Karte mit einer großen Fläche zeichnen wollen. Aber ich hoffe, das wird dich vorwärts bringen. Abschließend möchte ich ein paar SO-Benutzern danken, die große Fragen zu Polygonen gestellt haben. Ich konnte diese Antwort ohne sie nicht finden.

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"))