【翻译】Species distribution modeling 2 数据准备
2.准备数据
2.1 物种出现数据
向R中导入出现数据很简单。但收集、地理引用和交叉验证坐标数据是乏味的。对物种分布建模的讨论通常关注点在于对比建模方法，但当你要对付一些记录少或不确定的物种时，你的关注点恐怕应该放在提高出现数据的质量了。如果你的出现数据不偏不倚、没有错误并且你的记录数量相对较多，那么所有的方法都会做得更好。我们将为你演示你在R中会用到的一些有用的数据准备步骤，用好额外的工具是必要的。比如QGIS，一个在点数据集的交互式编辑中很有用的程序。
2.2 导入出现数据
大多数情况下你会拥有一个代表某个物种的已知分布的点数据文件。下面是一个使用read.table去读取存储在一个文本文档中的记录。#后为注释。
我们正在使用随dismo包安装的一个示例文件，因此我们构建这个文件名的方式稍显复杂，但你可以用你自己的文件名去取代它。（记得在文件名路径前使用正斜杠"/"！）。system.file将文件路径引入到dismo包安装的路径。如果你之前没有使用过paste函数，你应该去熟悉它（在命令行窗口输入?paste）。
# loads the dismo library library(dismo) ## Loading required package: raster ## Loading required package: sp file <- paste(system.file(package="dismo"), "/ex/bradypus.csv", sep="") # this is the file we will use: file ## [1] "C:/soft/R/R-devel/library/dismo/ex/bradypus.csv"
现在读取这个文件并检查其中的数据。
bradypus <- read.table(file, header=TRUE, sep=",") # or do: bradypus <- read.csv(file) # first rows head(bradypus) ## species lon lat ## 1 Bradypus variegatus -65.4000 -10.3833 ## 2 Bradypus variegatus -65.3833 -10.3833 ## 3 Bradypus variegatus -65.1333 -16.8000 ## 4 Bradypus variegatus -63.6667 -17.4500 ## 5 Bradypus variegatus -63.8500 -17.4000 ## 6 Bradypus variegatus -64.4167 -16.0000 # we only need columns 2 and 3: bradypus <- bradypus[,2:3] head(bradypus) ## lon lat ## 1 -65.4000 -10.3833 ## 2 -65.3833 -10.3833 ## 3 -65.1333 -16.8000 ## 4 -63.6667 -17.4500 ## 5 -63.8500 -17.4000 ## 6 -64.4167 -16.0000
你也可以从Excel或一个数据库中正确地读取这种数据（见RODBC包中的例子）。因为这是一个csv（comma speparated values）文件，我们还应该用到的函数是read.csv。不管你用什么方法，目的都是得到一个包括至少两列的矩阵（或一个data.frame）以包含物种发现位点的坐标。坐标通常表达为经纬度（i.e. angular)，但在UTM投影中也可以表达为Easting and Northing。或是其他平面坐标系（地图投影）。惯例是将这些坐标组织为经度在第一列，纬度在第二列（设想图中的x,y轴，经度是x，纬度是y）；它们经常被列为相反的数据，导致所不期待的结果。在很多情况下你需要额外的列，例如，当你在建模许多物种时需要一个标明物种的列；和一个标明这条记录是“presence”还是“absence”的列（一个通用的惯例是将presence表示为1而absence表示为0）。
如果你不持有任何物种分布数据，你可以在GBIF（Global Biodiversity Inventory Facility)中下载。在dismo包中有一个函数gbif可以用于此。下方的数据就是使用gbif函数下载的（而且保存为一个常设数据集以便使用）：
acaule <- gbif("solanum", "acaule*", geo=FALSE)
如果你想去理清gbif函数所包含这些参数的顺序或是找出其他你所能用到的参数，请查阅帮助文档（记住，帮助文档在library没有被加载之前是无法访问的），输入：?gbif或者是help(gbif)。注意“acaule*”中的星号的作用，不仅请求了Solanum acaule，也包含了类似全名的多种变化，比如Solanum acaule* Bitter,或者是亚种*Solanum acaule subsp.aemulans.
很多发现记录可能是没有地理坐标的。在本例中，超过1366条记录被GBIF所返回（2013.1），这其中有1082条记录包含了地理坐标（这里有699和54条是2010年3月提供的，一个巨大的进步！）
# load the saved S. acaule data data(acaule) # how many rows and colums? dim(acaule) ## [1] 1366 25 #select the records that have longitude and latitude data colnames(acaule) ## [1] "species" "continent" ## [3] "country" "adm1" ## [5] "adm2" "locality" ## [7] "lat" "lon" ## [9] "coordUncertaintyM" "alt" ## [11] "institution" "collection" ## [13] "catalogNumber" "basisOfRecord" ## [15] "collector" "earliestDateCollected" ## [17] "latestDateCollected" "gbifNotes" ## [19] "downloadDate" "maxElevationM" ## [21] "minElevationM" "maxDepthM" ## [23] "minDepthM" "ISO2" ## [25] "cloc" acgeo <- subset(acaule, !is.na(lon) & !is.na(lat)) dim(acgeo) ## [1] 1082 25 # show some values acgeo[1:4, c(1:5,7:10)] ## species continent country adm1 adm2 ## 1 Solanum acaule Bitter South America Argentina Jujuy Santa Catalina ## 2 Solanum acaule Bitter South America Peru Cusco Canchis ## 3 Solanum acaule f. acaule <NA> Argentina <NA> <NA> ## 4 Solanum acaule f. acaule <NA> Bolivia <NA> <NA> ## lat lon coordUncertaintyM alt ## 1 -21.9000 -66.1000 <NA> NaN ## 2 -13.5000 -71.0000 <NA> 4500 ## 3 -22.2666 -65.1333 <NA> 3800 ## 4 -18.6333 -66.9500 <NA> 3700
下方是以Solanum acaule为例的一个简单的制作发现位点地图的方法。制作这样的地图以保证这些点至少粗略的位于正确的位置很重要。
library(maptools) ## Checking rgeos availability: TRUE data(wrld_simpl) plot(wrld_simpl, xlim=c(-80,70), ylim=c(-60,10), axes=TRUE, col="light yellow") # restore the box around the map box() # add the points points(acgeo$lon, acgeo$lat, col='orange', pch=20, cex=0.75) # plot points again to add a border, for better visibility points(acgeo$lon, acgeo$lat, col='red', cex=0.75)
数据集wrld_simpl包含粗略的国界线。你可以使用其他的图形（或线或点的）数据集。举个例子，你可以使用raster包中的getData函数下载更高分辨率的国家和地方行政边界数据。你也可以使用raster包中的shapefile函数向R中读取你自己的shapefile数据。
2.3 数据清理
在处理来源于像是GBIF这种物种分布数据仓库的数据时，数据“清理”显得尤为重要。有些工作不是特别的为了物种分布建模而收集的数据，所以为了你的项目，你需要去理解这些数据并适当的清理它们。这里我们举一个例子。
Solanum acaule是一个在秘鲁南部、玻利维亚和阿根廷北部广泛出现的物种。你看到地图中的哪些错误了吗？
这里有几条标记在巴基斯坦南部海洋中的记录。想一想为什么会发生？这是个常见的错误，忽略了减号。这些坐标在（65.4，23.4）周围，但它们应该在阿根廷北部（-65.4，-23.4）周围（你可以使用“click”函数去查询地图中的坐标）。有两条记录（303行和886行）在地图上标记了位于南极洲（-76.3，-76.3）的同样的点。其地点描述表示它应该位于秘鲁首都利马附近的Huarochiri。所以这里的经度应该是正确的，而纬度的拷贝是错误的。有趣的是这条记录出现了两次。最初的来源是国际马铃薯中心，而另一个由“SINGER”提供的拷贝似乎已经将国家一栏“改正”为南极洲：
acaule[c(303,885),1:10] ## species continent country adm1 adm2 ## 303 solanum acaule acaule <NA> Antarctica <NA> <NA> ## 885 solanum acaule acaule BITTER <NA> Peru <NA> <NA> ## locality lat lon coordUncertaintyM alt ## 303 <NA> -76.3 -76.3 <NA> NaN ## 885 Lima P. Huarochiri Pacomanta -76.3 -76.3 <NA> 3800
一个巴西的点（记录在acaule[98,])本应在玻利维亚南部，所以这里可能有一个经度上的错误。同样地，有三条包含模糊纬度的记录，但经度则是完全错误，因为它们位于大西洋、西非南部。看起来它们的经度被标记为0了。在很多数据库中你会发现一些本意为“无数据”的地方的值被标记为“0”。函数gbif（在使用默认参数时）将坐标（0，0）设置为NA，但在这坐标中只有一个值为0的时候则不。让我们通过查询记录找出经度为0的记录吧。
我们来看看这些记录：
lonzero = subset(acgeo, lon==0) # show all records, only the first 13 columns lonzero[, 1:13] ## species continent country adm1 adm2 ## 1159 Solanum acaule Bitter subsp. acaule <NA> Argentina <NA> <NA> ## 1160 Solanum acaule Bitter subsp. acaule <NA> Bolivia <NA> <NA> ## 1161 Solanum acaule Bitter subsp. acaule <NA> Peru <NA> <NA> ## 1162 Solanum acaule Bitter subsp. acaule <NA> Peru <NA> <NA> ## 1163 Solanum acaule Bitter subsp. acaule <NA> Argentina <NA> <NA> ## 1164 Solanum acaule Bitter subsp. acaule <NA> Bolivia <NA> <NA> ## locality lat lon ## 1159 between Quelbrada del Chorro and Laguna Colorada -23.716667 0 ## 1160 Llave -16.083334 0 ## 1161 km 205 between Puno and Cuzco -6.983333 0 ## 1162 km 205 between Puno and Cuzco -6.983333 0 ## 1163 between Quelbrada del Chorro and Laguna Colorada -23.716667 0 ## 1164 Llave -16.083334 0 ## coordUncertaintyM alt institution collection catalogNumber ## 1159 <NA> 3400 IPK GB WKS 30027 ## 1160 <NA> 3900 IPK GB WKS 30050 ## 1161 <NA> 4250 IPK WKS 30048 304709 ## 1162 <NA> 4250 IPK GB WKS 30048 ## 1163 <NA> 3400 IPK WKS 30027 304688 ## 1164 <NA> 3900 IPK WKS 30050 304711
这些记录的地点来自于玻利维亚、秘鲁和阿根廷，证明坐标存在错误。另外，也有可能坐标是正确的，其指代的位于大西洋的可能是一条鱼被抓住的位置而不是采集S.acaule的地点。一条有着错误的物种名的记录是最难被纠正的（例如，区分熊与野人，Lozier et al., 2009）。一条位于厄瓜多尔的记录就是这样，这其中就有一些关于它到底是一个S.albicans的标本还是一个S.acaule的异常六倍体变种的争论。
2.4 重复记录
有趣的是，另一个数据质量问题被揭示：每条‘longzero’的记录出现两次。有可能是因为植物样本通常可以分割开来并送到多个标本馆。但在本例中似乎是IPK（The Leibniz Institute of Plant Genetics and Crop Plant Research)将这些数据向GBIF提供了两次（也许来自IPK的独立数据库？）。函数‘duplicated’有时被用于移除这些重复。
# which records are duplicates (only for the first 10 columns)? dups <- duplicated(lonzero[, 1:10]) # remove duplicates lonzero <- lonzero[dups, ] lonzero[,1:13] ## species continent country adm1 adm2 ## 1162 Solanum acaule Bitter subsp. acaule <NA> Peru <NA> <NA> ## 1163 Solanum acaule Bitter subsp. acaule <NA> Argentina <NA> <NA> ## 1164 Solanum acaule Bitter subsp. acaule <NA> Bolivia <NA> <NA> ## locality lat lon ## 1162 km 205 between Puno and Cuzco -6.983333 0 ## 1163 between Quelbrada del Chorro and Laguna Colorada -23.716667 0 ## 1164 Llave -16.083334 0 ## coordUncertaintyM alt institution collection catalogNumber ## 1162 <NA> 4250 IPK GB WKS 30048 ## 1163 <NA> 3400 IPK WKS 30027 304688 ## 1164 <NA> 3900 IPK WKS 30050 304711
另一种方法可能是检测数据中一些物种相同且坐标相同的重复，即使这些记录是被不同的人在不同的年份收集的。（在我们的例子中，使用物种是多余的，因为我们只有一个物种给的数据）。
# differentiating by (sub) species
# dups2 <- duplicated(acgeo[, c('species', 'lon', 'lat')])
# ignoring (sub) species and other naming variation
dups2 <- duplicated(acgeo[, c('lon', 'lat')])
# number of duplicates
sum(dups2)
## [1] 483
# keep the records that are _not_ duplicated
acg <- acgeo[!dups2, ]
我们把巴基斯坦附近的记录返还给阿根廷吧，然后去掉巴西、南极洲以及经度=0的记录。
i <- acg$lon > 0 & acg$lat > 0 acg$lon[i] <- -1 * acg$lon[i] acg$lat[i] <- -1 * acg$lat[i] acg <- acg[acg$lon < -50 & acg$lat > -50, ]
2.5 交叉检验
通过视觉和其他方式对坐标进行交叉检验是重要的。一种是使用坐标所意味的国家，通过比对由记录指定的地点所位于的国家（或更低级别的行政区划）（Hijmans et al., 1999）。在下方的例子中我们利用sp包中的coordinates函数去创建了一个SpatialPointsDataFrame，然后使用同样来自于sp包的over函数作了一个基于国家图形的图形中的点查询。
我们使用统计功能符号（波浪线）制作了一个SpatialPointsDataFrame。
library(sp) coordinates(acg) <- ~lon+lat crs(acg) <- crs(wrld_simpl) class(acg) ## [1] "SpatialPointsDataFrame" ## attr(,"package") ## [1] "sp"
现在我们可以使用坐标来对wrld_simpl（一个SpatialPolygonsDataFrame）中的图形进行空间查询。
class(wrld_simpl) ## [1] "SpatialPolygonsDataFrame" ## attr(,"package") ## [1] "sp" ovr <- over(acg, wrld_simpl)
对象'ovr'对于每一点都具有来自wrld_simpl的匹配记录。我们需要在wrld_simpl的data.frame中使用变量'NAME'。
head(ovr) ## FIPS ISO2 ISO3 UN NAME AREA POP2005 REGION SUBREGION LON ## 1 AR AR ARG 32 Argentina 273669 38747148 19 5 -65.167 ## 2 PE PE PER 604 Peru 128000 27274266 19 5 -75.552 ## 3 AR AR ARG 32 Argentina 273669 38747148 19 5 -65.167 ## 4 BL BO BOL 68 Bolivia 108438 9182015 19 5 -64.671 ## 5 BL BO BOL 68 Bolivia 108438 9182015 19 5 -64.671 ## 6 BL BO BOL 68 Bolivia 108438 9182015 19 5 -64.671 ## LAT ## 1 -35.377 ## 2 -9.326 ## 3 -35.377 ## 4 -16.715 ## 5 -16.715 ## 6 -16.715 cntr <- ovr$NAME
我们应该问这两个问题：（1）哪个点（由它们的记录编号确定）没有匹配任何一个国家（那就是它们在海洋里）？（没有（因为我们已经将映射在海洋中的点移除了））。（2）哪个点的坐标显示它们在与gbif记录中列出的country字段所不同的国家。
i <- which(is.na(cntr)) i ## integer(0) j <- which(cntr != acg$country) # for the mismatches, bind the country names of the polygons and points cbind(cntr, acg$country)[j,] ## cntr ## [1,] "27" "Argentina" ## [2,] "172" "Bolivia" ## [3,] "172" "Bolivia" ## [4,] "172" "Bolivia"
本例中所产生的不匹配可能是因为wrld_simple并不非常精确而这些地图上的记录地点过于靠近玻利维亚和它的邻国边界。
plot(acg) plot(wrld_simpl, add=T, border='blue', lwd=2) points(acg[j, ], col='red', pch=20, cex=2)
查看sp包中有关over函数的更多信息。我们在上述例子中使用的wrld_simpl图形不是很精确，而且可能不应该用在真正的分析中。在GADM上查看更详细的行政区划文件，或者使用raster包中的getData函数（例：getData('gadm', country='BOL', level=0)可以得到玻利维亚的国界；getData('countries')去获得所有国家边界）。
2.6地图配准
如果你有一些包含位置描述但无坐标信息的记录，你应该考虑的是对它们进行地图配准。不是所有的记录都能被地图配准的。有时甚至国家一栏都是unknown（country=="UNK"）。现在我们仅选择一些没有坐标但有位置描述的记录。
georef <- subset(acaule, (is.na(lon) | is.na(lat)) & ! is.na(locality) ) dim(georef) ## [1] 131 25 georef[1:3,1:13] ## species continent country adm1 adm2 ## 606 solanum acaule acaule BITTER <NA> Bolivia <NA> <NA> ## 607 solanum acaule acaule BITTER <NA> Peru <NA> <NA> ## 618 solanum acaule acaule BITTER <NA> Peru <NA> <NA> ## locality ## 606 La Paz P. Franz Tamayo Viscachani 3 km from Huaylapuquio to Pelechuco ## 607 Puno P. San Roman Near Tinco Palca ## 618 Puno P. Lampa Saraccocha ## lat lon coordUncertaintyM alt institution collection ## 606 NA NA <NA> 4000 PER001 CIP - Potato collection ## 607 NA NA <NA> 4000 PER001 CIP - Potato collection ## 618 NA NA <NA> 4100 PER001 CIP - Potato collection ## catalogNumber ## 606 CIP-762165 ## 607 CIP-761962 ## 618 CIP-762376
对于地图配准，你可以尝试使用dismo包中的geocode函数向Google API发送请求。我们在下方演示，但这种用法通常是不推荐的，因为对于精确的地图配准，你需要一个详尽的地图接口，在理想情况下，你可以捕获与每个地图配准相关的不确定性（Wieczorek et al., 2004）。
以下是其中一条经度=0的记录的示例，使用Google的地理编码服务。我们把这个函数放入一个‘try’函数中，去确保在电脑无法连接到互联网时有一个简洁的错误处理。注意我们使用了“cloc”（concatenated locality)域。
georef$cloc[4] ## [1] "Ayacucho P. Huamanga Minas Ckucho, Peru" b <- try( geocode(georef$cloc[4]) ) ## Loading required namespace: XML b ## originalPlace interpretedPlace longitude ## 1 Ayacucho P. Huamanga Minas Ckucho, Peru Ayacucho, Peru -74.22356 ## latitude xmin xmax ymin ymax uncertainty ## 1 -13.16387 -74.24456 -74.18004 -13.19747 -13.11926 4359
错误笔记：
> georef$cloc[4]
[1] "Ayacucho P. Huamanga Minas Ckucho, Peru"
> b <- try( geocode(georef$cloc[4]))
载入需要的名字空间：XML
No such file or directoryfailed to load external entity "http://maps.google.com/maps/api/geocode/xml?address=Ayacucho+P.+Huamanga+Minas+Ckucho,Peru&sensor=false"
Error : 1: No such file or directory2: failed to load external entity "http://maps.google.com/maps/api/geocode/xml?address=Ayacucho+P.+Huamanga+Minas+Ckucho,Peru&sensor=false"
Error : 找不到对象'doc'
尝试将ss设置为全局模式后提示：
> b <- try( geocode(georef$cloc[4]))
No such file or directoryfailed to load external entity "http://maps.google.com/maps/api/geocode/xml?address=Ayacucho+P.+Huamanga+Minas+Ckucho,Peru&sensor=false"
Error : 1: No such file or directory2: failed to load external entity "http://maps.google.com/maps/api/geocode/xml?address=Ayacucho+P.+Huamanga+Minas+Ckucho,Peru&sensor=false"
Error in .geocode(xx$place, oneRecord = oneRecord, extent = extent, progress = progress) :
找不到对象'doc'
在使用geocode函数之前最好将记录写入一个表中然后在电子表格中对数据进行“清理”。这种清理包括翻译、扩写缩写、纠正拼写错误以及让确保重复数据完全相同以便它们只能被地图配准一次。然后将这个表重新读入R中，创建唯一的位点，对其进行地图配准后合并到原始数据中。
2.7 抽样偏差
抽样偏差在发现记录中经常出现（Hijmans et al., 2001）。一种方式是可以尝试通过子采样记录来消除一些偏差，在下方有所演示。然而，二次抽样减少了记录数，而且无法更正那些根本没有采样的区域的数据。其遇到的另一个问题是地区密集记录可能实际上是相对适合栖息地的真实反映。如在SDM中的许多步骤，你需要了解有关你数据和物种的一些东西才能把它们实现好。参见Philips et al.(2009)的一种使用MaxEnt的方法去处理一组物种的发生记录的偏差。
# create a RasterLayer with the extent of acgeo r <- raster(acg) # set the resolution of the cells to (for example) 1 degree res(r) <- 1 # expand (extend) the extent of the RasterLayer a little r <- extend(r, extent(r)+1) # sample: acsel <- gridSample(acg, r, n=1) # to illustrate the method and show the result p <- rasterToPolygons(r) plot(p, border='gray') points(acg) # selected points in red points(acsel, cex=1, col='red', pch='x')
图片有误差，右边为官方原图。
需要注意的是使用gridSample函数你同样可以做‘棋盘’抽样。这可用于在‘training’和‘testing’集合中拆分数据。（参见model evaluation章节）。
此时，保存清理好的数据集是很有用的。比如使用函数write.table或是write.csv以便我们在以后使用它们。我们那样做了，这个保存好的文件就可以像这样通过dismo被检索：
file <- paste(system.file(package="dismo"), '/ex/acaule.csv', sep='') acsel <- read.csv(file)
在一个真正的研究项目中你将要在第一次数据清理到完成步骤中花费更多时间，部分在R中进行，也会用到其他程序。
2.8 练习
1. 使用gbif函数下载African elephant（或者其他的你偏爱的物种，尝试获得一个10-100条记录之间的数据集）。使用选项“geo=FALSE”也去获得没有（数字）地图配准的记录。
2. 总结这个数据：共有多少条记录，有多少记录包含坐标，有多少记录没有坐标但有一个文本形式的地理参考（位置描述）？
3. 使用函数‘geocode’对10条以上的无坐标记录进行地图配准。
4. 适用所有这些记录制作一个简单的地图，使用颜色和符号去区分来自gbif的坐标和来自Google的坐标（来自geocode函数）。使用‘gmap'去创建一个底图。
5. 你是否认为这些观察到的意见是对这些物种分布地（和生态位）的真实反映？
更进一步：
6. 使用‘rasterize’函数建立一个观察到的意见的数量的光栅并制作一张地图。国家边界使用maptools包中的“wrld_simpl”来获得。
7. 使用地图配准对不确定相关的数据进行绘制。从gbif返回的数据中有一些记录，你也可以从geocode函数返回的数据中提取它们。
