Географическая тепловая карта настраиваемого свойства в R с помощью ggmap

Мне кажется, что карта в прикрепленной вами ссылке была создана с использованием интерполяции. Имея это в виду, я подумал, смогу ли я добиться такого же аскетизма, наложив интерполированный растр на 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) 

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

К сожалению, я не мог понять, как изменить цвет или альфа-канал с помощью inset_raster ... вероятно, из-за моего незнания ggmap.

ИЗМЕНИТЬ 1

Это очень интересная проблема, от которой я ломаю голову. Интерполяция выглядела не совсем так, как я думал, применительно к реальным данным; многоугольник приближается самостоятельно, и jazzurro, безусловно, выглядит намного лучше!

Интересно, почему растровый подход выглядит таким неровным, я еще раз взглянул на карту, которую вы прикрепили, и заметил очевидный буфер вокруг точек данных ... Я подумал, могу ли я использовать некоторые инструменты rgeos, чтобы попытаться воспроизвести эффект:

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) 

P.S. Я обнаружил, что матрицу цветов нужно отправить в inset_raster для настройки наложения

Вот мой подход. geom_hex подход хорош. Когда это вышло, мне это очень понравилось. Я все еще делаю. Поскольку вы спросили кое-что еще, я попробовал следующее. Думаю, мой результат похож на результат с stat_density2d. Но я мог избежать проблем, которые были у вас. По сути, я сам создал шейп-файл и нарисовал многоугольники. Я разбил данные по ценовым зонам (price_cuts) и нарисовал полигоны от края до центра зоны. Этот подход находится в линии EDIT 1 и 2. Я думаю, что до вашей конечной цели еще есть некоторое расстояние, если вы хотите нарисовать карту с большой площадью. Но я надеюсь, что это позволит вам двигаться вперед. Наконец, я хотел бы поблагодарить пару пользователей SO, которые задали отличные вопросы, связанные с полигонами. Без них я не смог бы прийти к такому ответу.

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