4.1 栅格数据
在物种分布建模中,预测变量通常被组织为栅格(raster)(grid)形式的文件。每个预测变量应该是一个代表感兴趣变量的‘栅格’,变量可以包含气候、土壤、地形、植被、土地使用情况等等。这些数据通常被存储在某种GIS格式的文件中。几乎用到了所有的相关格式(包括ESRI grid,geoTiff,netCDF,IDRISI)。尽量避免使用ASCII文件,因为他们的处理速度相当慢。对于任何特定的研究,图层都应该有相同的空间范围、解析度、来源和投影。必要的话,可以使用raster包中的crop、extend、aggregate、resample和projectRaster等函数。有关准备预测参数数据的更多信息参见第二章的第8节。查看raster包的帮助文档和简介以获得有关这个的更多信息。
预测参数(rasters)的设置可以被用于制作RasterStack,一个‘RasterLayer’对象的集合(更多信息参见第二章第4节)。
我们制作了一个随dismo包安装的文件的名单然后从中创建了一个rasterStack,展示出每个图层的名称,最后将它们全部绘制出来。
首先得到我们所需文件的文件夹。在本例中是dismo包的“ex”目录。在你自己的文件上你无需操作此复杂的步骤(输入字符串就可以)。
path <- file.path(system.file(package="dismo"), 'ex')
现在获取这个文件夹中所有扩展名为“.grd”的文件名称。符号$声明了文件必须以字符‘grd’结尾。使用full.names=TRUE,返回完整的路径名。
library(dismo)
files <- list.files(path, pattern='grd$', full.names=TRUE )
files
## [1] "C:/soft/R/R-3.3.1/library/dismo/ex/bio1.grd"
## [2] "C:/soft/R/R-3.3.1/library/dismo/ex/bio12.grd"
## [3] "C:/soft/R/R-3.3.1/library/dismo/ex/bio16.grd"
## [4] "C:/soft/R/R-3.3.1/library/dismo/ex/bio17.grd"
## [5] "C:/soft/R/R-3.3.1/library/dismo/ex/bio5.grd"
## [6] "C:/soft/R/R-3.3.1/library/dismo/ex/bio6.grd"
## [7] "C:/soft/R/R-3.3.1/library/dismo/ex/bio7.grd"
## [8] "C:/soft/R/R-3.3.1/library/dismo/ex/bio8.grd"
## [9] "C:/soft/R/R-3.3.1/library/dismo/ex/biome.grd"
现在创建一个预测参数的RasterStack。
predictors <- stack(files)
predictors
## class : RasterStack
## dimensions : 192, 186, 35712, 9 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -125, -32, -56, 40 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## names : bio1, bio12, bio16, bio17, bio5, bio6, bio7, bio8, biome
## min values : -23, 0, 0, 0, 61, -212, 60, -66, 1
## max values : 289, 7682, 2458, 1496, 422, 242, 461, 323, 14
names(predictors)
## [1] "bio1" "bio12" "bio16" "bio17" "bio5" "bio6" "bio7" "bio8" "biome"
plot(predictors)
我们也可以在RasterStack中绘制一个单独的图层,然后在其上绘制一些附加数据。
library(maptools)
## Checking rgeos availability: TRUE
data(wrld_simpl)
file <- paste(system.file(package="dismo"), "/ex/bradypus.csv", sep="")
bradypus <- read.table(file, header=TRUE, sep=',')
# we do not need the first column
bradypus <- bradypus[,-1]
然后现在绘制:
# first layer of the RasterStack
plot(predictors, 1)
# note the "add=TRUE" argument with plot
plot(wrld_simpl, add=TRUE)
# with the points function, "add" is implicit
points(bradypus, col='blue')
上方的例子使用了WorldClim database中的数据代表‘bioclimatic variables’(Hijmans et al., 2004)和WWF中的‘terrestrial biome’数据(Olsen et al., 2001)。如果你想要更高分辨率的数据可以去网站中获得。你同样可以使用raster包中的getData函数下载WorldClim气候数据。
预测参数的选择可能很重要,尤其在研究的目的是解释的情况下。参加例子Austin and Smith(1987),Austin(2002),Mellert et al., (2001)。早期的物种建模应用倾向于关注解释(Elith and Leathwick 2009)。现今,SDM的目的倾向于预测。在相同地理区域内的预测,参数选择可能相对不那么重要,但对于多数预测任务来说(例如对于新的时间,新的地点,如下)参数选择是至关重要的。在所有的情况下使用与物种生态位相对应的参数(而不是你能在网络上找到的第一个数据集)是重要的。有些情况下从现有的数据中开发新的、与生态更加相关的预测参数是很有用的。举个例子,可以使用土地覆盖数据和raster包中的focal函数去建立一个新的参数去表示在网格单元的x km内有多少森林面积是可用的,对于一个物种来说x可能就是家庭范围。
4.2 从栅格中提取值
我们现在拥有一组预测变量(rasters)和出现点。下一步就是从点的位置提取预测参数的值。(对于已经在dismo包中实现的建模方法,这一步可以省略)。使用raster包中的‘extract'方法让这件事变得十分简单。在下方的例子中我们首先使用这个函数来提取Bradypus的发生点,然后提取500个随机的背景点。我们把这些结合到一个单独的data.frame里,其第一列(变量‘pb’)标明这是一个存在点还是一个背景点。明确地定义‘biome’是一个分类变量(在R中叫做因子)是重要的,因此它不会被看作为任何其他数字变量来对待。
presvals <- extract(predictors, bradypus)
# setting random seed to always create the same
# random set of points for this example
set.seed(0)
backgr <- randomPoints(predictors, 500)
absvals <- extract(predictors, backgr)
pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals)))
sdmdata <- data.frame(cbind(pb, rbind(presvals, absvals)))
sdmdata[,'biome'] = as.factor(sdmdata[,'biome'])
head(sdmdata)
## pb bio1 bio12 bio16 bio17 bio5 bio6 bio7 bio8 biome
## 1 1 263 1639 724 62 338 191 147 261 1
## 2 1 263 1639 724 62 338 191 147 261 1
## 3 1 253 3624 1547 373 329 150 179 271 1
## 4 1 243 1693 775 186 318 150 168 264 1
## 5 1 243 1693 775 186 318 150 168 264 1
## 6 1 252 2501 1081 280 326 154 172 270 1
tail(sdmdata)
## pb bio1 bio12 bio16 bio17 bio5 bio6 bio7 bio8 biome
## 611 0 253 2341 928 166 318 178 141 256 1
## 612 0 251 1191 528 74 329 160 169 271 9
## 613 0 175 161 67 7 335 35 300 179 13
## 614 0 263 2789 943 418 318 211 107 263 1
## 615 0 269 1859 840 32 351 184 168 269 1
## 616 0 117 1136 341 243 297 -67 365 201 4
summary(sdmdata)
## pb bio1 bio12 bio16
## Min. :0.0000 Min. : -3.0 Min. : 3.0 Min. : 2.0
## 1st Qu.:0.0000 1st Qu.:177.0 1st Qu.: 752.2 1st Qu.: 327.5
## Median :0.0000 Median :242.0 Median :1489.5 Median : 655.5
## Mean :0.1883 Mean :214.3 Mean :1574.5 Mean : 643.2
## 3rd Qu.:0.0000 3rd Qu.:261.0 3rd Qu.:2250.8 3rd Qu.: 911.8
## Max. :1.0000 Max. :286.0 Max. :7682.0 Max. :2458.0
##
## bio17 bio5 bio6 bio7
## Min. : 0.0 Min. : 67.0 Min. :-207.00 Min. : 68.0
## 1st Qu.: 37.0 1st Qu.:302.8 1st Qu.: 42.75 1st Qu.:117.0
## Median : 98.5 Median :321.0 Median : 160.00 Median :159.5
## Mean : 159.2 Mean :309.8 Mean : 119.32 Mean :190.4
## 3rd Qu.: 230.5 3rd Qu.:333.2 3rd Qu.: 201.25 3rd Qu.:251.5
## Max. :1496.0 Max. :405.0 Max. : 232.00 Max. :443.0
##
## bio8 biome
## Min. :-33.0 1 :280
## 1st Qu.:218.8 13 : 79
## Median :250.5 7 : 63
## Mean :225.0 8 : 52
## 3rd Qu.:262.0 2 : 46
## Max. :316.0 (Other): 95
## NA's : 1
在这里可能有其他的方法。举个例子,作为一个潜在的手段,在处理位置精度和网格单元尺寸之间的不匹配中,可以提取某个半径中的多个点。如果在这个半径中可以构造代表10个同样有效的环境“样本(samples)”的10个数据集,那就可以用来匹配10个模型来探索位置不确定性的影响。
为了在环境数据中将调查共线性(colinearity)可视化(visually)(在存在和背景点上)你可以使用对绘制(pairs plot)。参见Dormann et al.(2013)中对于移除共线性的方法的讨论。
pairs(sdmdata[,2:5], cex=0.1)
在Bradypus出现地点的一个关于气候数据的值的对绘制(pairs plot)。
为了在下一章使用sdmdata和presvals函数,我们将其存盘。
saveRDS(sdmdata, "sdm.Rds")
saveRDS(presvals, "pvals.Rds")