背景

在绘制地图时候，我们经常会用到热图，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

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")