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

快速计算点密度的度量并将其显示在地图上通常很有用。在本教程中，我们将使用 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))