518 lines
11 KiB
Plaintext
518 lines
11 KiB
Plaintext
---
|
||
title: "空间分析"
|
||
subtitle: 《区域水环境污染数据分析实践》<br>Data analysis practice of regional water environment pollution
|
||
author: 苏命、王为东<br>中国科学院大学资源与环境学院<br>中国科学院生态环境研究中心
|
||
date: today
|
||
lang: zh
|
||
format:
|
||
revealjs:
|
||
theme: dark
|
||
slide-number: true
|
||
chalkboard:
|
||
buttons: true
|
||
preview-links: auto
|
||
lang: zh
|
||
toc: true
|
||
toc-depth: 1
|
||
toc-title: 大纲
|
||
logo: ./_extensions/inst/img/ucaslogo.png
|
||
css: ./_extensions/inst/css/revealjs.css
|
||
pointer:
|
||
key: "p"
|
||
color: "#32cd32"
|
||
pointerSize: 18
|
||
revealjs-plugins:
|
||
- pointer
|
||
filters:
|
||
- d2
|
||
knitr:
|
||
opts_chunk:
|
||
dev: "svg"
|
||
retina: 3
|
||
execute:
|
||
freeze: auto
|
||
cache: true
|
||
echo: true
|
||
fig-width: 5
|
||
fig-height: 6
|
||
---
|
||
|
||
|
||
## 目录
|
||
1. 空间数据基础概念
|
||
2. Simple Features (SF) 核心规范
|
||
3. SFEXT扩展功能实战
|
||
4. 空间关系与几何操作
|
||
5. 空间统计分析方法
|
||
6. 综合实战案例
|
||
|
||
|
||
|
||
## 1. 空间数据基础概念
|
||
|
||
### 空间数据类型体系
|
||
```r
|
||
library(sf)
|
||
# 创建示例几何类型
|
||
st_point(c(0, 0)) # 点
|
||
st_linestring(rbind(c(0, 0), c(1, 1))) # 线
|
||
st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0)))) # 面
|
||
```
|
||
|
||
### 坐标参考系(CRS)
|
||
```r
|
||
# 常用坐标系定义
|
||
wgs84 <- st_crs(4326) # 经纬度坐标
|
||
utm50n <- st_crs(32650) # UTM 50N投影
|
||
cat(st_crs(wgs84)$proj4string)
|
||
```
|
||
|
||
|
||
|
||
## 2. Simple Features (SF) 核心规范
|
||
|
||
### 创建SF对象
|
||
```r
|
||
# 从数据框转换
|
||
df <- data.frame(
|
||
city = c("Beijing", "Shanghai"),
|
||
pop = c(2171, 2424),
|
||
geometry = st_sfc(
|
||
st_point(c(116.4, 39.9)),
|
||
st_point(c(121.5, 31.2))
|
||
)
|
||
)
|
||
sf_df <- st_as_sf(df, crs = 4326)
|
||
```
|
||
|
||
### 数据IO操作
|
||
```r
|
||
# 写入/读取空间文件
|
||
st_write(sf_df, "cities.shp")
|
||
new_sf <- st_read("cities.shp")
|
||
|
||
# 从内置数据集加载
|
||
nc <- st_read(system.file("shape/nc.shp", package = "sf"))
|
||
```
|
||
|
||
|
||
|
||
## 3. SFEXT扩展功能实战
|
||
|
||
### 高级几何生成
|
||
```r
|
||
remotes::install_github("r-spatial/sfext")
|
||
library(sfext)
|
||
|
||
# 生成渔网网格
|
||
grid <- st_regular_grid(
|
||
xmin = 0,
|
||
ymin = 0,
|
||
xmax = 100,
|
||
ymax = 100,
|
||
cell_size = 10
|
||
)
|
||
plot(grid)
|
||
```
|
||
|
||
### 空间索引加速
|
||
```r
|
||
# 创建空间索引
|
||
system.time(
|
||
without_idx <- st_intersects(nc, nc)
|
||
)
|
||
|
||
system.time(
|
||
with_idx <- st_intersects(st_make_valid(nc), sparse = TRUE)
|
||
)
|
||
```
|
||
|
||
|
||
|
||
## 4. 空间关系与操作
|
||
|
||
### 拓扑关系判断
|
||
```r
|
||
p1 <- st_point(c(0, 0))
|
||
p2 <- st_point(c(1, 1))
|
||
st_distance(p1, p2) # 距离计算
|
||
|
||
poly <- st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))
|
||
st_contains(poly, p1) # 包含关系
|
||
```
|
||
|
||
### 几何变换操作
|
||
```r
|
||
# 缓冲区分析
|
||
buf <- st_buffer(nc[1, ], dist = 0.1)
|
||
plot(buf)
|
||
plot(nc[1, ]$geometry, add = TRUE, col = 'red')
|
||
|
||
# 凸包计算
|
||
convex <- st_convex_hull(nc)
|
||
plot(convex)
|
||
```
|
||
|
||
|
||
|
||
## 5. 空间统计分析方法
|
||
|
||
### 点模式分析
|
||
```r
|
||
library(spatstat)
|
||
# 创建点模式对象
|
||
ppp <- ppp(x = runif(100), y = runif(100), window = owin(c(0, 1), c(0, 1)))
|
||
|
||
# K函数分析
|
||
K <- Kest(ppp)
|
||
plot(K, main = "Ripley's K函数")
|
||
```
|
||
|
||
### 空间自相关检验
|
||
```r
|
||
library(spdep)
|
||
# Moran's I检验
|
||
nb <- poly2nb(nc)
|
||
lw <- nb2listw(nb)
|
||
moran.test(nc$BIR74, lw)
|
||
```
|
||
|
||
|
||
|
||
## 6. 综合实战案例
|
||
|
||
### 城市公园服务区分析
|
||
```r
|
||
library(osmdata)
|
||
library(sfnetworks)
|
||
|
||
# 获取OSM数据
|
||
parks <- opq("New York") %>%
|
||
add_osm_feature("leisure", "park") %>%
|
||
osmdata_sf()
|
||
|
||
roads <- opq("New York") %>%
|
||
add_osm_feature("highway") %>%
|
||
osmdata_sf()
|
||
|
||
# 构建网络分析
|
||
net <- as_sfnetwork(roads$osm_lines) %>%
|
||
st_transform(32618)
|
||
|
||
# 计算服务范围（示例）
|
||
service_area <- st_network_blend(net, parks$osm_polygons) %>%
|
||
st_buffer(500) # 500米服务区
|
||
```
|
||
|
||
|
||
|
||
## 可视化增强技巧
|
||
|
||
### 分层绘图
|
||
|
||
```r
|
||
library(ggplot2)
|
||
ggplot() +
|
||
geom_sf(data = nc, aes(fill = BIR74)) +
|
||
geom_sf(data = service_area, alpha = 0.3) +
|
||
scale_fill_viridis_c() +
|
||
theme_minimal()
|
||
```
|
||
|
||
### 交互可视化
|
||
|
||
```r
|
||
library(mapview)
|
||
mapview(nc, zcol = "BIR74", burst = TRUE) +
|
||
mapview(service_area, alpha.regions = 0.2)
|
||
```
|
||
|
||
# 实践案例
|
||
|
||
|
||
## 中国地图
|
||
|
||
::: panel-tabset
|
||
### Code
|
||
|
||
|
||
```{r}
|
||
#| echo: true
|
||
#| eval: false
|
||
#| out-width: 50%
|
||
library(tidyverse)
|
||
library(sf)
|
||
|
||
# 更稳定的预处理数据源（含港澳台）
|
||
# china <- st_read("https://geo.datav.aliyun.com/areas_v3/bound/100000.json")
|
||
|
||
china <- st_read("./中华人民共和国.json", quiet = TRUE)
|
||
|
||
library(ggplot2)
|
||
|
||
ggplot(china) +
|
||
geom_sf(fill = "#F6F6F6", color = "#666666", size = 0.3) +
|
||
theme_void() +
|
||
labs(title = "China")
|
||
|
||
# 添加模拟经济数据（实际应用可连接统计年鉴）
|
||
set.seed(123)
|
||
china$gdp <- runif(nrow(china), 100, 1000) # 模拟GDP数据（亿元）
|
||
|
||
china |>
|
||
ggplot() +
|
||
geom_sf(aes(fill = gdp), color = "gray90", size = 0.2) +
|
||
scale_fill_gradientn(
|
||
colours = c("#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c"),
|
||
name = NULL
|
||
) +
|
||
theme_void() +
|
||
theme(legend.position = c(0.85, 0.3))
|
||
```
|
||
|
||
### Output
|
||
|
||
```{r}
|
||
#| echo: false
|
||
#| out-width: 80%
|
||
#| fig-width: 6
|
||
#| fig-height: 3
|
||
library(tidyverse)
|
||
library(sf)
|
||
|
||
# 更稳定的预处理数据源（含港澳台）
|
||
# china <- st_read("https://geo.datav.aliyun.com/areas_v3/bound/100000.json")
|
||
|
||
china <- st_read("./中华人民共和国.json", quiet = TRUE)
|
||
|
||
library(ggplot2)
|
||
|
||
ggplot(china) +
|
||
geom_sf(fill = "#F6F6F6", color = "#666666", size = 0.3) +
|
||
theme_void() +
|
||
labs(title = "China")
|
||
|
||
# 添加模拟经济数据（实际应用可连接统计年鉴）
|
||
set.seed(123)
|
||
china$gdp <- runif(nrow(china), 100, 1000) # 模拟GDP数据（亿元）
|
||
|
||
china |>
|
||
ggplot() +
|
||
geom_sf(aes(fill = gdp), color = "gray90", size = 0.2) +
|
||
scale_fill_gradientn(
|
||
colours = c("#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c"),
|
||
name = NULL
|
||
) +
|
||
theme_void() +
|
||
theme(legend.position = c(0.85, 0.3))
|
||
```
|
||
|
||
:::
|
||
|
||
|
||
## Kriging
|
||
|
||
::: panel-tabset
|
||
### Code
|
||
|
||
```{r}
|
||
#| echo: true
|
||
#| eval: false
|
||
#| out-width: 50%
|
||
require(sf)
|
||
require(sp)
|
||
require(gstat)
|
||
require(raster)
|
||
require(ggplot2)
|
||
# 读取太湖边界
|
||
lake_boundary <- st_read("../../data/taihu.shp") |>
|
||
sf::st_make_valid()
|
||
main_water <- st_difference(
|
||
lake_boundary[lake_boundary$Name == "boundary", ],
|
||
st_union(lake_boundary[lake_boundary$Name != "boundary", ])
|
||
)
|
||
|
||
get_kriging <- function(indf, Yvarname, grid = NULL, boundsf = main_water) {
|
||
insf <- indf |>
|
||
sfext::df_to_sf(crs = 4326, coords = c("long", "lat")) |>
|
||
dplyr::select(-c("long", "lat")) |>
|
||
dplyr::rename(Y = tidyselect::all_of(Yvarname))
|
||
if (is.null(grid)) {
|
||
# 1. 将sf边界转为terra格式
|
||
v <- terra::vect(boundsf)
|
||
# 2. 创建基础网格
|
||
base_grid <- terra::rast(v, resolution = 0.002) # 基础低分辨率
|
||
# 3. 创建高分辨率区域（例如边界附近）
|
||
buffer_zone <- buffer(v, width = 0.005) # 边界附近创建缓冲区
|
||
hi_res_grid <- terra::rast(buffer_zone, resolution = 0.001)
|
||
# 4. 合并网格
|
||
final_grid <- merge(base_grid, hi_res_grid)
|
||
# 5. 转为点并裁剪
|
||
grid <- terra::as.points(final_grid) |>
|
||
sf::st_as_sf() |>
|
||
sf::st_filter(main_water)
|
||
}
|
||
insp <- sf::as_Spatial(insf)
|
||
fit <- automap::autofitVariogram(Y ~ 1, insp)
|
||
# 克里金插值
|
||
m <- gstat::gstat(
|
||
formula = Y ~ 1,
|
||
data = insp,
|
||
model = fit$var_model
|
||
)
|
||
predsf <- predict(m, grid) |>
|
||
sf::st_as_sf()
|
||
outsf <- predsf |>
|
||
st_coordinates() |>
|
||
as.data.frame() |>
|
||
cbind(pred = predsf$var1.pred)
|
||
return(outsf)
|
||
}
|
||
|
||
|
||
wqdf <- readxl::read_xlsx("../../data/wqdata.xlsx")
|
||
|
||
if (FALSE) {
|
||
wqsf <- wqdf |>
|
||
nest(datedf = -date) |>
|
||
dplyr::mutate(
|
||
krigingsf = purrr::map(datedf, \(x) {
|
||
get_kriging(x, "conc", boundsf = main_water)
|
||
})
|
||
)
|
||
saveRDS(wqsf, "./wqsf.rds")
|
||
}
|
||
wqsf <- readRDS("./wqsf.rds")
|
||
|
||
wqsf |>
|
||
unnest(krigingsf) |>
|
||
ggplot(aes(X, Y)) +
|
||
geom_contour_filled(aes(z = pred), bins = 20) +
|
||
geom_sf(
|
||
data = main_water,
|
||
aes(x = NULL, y = NULL),
|
||
fill = NA,
|
||
colour = "black",
|
||
linewidth = 1
|
||
) +
|
||
# scale_fill_gradientn(colors = hcl.colors(100, "RdYlBu")) + # 设置颜色渐变
|
||
# ggsci::scale_fill_aaas() +
|
||
coord_sf() +
|
||
theme_void()
|
||
```
|
||
|
||
### Output
|
||
|
||
```{r}
|
||
#| echo: false
|
||
#| out-width: 80%
|
||
#| message: false
|
||
#| fig-width: 6
|
||
#| fig-height: 3
|
||
require(sf)
|
||
require(sp)
|
||
require(gstat)
|
||
require(raster)
|
||
require(ggplot2)
|
||
# 读取太湖边界
|
||
lake_boundary <- st_read("../../data/taihu.shp", quiet = TRUE) |>
|
||
sf::st_make_valid()
|
||
main_water <- st_difference(
|
||
lake_boundary[lake_boundary$Name == "boundary", ],
|
||
st_union(lake_boundary[lake_boundary$Name != "boundary", ])
|
||
)
|
||
|
||
get_kriging <- function(indf, Yvarname, grid = NULL, boundsf = main_water) {
|
||
insf <- indf |>
|
||
sfext::df_to_sf(crs = 4326, coords = c("long", "lat")) |>
|
||
dplyr::select(-c("long", "lat")) |>
|
||
dplyr::rename(Y = tidyselect::all_of(Yvarname))
|
||
if (is.null(grid)) {
|
||
# 1. 将sf边界转为terra格式
|
||
v <- terra::vect(boundsf)
|
||
# 2. 创建基础网格
|
||
base_grid <- terra::rast(v, resolution = 0.002) # 基础低分辨率
|
||
# 3. 创建高分辨率区域（例如边界附近）
|
||
buffer_zone <- buffer(v, width = 0.005) # 边界附近创建缓冲区
|
||
hi_res_grid <- terra::rast(buffer_zone, resolution = 0.001)
|
||
# 4. 合并网格
|
||
final_grid <- merge(base_grid, hi_res_grid)
|
||
# 5. 转为点并裁剪
|
||
grid <- terra::as.points(final_grid) |>
|
||
sf::st_as_sf() |>
|
||
sf::st_filter(main_water)
|
||
}
|
||
insp <- sf::as_Spatial(insf)
|
||
fit <- automap::autofitVariogram(Y ~ 1, insp)
|
||
# 克里金插值
|
||
m <- gstat::gstat(
|
||
formula = Y ~ 1,
|
||
data = insp,
|
||
model = fit$var_model
|
||
)
|
||
predsf <- predict(m, grid) |>
|
||
sf::st_as_sf()
|
||
outsf <- predsf |>
|
||
st_coordinates() |>
|
||
as.data.frame() |>
|
||
cbind(pred = predsf$var1.pred)
|
||
return(outsf)
|
||
}
|
||
|
||
|
||
wqdf <- readxl::read_xlsx("../../data/wqdata.xlsx")
|
||
|
||
if (FALSE) {
|
||
wqsf <- wqdf |>
|
||
nest(datedf = -date) |>
|
||
dplyr::mutate(
|
||
krigingsf = purrr::map(datedf, \(x) {
|
||
get_kriging(x, "conc", boundsf = main_water)
|
||
})
|
||
)
|
||
saveRDS(wqsf, "./wqsf.rds")
|
||
}
|
||
wqsf <- readRDS("./wqsf.rds")
|
||
|
||
wqsf |>
|
||
unnest(krigingsf) |>
|
||
ggplot(aes(X, Y)) +
|
||
geom_contour_filled(aes(z = pred), bins = 20) +
|
||
geom_sf(
|
||
data = main_water,
|
||
aes(x = NULL, y = NULL),
|
||
fill = NA,
|
||
colour = "black",
|
||
linewidth = 1
|
||
) +
|
||
# scale_fill_gradientn(colors = hcl.colors(100, "RdYlBu")) + # 设置颜色渐变
|
||
# ggsci::scale_fill_aaas() +
|
||
coord_sf() +
|
||
theme_void()
|
||
```
|
||
|
||
:::
|
||
|
||
## 学习资源
|
||
|
||
1. SF官方文档：https://r-spatial.github.io/sf/
|
||
2. 《Geocomputation with R》在线书
|
||
3. R-spatial教程：https://github.com/r-spatial/workshops
|
||
```
|
||
|
||
### 使用说明：
|
||
1. 需要R 4.0+环境
|
||
2. 推荐安装依赖：
|
||
```r
|
||
install.packages(c("sf", "spatstat", "spdep", "osmdata", "sfnetworks", "ggplot2", "mapview"))
|
||
remotes::install_github("r-spatial/sfext")
|
||
```
|
||
3. 实际案例数据可替换为本地shapefile
|
||
4. 交互可视化部分需要浏览器支持
|
||
|
||
|
||
## 欢迎讨论！{.center}
|
||
|
||
|
||
`r rmdify::slideend(wechat = FALSE, type = "public", tel = FALSE, thislink = "https://drc.drwater.net/course/public/RWEP/PUB/SD/")`
|