快速计算点密度的度量并将其显示在地图上通常很有用。在本教程中，我们将使用 ggmap R 包中包含的德克萨斯州休斯顿的犯罪数据来演示这一点。
目标
- 计算点的二维空间密度
- 用 ggplot2 绘制密度表面
我们将从加载库开始。请注意，由于 Google 提供地图的方式发生了变化，本课程中不再使用 ggmap 包来生成底图，但本教程中使用的数据包含在 ggmap 包中。
library(ggplot2) library(ggmap)
然后，我们可以加载德克萨斯州休斯顿的内置犯罪数据集。
data(crime) # 删除任何有缺失数据的行 crime <- crime[complete.cases(crime), ] # 看一下犯罪数据的结构 str(crime) ## 'data.frame': 81803 obs. of 17 variables: ## $ time : POSIXct, format: "2010-01-01 06:00:00" "2010-01-01 06:00:00" ... ## $ date : chr "1/1/2010" "1/1/2010" "1/1/2010" "1/1/2010" ... ## $ hour : int 0 0 0 0 0 0 0 0 0 0 ... ## $ premise : chr "18A" "13R" "20R" "20R" ... ## $ offense : Factor w/ 7 levels "aggravated assault",..: 4 6 1 1 1 3 3 3 3 3 ... ## $ beat : chr "15E30" "13D10" "16E20" "2A30" ... ## $ block : chr "9600-9699" "4700-4799" "5000-5099" "1000-1099" ... ## $ street : chr "marlive" "telephone" "wickview" "ashland" ... ## $ type : chr "ln" "rd" "ln" "st" ... ## $ suffix : chr "-" "-" "-" "-" ... ## $ number : int 1 1 1 1 1 1 1 1 1 1 ... ## $ month : Ord.factor w/ 8 levels "january"<"february"<..: 1 1 1 1 1 1 1 1 1 1 ... ## $ day : Ord.factor w/ 7 levels "monday"<"tuesday"<..: 5 5 5 5 5 5 5 5 5 5 ... ## $ location: chr "apartment parking lot" "road / street / sidewalk" "residence / house" "residence / house" ... ## $ address : chr "9650 marlive ln" "4750 telephone rd" "5050 wickview ln" "1050 ashland st" ... ## $ lon : num -95.4 -95.3 -95.5 -95.4 -95.4 ... ## $ lat : num 29.7 29.7 29.6 29.8 29.7 ...
让我们用 ggplot2 绘制犯罪地点。
ggplot(crime, aes(x = lon, y = lat)) + geom_point() + coord_equal() + xlab('Longitude') + ylab('Latitude')
似乎有相当多的过度绘图。
让我们来绘制一个密度估计。计算密度的方法有很多种，如果密度估计的机制对您的应用程序很重要，那么研究专门用于点模式分析的软件包（例如spatstat）是值得的。另一方面，如果为了探索性数据分析的目的，您正在寻找快速而肮脏的实现，您还可以使用 ggplot's stat_density2d，它MASS::kde2d在后端使用二元正态核来估计密度。
ggplot(crime, aes(x = lon, y = lat)) + coord_equal() + xlab('Longitude') + ylab('Latitude') + stat_density2d(aes(fill = ..level..), alpha = .5, geom = "polygon", data = crime) + scale_fill_viridis_c() + theme(legend.position = 'none')
您可以通过对kde2d的调用传递参数stat_density2d。在这种情况下，我们改变参数h，它是与密度估计的空间范围或平滑度相关的带宽参数。
ggplot(crime, aes(x = lon, y = lat)) + coord_equal() + xlab('Longitude') + ylab('Latitude') + stat_density2d(aes(fill = ..level..), alpha = .5, h = .02, n = 300, geom = "polygon", data = crime) + scale_fill_viridis_c() + theme(legend.position = 'none')
作为替代方案，我们可能会考虑使用 alpha 透明度绘制原始数据点，以便我们可以看到实际数据，而不仅仅是数据模型。我们还将设置坐标以用作限制以专注于休斯顿市中心。
ggplot(crime, aes(x = lon, y = lat)) + geom_point(size = 0.1, alpha = 0.05) + coord_equal() + xlab('Longitude') + ylab('Latitude') + coord_cartesian(xlim = c(-95.1, -95.7), ylim = c(29.5, 30.1))
