r 热图 - stat_density2d（ggmap）与 addHeatmap（闪亮的传单）

Question

我用 library(ggmap) 和 stat_density2d() 函数制作了静态热图。想要在动态 leaflet 地图上的闪亮应用程序中重新创建它，我发现了 addHeatmap()。然而，生成的图像是不同的，ggmap 版本似乎提供了正确的结果。

GGMAP

传单

造成这种差异的原因是什么？

对于运行下面两个可重现的例子，你可以下载我放在这里的一些数据（csv文件）。 https://drive.google.com/drive/folders/0B8_GTHBuoKSRR1VIRmhOUTJKYU0?usp=sharing

请注意，leaflet 结果因缩放级别而异，但与 ggmap 结果不匹配（例如，在最大热量位置方面）。

这是 ggmap 代码。

library(ggmap)
data <- read.csv("DATA.csv", sep=";")
data <- subset(data, !is.na(CrdLatDeg))
xmin <- min(data$CrdLonDeg)
xmax <- max(data$CrdLonDeg)
ymin <- min(data$CrdLatDeg)
ymax <- max(data$CrdLatDeg)
lon <- c(xmin,xmax)
lat <- c(ymin,ymax)
map <- get_map(location = c(lon = mean(lon), lat = mean(lat)), zoom = 17,
               maptype = "satellite", source = "google")
ggmap(map) + 
  labs(x="longitude", y="latitude") + 
  stat_density2d(data=data, aes(x=CrdLonDeg, y=CrdLatDeg, alpha= ..level.., fill= ..level..), colour=FALSE,
                 geom="polygon", bins=100) + 
  scale_fill_gradientn(colours=c(rev(rainbow(100, start=0, end=.7)))) + scale_alpha(range=c(0,.8)) + 
  guides(alpha=FALSE,fill=FALSE)

这是 leaflet 代码。

library(leaflet.extras)
data <- read.csv("DATA.csv", sep=";")
data <- subset(data, !is.na(CrdLatDeg))
leaflet(data) %>%
  addTiles(group="OSM") %>%
  addHeatmap(group="heat", lng=~CrdLonDeg, lat=~CrdLatDeg, max=.6, blur = 60)

Answer 1

两者使用不同的算法。您需要调整 addHeatmap 的 radius 和 blur 参数以及 stat_density2d 的 h 参数以获得有些相似的结果。

Answer 2

由于算法不同，图像看起来不同。

stat_density2d() 从离散数据中推断出 probability density function。

热图的 Leaflet 实现依赖于像 simpleheat, heatmap.js or webgl-heatmap 这样的库。这些热图不依赖于概率密度。（我不完全确定 r-leaflet 的 addHeatmap 使用了其中的哪一个）。

相反，这些热图的工作原理是为每个点绘制一个模糊的圆圈，将每个像素的值提高一个与点的强度成正比的量（在您的情况下为常数），并与点之间的距离成反比点和圆。每个数据点在热图中都显示为一个圆圈。您可以通过在 the heatmap.js webpage 中移动鼠标光标或通过查看图像右上角的这个孤点来看到这一点：

将热图想象成函数的可视化

f(pixel) = ∑ ( max( 0, radius - distance(pixel, point) ) · intensity(point) )

可以调整热图的半径和强度，但结果永远不会与统计密度估计相同。

Answer 3

我在 GIS 上找到了 this answer，我试图创建一个函数并将其应用到这个案例中。我还没有想出如何微调颜色渐变方案，但它似乎是一个很好的开始：

library(leaflet)
library(rlang)

addHeatMap <- function(data, lon, lat, intensity, show.legend, ...) {
  df <- data.table::as.data.table(data)
  df_expanded <- dplyr::slice(df, rep(1:dplyr::n(), times = !! enquo(intensity)))

  lon_var <- dplyr::pull(df_expanded, !! enquo(lon))
  lat_var <- dplyr::pull(df_expanded, !! enquo(lat))

  lon_bw <- MASS::bandwidth.nrd(lon_var)
  lat_bw <- MASS::bandwidth.nrd(lat_var)

  lon_lat_df <- dplyr::select(df_expanded, !! enquo(lon), !! enquo(lat))

  kde <- KernSmooth::bkde2D(lon_lat_df, bandwidth = c(lon_bw, lat_bw))
  CL <- contourLines(kde$x1 , kde$x2 , kde$fhat)
  LEVS <- as.factor(sapply(CL, `[[`, "level"))
  NLEV <- nlevels(LEVS)
  pgons <- lapply(1:length(CL), function(i)
  sp::Polygons(list(sp::Polygon(cbind(CL[[i]]$x, CL[[i]]$y))), ID = i))
  spgons <- sp::SpatialPolygons(pgons)

  if (show.legend) {
    leaflet::addPolygons(data = spgons, color = heat.colors(NLEV, NULL)[LEVS], stroke = FALSE, ...) %>% 
      leaflet::addLegend(colors = heat.colors(NLEV, NULL)[LEVS], labels = LEVS)
    } else {
    leaflet::addPolygons(data = spgons, color = heat.colors(NLEV, NULL)[LEVS], stroke = FALSE, ...)
    }
}

mydata <- read.csv("DATA.csv", sep=";")
mydata <- subset(mydata, !is.na(CrdLatDeg))

leaflet() %>%
  addTiles(group = "OSM") %>%
  addHeatMap(data = mydata, lon = CrdLonDeg, lat = CrdLatDeg, intensity = FsmIdf, show.legend = TRUE)

r 热图 - stat_density2d（ggmap）与 addHeatmap（闪亮的传单）

r heatmap - stat_density2d (ggmap) vs. addHeatmap (shiny leaflet)

r

heatmap

leaflet

ggmap

shiny