Google Earth Engine——在 R 中计算和绘制二维空间点密度

快速计算点密度的度量并将其显示在地图上通常很有用。在本教程中,我们将使用 ggmap R 包中包含的德克萨斯州休斯顿的犯罪数据来演示这一点。

我们将从加载库开始。请注意,由于 Google 提供地图的方式发生了变化,本课程中不再使用 ggmap 包来生成底图,但本教程中使用的数据包含在 ggmap 包中。

代码语言:javascript
复制
library(ggplot2)
library(ggmap)

然后,我们可以加载德克萨斯州休斯顿的内置犯罪数据集。

代码语言:javascript
复制
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 绘制犯罪地点。

代码语言:javascript
复制
ggplot(crime, aes(x = lon, y = lat)) +

geom_point() +
coord_equal() +
xlab('Longitude') +
ylab('Latitude')

块基本图的图

似乎有相当多的过度绘图。

让我们来绘制一个密度估计。计算密度的方法有很多种,如果密度估计的机制对您的应用程序很重要,那么研究专门用于点模式分析的软件包(例如spatstat)是值得的。另一方面,如果为了探索性数据分析的目的,您正在寻找快速而肮脏的实现,您还可以使用 ggplot's stat_density2d,它MASS::kde2d在后端使用二元正态核来估计密度。

代码语言:javascript
复制
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')
块 stat_density2d 的图

您可以通过对kde2d的调用传递参数stat_density2d。在这种情况下,我们改变参数h,它是与密度估计的空间范围或平滑度相关的带宽参数。

代码语言:javascript
复制
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')
块 stat_density2d-bw 的图

作为替代方案,我们可能会考虑使用 alpha 透明度绘制原始数据点,以便我们可以看到实际数据,而不仅仅是数据模型。我们还将设置坐标以用作限制以专注于休斯顿市中心。

代码语言:javascript
复制
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))
块小点图