背景
在绘制地图时候,我们经常会用到热图,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')
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')
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")
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")
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")
结语
Leaflet为交互式地图提供了极大的方便,但是这里产生的热图,只是根据经纬度生成的。如何根据点上对应的value值,生成热图。这是很关键的一步。
参考
1.2d density plot with ggplot2
2.JOHN SNOW’S CHOLERA
3.How to build heatmap with the leaflet package