如何使用 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')

虽然您可以看到 "heatmap" 的部分内容比较清楚，但这显然不是最佳解决方案。

如果我尝试使用 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

我实际上会将其转换为栅格，为此，首先将您的数据框转换为空间点对象

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)

Answer 2

我找到了另一种使用 ggmap 将结果转换为栅格的方法

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

虽然

看起来与geom_tile没有太大区别

如何使用 ggmap 将手动计算的密度绘制到地图上？

How can I plot manually-calculated densities onto a map with ggmap?

plot

r

geospatial

heatmap

ggplot2