Как я могу нарисовать вручную рассчитанные плотности на карте с помощью ggmap?

Question

Я хочу построить тепловую карту, показывающую близость автобусных остановок в Чикаго, используя расстояние Манхэттена (L1) вместо евклидова и обратно-квадратный вес, используя среднесуточное пассажирское судно в качестве масштабирующего фактора.

Для загрузки карты, я использую ggmap

require(ggmap) 
chicago_map = get_map(location = c(lon=-87.64,lat=41.8787),zoom=14) 

Вот пример работы с:

require(foreach) 
require(geosphere) 
bounding_box = bb = list(c(41.96,-87.75),c(41.80,-87.6)) 
meter_per_lat = distVincentyEllipsoid(rev(bb[[1]]), 
             rev(bb[[1]] + c(1,0))) 
meter_per_lon = distVincentyEllipsoid(rev(bb[[1]]), 
             rev(bb[[1]] + c(0,1))) 
dist_converter = c(meter_per_lon, meter_per_lat) 
dist_manhattan <- function(p1, p2){ 
    #assume 20 meters is the smallest distance between points 
    pmax(20,abs(sweep(p2,2, p1)) %*% dist_converter) 
} 
set.seed(0) 
sample_data = data.frame(longitude=runif(20,bb[[1]][2],bb[[2]][2]), 
         latitude =runif(20,bb[[2]][1],bb[[1]][1]), 
         mean_riders=4*rpois(20,500)) 
chicago_grid = expand.grid(latitude=seq(bb[[1]][1], bb[[2]][1],length.out=200), 
          longitude=seq(bb[[1]][2],bb[[2]][2],length.out=400)) 
mean_values = foreach(i=1:20, .combine='+') %do% { 
    row = sample_data[i,] 
    dists = dist_manhattan(c(row$longitude,row$latitude), 
      as.matrix(chicago_grid[,c('longitude','latitude')])) 
    return(row$mean_riders/dists^2) 

}

chicago_data = chicago_grid 
chicago_data$rider_weight = mean_values 

И для черчения, я приблизительно с geom_point()

ggmap(chicago_map, extent='device') + geom_point(data=chicago_data, 
     aes(x=longitude,y=latitude,color=rider_weight,alpha=0.2,shape='15',size=4)) + 
     scale_alpha_identity() + scale_color_gradient(low='blue',high='red') + 
     scale_size_identity() + guides(shape='none') + ggtitle('Chicago Station Example Map') 

Хотя вы можете видеть части «тепловой карты» несколько четко, это, очевидно, не оптимальное решение.

Если я пытаюсь использовать geom_tile, я могу получить ладно выглядящий карту, но это занимает гораздо больше времени, чтобы произвести (что не желательно)

ggmap(chicago_map, extent='device') + geom_tile(data=chicago_data, 
     aes(x=longitude,y=latitude,fill=rider_weight,alpha=0.6)) + 
     scale_alpha_identity() + scale_fill_gradient(low='blue',high='red') + 
     scale_size_identity() + guides(shape='none') + ggtitle('Chicago Station Example Map') 

Я также могу заменить geom_tile на geom_raster, но geom_raster не работает вне декартовых координат. В частности, ошибка,

Ошибка: geom_raster работает только с декартовыми координатами

Есть ли лучший способ выполнить эту задачу?

Answer 1

я нашел еще один способ с ggmap превратить свои результаты в растре в нем

ggmap(chicago_map) + coord_cartesian() +geom_raster(data=chicago_data, aes(x=longitude,y=latitude,fill=rider_weight,alpha=0.2)) 

Это не выглядит очень отличается от geom_tile хотя

Answer 2

Я бы на самом деле превратить его в растр для того, чтобы сделать это, первый или обработка кадра данных в объект через пространственные точки

library(sp) 
library(rgdal) 
coordinates(chicago_data)=~longitude+latitude 

Дайте ему правильную проекцию

proj4string(chicago_data)=CRS("+init=epsg:4326") 
chicago_data = spTransform(chicago_data,CRS("+proj=longlat +datum=WGS84")) 
gridded(chicago_data) = TRUE 

И поворот это в растровый

r = raster(chicago_data) 
projection(r) = CRS("+proj=longlat +datum=WGS84") 

Я не уверен, что вы можете использовать это в ggplot2, но вы можете построить его как растр

plot(r)