analyze-point-pattern.qmd

# 点模式数据分析 {#sec-analyze-point-pattern}

```{r}
#| echo: false

source("_common.R")
```

本章以斐济地震数据集 quakes 为例介绍空间点模式数据的操作、探索和分析。

```{r}
#| message: false

library(spatstat)
library(sf)
library(ggplot2)
```

[**spatstat**](https://github.com/spatstat/spatstat/) 是一个伞包，囊括 8 个子包，构成一套完整的空间点模式分析工具。

1.  spatstat.utils 基础的辅助分析函数
2.  spatstat.data 点模式分析用到的示例数据集
3.  spatstat.sparse 稀疏数组
4.  spatstat.geom 空间数据类和几何操作
5.  spatstat.random 生成空间随机模式
6.  spatstat.explore 空间数据的探索分析
7.  spatstat.model 空间数据的参数建模和推理
8.  spatstat.linnet 线性网络上的空间分析

**sf** 包是一个专门用于空间矢量数据操作的 R 包。**ggplot2** 包提供的几何图层函数 `geom_sf()` 和坐标参考系图层函数 `coord_sf()` 支持可视化空间点模式数据。

## 数据操作

### 类型转化

先对斐济地震数据 quakes 数据集做一些数据类型转化，从 data.frame 转 Simple feature 对象。

```{r}
library(sf)
quakes_sf <- st_as_sf(quakes, coords = c("long", "lat"), crs = st_crs(4326))
quakes_sf
```

### 坐标转化

如果知道两个投影坐标系的 EPSG 代码，输入坐标就可以完成转化。如将坐标系 `EPSG:4326` 下的坐标 $(2,49)$ 投影到另一个坐标系 `EPSG:3857` 。

```{r}
st_transform(
  x = st_sfc(st_point(x = c(2, 49)), crs = 4326), crs = 3857
)
```

| 名称  | EPSG | 赤道半径    | 半轴             | 发明者               |
|-------|------|-------------|------------------|----------------------|
| GRS80 | 3857 | a=6378137.0 | rf=298.257222101 | GRS 1980(IUGG, 1980) |
| WGS84 | 4326 | a=6378137.0 | rf=298.257223563 | WGS 84               |

函数 `st_crs()` 查看坐标参考系的信息，比如 EPSG 代码为 4326 对应的坐标参考系统信息。我们也可以通过[网站](https://epsg.io/3832)查询 EPSG 代码对应的坐标参考系统的详细介绍。

```{r}
st_crs("EPSG:4326")
```

地球看作一个椭球体 ELLIPSOID，长半轴 6378137 米，短半轴 298.257223563 米，椭圆形的两个轴，纬度单位 0.0174532925199433， 经度单位 0.0174532925199433 。

地球是一个不规则的球体，不同的坐标参考系对地球的抽象简化不同，会体现在坐标原点、长半轴、短半轴等属性上。为了方便在平面上展示地理信息，需要将地球表面投影到平面上，墨卡托投影是其中非常重要的一种投影方式，墨卡托投影的详细介绍见 [PROJ 网站](https://proj.org/operations/projections/merc.html) 。WGS 84 / Pseudo-Mercator 投影主要用于网页上的地理可视化，UTM 是 Universal Transverse Mercator 的缩写。360 度对应全球 60 个时区，每个时区横跨 6 经度。

```{r}
st_transform(
  x = st_sfc(st_point(x = c(2, 49)), crs = 4326),
  crs = st_crs("+proj=utm +zone=32 +ellps=GRS80")
)
```

快速简单绘图，可采用图层 `geom_sf()`，它相当于统计图层 `stat_sf()` 和坐标映射图层 `coord_sf()` 的叠加，`geom_sf()` 支持点、线和多边形等数据数据对象，可以混合叠加。 `coord_sf()` 有几个重要的参数：

1.  `crs`：在绘图前将各个 `geom_sf()` 图层中的**数据**映射到该坐标参考系。

2.  `default_crs`：将非 sf 图层（没有携带 CRS 信息）的数据映射到该坐标参考系，默认使用 `crs` 参数的值，常用设置 `default_crs = sf::st_crs(4326)` 将非 sf 图层中的横纵坐标转化为经纬度，采用 World Geodetic System 1984 (WGS84)。

3.  `datum`：经纬网线的坐标参考系，默认值 `sf::st_crs(4326)`。

下图的右子图将 quakes_sf 数据集投影到坐标参考系统[EPSG:3460](https://epsg.io/3460)。

```{r}
#| label: fig-quakes-ggplot2-grid
#| fig-cap: 斐济地震的空间分布
#| fig-subcap: 
#| - 坐标参考系 4326（默认）
#| - 坐标参考系 3460
#| fig-width: 4
#| fig-height: 4
#| fig-showtext: true
#| layout-ncol: 2

library(ggplot2)
ggplot() +
  geom_sf(data = quakes_sf, aes(color = mag))
ggplot() +
  geom_sf(data = quakes_sf, aes(color = mag)) +
  coord_sf(crs = 3460)
```

数据集 quakes_sf 已经准备了坐标参考系统，此时，`coord_sf()` 就会采用数据集相应的坐标参考系统，即 `sf::st_crs(4326)`。上图的左子图相当于：

```{r}
#| eval: false

ggplot() +
  geom_sf(data = quakes_sf, aes(color = mag)) +
  coord_sf(
    crs = 4326, datum = sf::st_crs(4326),
    default_crs = sf::st_crs(4326)
  )
```

### 凸包操作

```{r}
quakes_sf <- st_transform(quakes_sf, crs = 3460)
# 组合 POINT 构造 POLYGON
quakes_sfp <- st_cast(st_combine(st_geometry(quakes_sf)), "POLYGON")
# 构造 POLYGON 的凸包
quakes_sfp_hull <- st_convex_hull(st_geometry(quakes_sfp))
```

```{r}
#| label: fig-convex-hull
#| fig-cap: 凸包
#| fig-subcap: 
#| - 凸包（base R）
#| - 凸包（ggplot2）
#| fig-showtext: true
#| fig-width: 4
#| fig-height: 4
#| layout-ncol: 2
#| par: true

# 绘制点及其包络
plot(st_geometry(quakes_sf))
# 添加凸包曲线
plot(quakes_sfp_hull, add = TRUE)

ggplot() +
  geom_sf(data = quakes_sf) +
  geom_sf(data = quakes_sfp_hull, fill = NA) +
  coord_sf(crs = 3460, xlim = c(569061, 3008322), ylim = c(1603260, 4665206))
```

## 数据探索

### 核密度估计

给定边界内的[核密度估计与绘制热力图](https://stackoverflow.com/questions/68643517)

```{r}
# spatial point pattern ppp 类型
quakes_ppp <- spatstat.geom::as.ppp(st_geometry(quakes_sf))
# 限制散点在给定的窗口边界内平滑
spatstat.geom::Window(quakes_ppp) <- spatstat.geom::as.owin(quakes_sfp_hull)

# quakes_ppp <- spatstat.geom::as.ppp(st_geometry(quakes_sf), W = spatstat.geom::as.owin(quakes_sfp_hull))

# 密度估计
density_spatstat <- spatstat.explore::density.ppp(quakes_ppp, dimyx = 256)
# 转化为 stars 对象 栅格数据
density_stars <- stars::st_as_stars(density_spatstat)
# 设置坐标参考系
density_sf <- st_set_crs(st_as_sf(density_stars), 3460)
```

### 绘制热力图

```{r}
#| label: fig-kernel-heatmap
#| fig-cap: 热力图
#| fig-subcap: 
#| - 核密度估计
#| - 核密度估计（原始数据）
#| fig-showtext: true
#| fig-width: 4
#| fig-height: 4
#| layout-ncol: 2

ggplot() +
  geom_sf(data = density_sf, aes(fill = v), col = NA) +
  scale_fill_viridis_c() +
  geom_sf(data = st_boundary(quakes_sfp_hull))

ggplot() +
  geom_sf(data = density_sf, aes(fill = v), col = NA) +
  scale_fill_viridis_c() +
  geom_sf(data = st_boundary(quakes_sfp_hull)) +
  geom_sf(data = quakes_sf, size = 1, col = "black")
```