504 lines
10 KiB
Plaintext
504 lines
10 KiB
Plaintext
---
|
|
title: "空间分析"
|
|
format:
|
|
dwsd-revealjs:
|
|
logo: _extensions/drwater/dwsd/inst/ucaslogo.png
|
|
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}
|
|
#| results: 'asis'
|
|
rmdify::slideend(
|
|
wechat = FALSE,
|
|
type = "public",
|
|
tel = FALSE,
|
|
thislink = "../"
|
|
)
|
|
```
|