gstat | 空间插值（一）——反距离权重插值；使用ggplot2绘制地图

本篇既是空间插值系列的第一篇推文，也是ggplot2工具包系列推文中的一篇。空间插值使用的工具包是gstat，该工具包主要用于地统计分析。

library(gstat)

示例数据来自HSAR工具包：

library(HSAR)
data("Beijingdistricts")
data("landSPDF")
data("landprice")

• Beijingdistricts：北京五环内乡镇街道层次的矢量文件，sp格式；

• landSPDF：样本点位置信息的矢量文件，sp格式；

• landprice：样本点的土地价格及其他属性信息，数据框格式。

• 以上数据都是真实数据经过随机扰动后得到的。

首先将sp格式的矢量文件转换成sf格式，并进行相应的数据连接，然后使用plot函数绘制样本点土地价格的空间分布图：

library(sf)
library(RColorBrewer)
library(tidyverse)

Beijingdistricts <- st_as_sf(Beijingdistricts)
bjsf <- landSPDF %>% st_as_sf() %>%
  left_join(landprice, by = "obs")
bj.border <- st_boundary(Beijingdistricts)

plot(st_geometry(Beijingdistricts), reset = F)
plot(bjsf["lnprice"], add = T)

1 反距离权重插值

反距离权重插值（IDW）的原理是，已知点属性对插值点属性的权重与两者之间距离的幂成反比。写成公式如下：

其中 为距离衰减参数，取值主要依靠经验。

gstat工具包进行反距离权重插值的函数为idw：

idw(formula, locations, data, newdata, nmax = Inf, nmin = 0, 
    omax = 0, maxdist = Inf, block = numeric(0),
    na.action = na.pass, idp = 2.0, debug.level = 1)

主要的参数有：

• formula：插值表达式，一般为var ~ 1；

• locations：已知点对象；

• data：插值表达式中变量所在的数据框，若locations中的属性已经包含了相关变量，该参数应该省略；

• newdata：插值点对象；

• nmax、nmin：插值任意位置的属性值至多和至少所需的已知点数目；

• maxdist：距离阈值；与插值点超过这个距离的已知点不参与该插值点属性值的计算；

• idp：距离衰减参数，默认取值为2。

创建需要进行属性插值的网格对象：

# 创建网格
bjsf.grid <- st_make_grid(Beijingdistricts, n = c(100, 100))

plot(st_geometry(bjsf.grid), reset = F)
plot(st_geometry(bj.border), add = T)

分别将距离衰减参数设置为1、2、3进行插值。此外，由于空间插值比较费时，为避免重复工作，建议及时将插值出的数据保存到本地。

sf.idw1 <- idw(lnprice ~ 1, bjsf, bjsf.grid, idp = 1)
save(sf.idw1, file = "28-1.Rdata")

sf.idw2 <- idw(lnprice ~ 1, bjsf, bjsf.grid, idp = 2)
save(sf.idw2, file = "28-2.Rdata")

sf.idw3 <- idw(lnprice ~ 1, bjsf, bjsf.grid, idp = 3)
save(sf.idw3, file = "28-3.Rdata")

插值结果如下：

head(sf.idw2)

## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 426987.3 ymin: 4403559 xmax: 429443.4 ymax: 4403956
## CRS:           +proj=tmerc +lat_0=0 +lon_0=117 +k=1 +x_0=500000 +y_0=0 +ellps=krass +units=m +no_defs
##            x       y var1.pred var1.var                       geometry
## ID1 427192.0 4403757  7.337970       NA POLYGON ((426987.3 4403559,...
## ID2 427601.3 4403757  7.336139       NA POLYGON ((427396.7 4403559,...
## ID3 428010.7 4403757  7.334229       NA POLYGON ((427806 4403559, 4...
## ID4 428420.0 4403757  7.332238       NA POLYGON ((428215.3 4403559,...
## ID5 428829.3 4403757  7.330163       NA POLYGON ((428624.7 4403559,...
## ID6 429238.7 4403757  7.328003       NA POLYGON ((429034 4403559, 4...

• 变量var1.pred即插值结果；

• 反距离权重方法无法计算插值结果的方差和置信区间。

2 可视化空间插值结果

这里使用ggplot2工具包对空间插值结果进行可视化。首先从本地加载已经插值出的数据：

load("28-1.Rdata")
load("28-2.Rdata")
load("28-3.Rdata")

提前准备好ggplot2中theme函数的相关设置：

ggplot() +
  theme(text = element_text(family = "mono"),
        axis.line = element_line(color = "black", size = 0.2),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 10),
        axis.ticks = element_line(colour = "grey"),
        
        legend.position = c(1,0),
        legend.justification = c(1,0),
        legend.background = element_blank(),
        legend.title = element_text(face = "bold"),
        
        panel.background = element_blank(),
        panel.border = element_rect(fill = NA, size = 0.3),
        panel.grid.major = element_line(linetype = 2, colour = "grey", size = 0.5),
        
        plot.title = element_text(size = 12, hjust = 0.5,
                                  face = "bold")) -> p

在进行可视化之前，可以先对插值网格进行空间裁剪，使其范围与Beijingdistricts一致。矢量对象的裁剪方法在前面推文中已有介绍：

• sf | 空间矢量对象的几何信息处理方法

sf.idw <- st_interp(sf.idw1, st_union(Beijingdistricts))

ggplot2中绘制矢量地图的函数为geom_sf：

p + geom_sf(data = sf.idw, aes(fill = var1.pred), col = NA) +
  geom_sf(data = bj.border) +
  labs(title = "idp = 1") -> p2
p2

使用scale_fill_continuous函数对填充颜色进行调整。该函数只需指定最低（low）和最高（high）状态下的颜色，然后就会自动计算出连续的颜色条。

p2 + 
  scale_fill_continuous(low = "white", high = "red",
                        space = "rgb", name = "lnprice")

sf.idw <- st_interp(sf.idw2, st_union(Beijingdistricts))

p + geom_sf(data = sf.idw, aes(fill = var1.pred), col = NA) +
  geom_sf(data = bj.border) +
  scale_fill_continuous(low = "white", high = "red",
                        space = "rgb", name = "lnprice")+
  labs(title = "idp = 2")

sf.idw <- st_interp(sf.idw3, st_union(Beijingdistricts))

p + geom_sf(data = sf.idw, aes(fill = var1.pred), col = NA) +
  geom_sf(data = bj.border) +
  scale_fill_continuous(low = "white", high = "red",
                        space = "rgb", name = "lnprice")+
  labs(title = "idp = 3")

此外，还可以将插值的结果重新聚合到乡镇街道层次上，这里使用按面积加权的方法进行聚合操作。

相应的函数方法见推文sf | 空间矢量对象的“聚合”操作。

street <- st_interpolate_aw(sf.idw2, Beijingdistricts, extensive = F)

p + geom_sf(data = street, aes(fill = var1.pred), col = NA) +
  geom_sf(data = bj.border) +
  scale_fill_continuous(low = "white", high = "red",
                        space = "rgb", name = "lnprice")

---