大样地海拔数据的克里格插值与地形图绘制

最近收到几位同学求助， 咨询如何通过全站仪测得的海拔数据绘制大样地的地形图以及计算地形因子，本文作简要说明。本文中，“样方”是大样地中边长为20m的样方。前面的博文 （http://blog.sciencenet.cn/blog-255662-1081363.html ）简要介绍了大样地中样方的命名方法以及物种调查数据的整理过程，可供参考。

大样地研究除了收集物种分布数据，还需要获得每个样方对应的环境因子才能进行进一步分析。环境因子又包括地形变量和土壤因子变量，这里只说地形数据。地形数据一般包括某尺度下（如5m, 10m, 20m, 50m, 100m等）每个样方中心点的海拔、坡度、坡向以及凹凸度等。其中坡度、坡向和凹凸度又都是通过样方中心点海拔矩阵计算得出的。

每个样方中心点的海拔有两种获取方式：（1） 用全站仪直接测量，但是这种情况极少见；（2）对全站仪测得的海拔数据进行克里格插值之后，再在海拔模型上从中心点坐标提取。

为什么不全都用全站仪在每个20m样方的中心点直接测量呢？ 这是由于全站仪在打点时一般都是测量各样方顶点的海拔。有些样方的顶点正好位于深沟内或巨石上，难以获得，此时只能在附近的点测量。因此大多数情况下，全站仪所测得的数据并总不出现在样方顶点。所以第二种处理方式就很常见了。

有些研究中，样方的海拔没有经过克里格插值之后再提取中心点，而是直接用四个顶点海拔的平均值，这也是一种近似的处理方法。仍然是上述原因，相应条件不容易满足。

下面以西双版纳补蚌大样地的海拔数据为例， 介绍海拔数据读取，克里格插值，并提取20m中心点的海拔。获得样方中心点海拔之后，即可计算出坡度、坡向和凹凸度等（http://blog.sciencenet.cn/blog-255662-690600.html ）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
setwd("C:/Users/jlzhang/Desktop/bubeng plot")

library(gstat)

library(sp)

library(raster)

library(openxlsx)

library(lattice)

 

#### 边长为20m米的点阵

bubeng.grid <- expand.grid(x = seq(0, 400, 20), y = seq(0, 500, 20))

coordinates(bubeng.grid) <- ~ x + y

gridded(bubeng.grid) <- TRUE

 

## 数据来源 https://doi.org/10.1371/journal.pone.0108450.s008

bubeng.topo20160318 <- read.xlsx("bubeng_topo.xlsx")

 

## 指定横纵坐标

coordinates(bubeng.topo20160318) = ~ x + y

 

## 变异函数

bubeng_variogram <- variogram(alt ~ 1, data = bubeng.topo20160318)

plot(bubeng_variogram)

 

### 拟合高斯模型， 初始值选择参见 （How Kriging works—Help）

bubeng_variogram.fit = fit.variogram(bubeng_variogram, model = vgm(psill = 10, "Gau", range = 200, nugget = .1))

plot(bubeng_variogram, bubeng_variogram.fit)

 

### 进行克里格插值

bubeng.krig <- krige(alt ~ 1, bubeng.topo20160318, bubeng.grid, model = bubeng_variogram.fit)

 

### 绘图

levelplot(bubeng.krig$var1.pred ~ bubeng.krig$x+bubeng.krig$y, cuts = 50,

          main="Ordinary kriging predictions",contour = TRUE,

          labels = TRUE, col.regions = terrain.colors(255),  pretty = TRUE)

 

## 三维图

wireframe(bubeng.krig$var1.pred ~ bubeng.krig$x+bubeng.krig$y,

          main="Ordinary kriging predictions",

          labels = TRUE, col.regions = terrain.colors(255),  pretty = TRUE, drape = TRUE)

 

## 转换为raster

dat <- data.frame(bubeng.krig$var1.pred, bubeng.krig$x, bubeng.krig$y)

 

### Convert the object x to a raster object.

rdat <- raster(bubeng.krig)

 

### plot the raster

plot(rdat, col = terrain.colors(255))

contour(rdat, add = TRUE, nlevel = 15)

 

### 获得每个20m样方的中心点

x10 <- seq(10, 400-10, by = 20)

y10 <- seq(10, 500-10, by = 20)

grids10 <- expand.grid(x10, y10)

plot(Var2 ~ Var1, data = grids10)

coordinates(grids10) <- c(1, 2)

 

###从raster图层中提取海拔

alt_center <- extract(x = rdat, y = grids10)

 

### 20m样方的中心点海拔数据

res_alt_center <- data.frame(grids10, alt_center)

参考文献

How Kriging works—Help | ArcGIS for Desktop. (2018). Available at: http://desktop.arcgis.com/en/arcmap/10.3/tools/3d-analyst-toolbox/how-kriging-works.htm. Last accessed 24 April 2018.

Hu, Y.-H., Kitching, R.L., Lan, G.-Y., Zhang, J.-L., Sha, L.-Q. & Cao, M. (2014). Size-Class Effect Contributes to Tree Species Assembly through Influencing Dispersal in Tropical Forests. PLoS ONE, 9, e108450.