【翻译】Species distribution modeling 4 环境数据
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")
