如何使用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_raster代替geom_tile，但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)