张金龙的博客分享 http://blog.sciencenet.cn/u/zjlcas 物种适应性、分布与进化

博文

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

已有 5679 次阅读 2018-4-24 20:11 |系统分类:科研笔记|关键词:学者

 

最近收到几位同学求助, 咨询如何通过全站仪测得的海拔数据绘制大样地的地形图以及计算地形因子,本文作简要说明。本文中,“样方”是大样地中边长为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 )


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.




https://m.sciencenet.cn/blog-255662-1110708.html

上一篇:网络时光与网络时光机
下一篇:用R绘制采样点地图——带指北针、比例尺、图例

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2022-11-29 15:42

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部