【翻译】Species distribution modeling 4 环境数据

4.1 栅格数据

  在物种分布建模中,预测变量通常被组织为栅格(raster)(grid)形式的文件。每个预测变量应该是一个代表感兴趣变量的‘栅格’,变量可以包含气候、土壤、地形、植被、土地使用情况等等。这些数据通常被存储在某种GIS格式的文件中。几乎用到了所有的相关格式(包括ESRI grid,geoTiff,netCDF,IDRISI)。尽量避免使用ASCII文件,因为他们的处理速度相当慢。对于任何特定的研究,图层都应该有相同的空间范围、解析度、来源和投影。必要的话,可以使用raster包中的cropextendaggregateresampleprojectRaster等函数。有关准备预测参数数据的更多信息参见第二章的第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")

  

posted @ 2017-12-28 19:05 钱天宇 阅读(...) 评论(...) 编辑 收藏