在Leaflet如何画热图-R
背景
在绘制地图时候，我们经常会用到热图，Density map，在ggplot2中可根据坐标产生相应的密度图，2d distribution， 那么在交互式地图中，如何制作Density Map，
本次文章，主要介绍如何在Leaflet中，如何绘制热图。
英国伦敦霍乱地图
在该例子中，我们使用英国伦敦霍乱的数据来展示，在Leaflet中绘制Density map，
约翰·斯诺（John Snow）于1854年制作了一张著名的地图，显示了伦敦苏活区霍乱疫情造成的死亡以及该地区水源的位置。通过这样溯源的方法，他发现某个水源周围有大量霍乱死亡病例，并根据此推断阻止了疫情的爆发。（现在空间流行病学起源）
数据来源：Download
1.1 读取数据
我们从shp文件中读取Cholera数据，然后转换成经纬度坐标。
然后在Leaflet上显示出来
# read data
library(sf)
library(rgdal)
library(leaflet)
setwd()# set with your directory
deaths = read_sf("./SnowGIS_SHP/Cholera_Deaths.shp") %>%
st_transform(crs = 4326) # change to 4326
df_deaths = st_coordinates(deaths) %>% as.data.frame()
# leaflet map
leaflet(df_deaths) %>%
addTiles() %>%
addCircles(df_deaths$X,df_deaths$Y, radius = 0.5,opacity=0.6,col='blue')
image.png
1.2 点生成热图
这里我们主要利用的一个函数是bkde2D，将点转换成密度数据，然后
使用contourLines，将生成的2D转成polygons。
我们看一下效果。
# heat map
X=cbind(lng,lat)
kde2d = bkde2D(X, bandwidth=c(bw.ucv(X[,1]),bw.ucv(X[,2])),gridsize = c(100,100))
# to contourlines
CL = contourLines(kde2d$x1 , kde2d$x2 , kde2d$fhat)
## EXTRACT CONTOUR LINE LEVELS
LEVS = as.factor(sapply(CL, `[[`, "level"))
NLEV = length(levels(LEVS))
## CONVERT CONTOUR LINES TO POLYGONS
pgons = lapply(1:length(CL), function(i)
Polygons(list(Polygon(cbind(CL[[i]]$x, CL[[i]]$y))), ID=i))
spgons = SpatialPolygons(pgons)
## Leaflet map with polygons
leaflet(spgons) %>% addTiles() %>%
addPolygons(color = heat.colors(NLEV, NULL)[LEVS]) %>%
addCircles(df_deaths$X,df_deaths$Y, radius = 0.5,opacity=0.6,col='blue') 
image.png
1.3 密度图转换成Raster
上述的结果可以看到，有10层轮廓，我们进一步优化。将contourLines转换成Raster。
# Create Raster from Kernel Density output
KernelDensityRaster = raster(list(x=kde2d$x1 ,y=kde2d$x2 ,z = kde2d$fhat))
#create pal function for coloring the raster
palRaster = colorNumeric("Spectral", domain = KernelDensityRaster@data@values)
## Leaflet map with raster
leaflet() %>% addTiles() %>%
addRasterImage(KernelDensityRaster,
colors = palRaster,
opacity = .8) %>%
addLegend(pal = palRaster,
values = KernelDensityRaster@data@values,
title = "Kernel Density of Points")
image.png
1.4去除Raster边框
上面的图，我们看到边缘值很小，我们设定一个阈值。去除边缘值。 然后增加图例。根据value值分成不同颜色段。
#set low density cells as NA so we can make them transparent with the colorNumeric function
KernelDensityRaster@data@values[which(KernelDensityRaster@data@values < 1000)] = NA
#create pal function for coloring the raster
palRaster = colorNumeric("Spectral", domain = KernelDensityRaster@data@values, na.color = "transparent")
## Redraw the map
leaflet() %>% addTiles() %>%
addRasterImage(KernelDensityRaster,
colors = palRaster,
opacity = .8) %>%
addLegend(pal = palRaster,
values = KernelDensityRaster@data@values,
title = "Kernel Density of Points")
# Set legend
palRaster = colorBin("Spectral", bins = 7, domain = KernelDensityRaster@data@values, na.color = "transparent")
## Leaflet map with raster
leaflet() %>% addTiles() %>%
addRasterImage(KernelDensityRaster,
colors = palRaster,
opacity = .8) %>%
addLegend(pal = palRaster,
values = KernelDensityRaster@data@values,
title = "Kernel Density of Points")
image.png
image.png
1.5 提高热图像素
為了使其更平滑，只需在bkde2D函數中增加gridsize即可。這樣可以提高生成的柵格的分辨率。
## more smooth
X=cbind(lng,lat)
kde2d = bkde2D(X, bandwidth=c(bw.ucv(X[,1]),bw.ucv(X[,2])),
gridsize = c(1000,1000))
# Create Raster from Kernel Density output
KernelDensityRaster = raster(list(x=kde2d$x1 ,y=kde2d$x2 ,z = kde2d$fhat))
palRaster = colorBin("Spectral", bins = 7, domain = KernelDensityRaster@data@values, na.color = "transparent")
## Leaflet map with raster
leaflet() %>% addTiles() %>%
addRasterImage(KernelDensityRaster,
colors = palRaster,
opacity = .8) %>%
addLegend(pal = palRaster,
values = KernelDensityRaster@data@values,
title = "Kernel Density of Points")
image.png
结语
Leaflet为交互式地图提供了极大的方便，但是这里产生的热图，只是根据经纬度生成的。如何根据点上对应的value值，生成热图。这是很关键的一步。
参考
1.2d density plot with ggplot2 2.JOHN SNOW’S CHOLERA 3.How to build heatmap with the leaflet package
腾讯云开发者
