Species Distribution Models
1.写在前面
我为什么想要写这个内容呢? 其实我已经想了很久了,但是一直没有把所有的代码给理清楚,此外这将是我之后研究的主要方向。我对这个方向已经有了初步的了解,看到的满是绝望和不甘。就我个人看法,这个方向已经过了它的热潮,目前正处于低谷期,无人问暇。因此,今年上半年我曾经有一段时间非常痛苦,辗转反侧,我不知道方向还有研究的价值或者必要吗?已经有很多人做过这个方向了,我再来做还有用吗?能提升吗?实际上,可能大多数人只是把这个作为一个应用的工具(水论文),而不是深入的探究,因此我认为只是做应用层面的,这个方向已经饱和了,而我们想要继续深挖下去,就得更加了解其原理和意义,与实际情况相结合才能用更好的产出(如生物多样性保护等)。目前来说物种分布模型(Species Distribution Models,SDMs)已经有很多人在研究了,包括国内的中科院动物所乔慧捷老师、中科院青藏高原所的郭彦龙老师、以及法国科学院的Wilfried Thuiller教授(Biomod2包的发明者)。
现在,我把我的所有SDMs的构建流程进行记录,希望对相关研究方向的学者有所助益。
2.物种分布模型介绍
目前常用的SDMs有单一模型和集成模型,考虑到单个模型的优点和缺点,使用整体模型可能更安全。其中,Biomod2平台是目前常用的集成平台,它包括了十个不同类型的模型:
- 三个回归模型:广义线性模型(generalized linear model, GLM);多元自适应回归样条(multivariate adaptive regression splines, MARS),广义加法模型(generalized additive models, GAM);
- 五个机器学习模型:人工神经网络(artificial neural network, ANN),最大熵(maximum entropy, MaxEnt),随机森林(random forest, RF)和广义助推模型(generalized boosting model, GBM),分类树分析(classification tree analysis, CTA)、
- 一个分类模型:灵活判别分析(flexible discriminant analysis, FDA)
- 一个范围包络:表面分布区分室模型(surface range envelope, SRE)
如果你想继续了解SDMs的研究方向和基本知识,我推荐你阅读以下几篇 综诉类文献:
[1] 郭彦龙,赵泽芳,乔慧捷,等.物种分布模型面临的挑战与发展趋势[J].地球科学进展,2020,35(12):1292-1305.
[2] 刘晓彤,袁泉,倪健.中国植物分布模拟研究现状[J].植物生态学报,2019,43(04):273-283.
[3] 许仲林,彭焕华,彭守璋.物种分布模型的发展及评价方法[J].生态学报,2015,35(02):557-567.
3.输入数据准备及预处理
首先我们需要知道构建SDMs需要哪些要素:物种分布点数据(Species occurrence data)和环境变量数据(Environmental data),分布点数据可从几个地方进行获取,首先是全球生物多样性网络(GBIF)、中国数字植物标本馆(CVH)、NSII-中国国家标本资源平台以及自己的野外调查数据。
3.1.如何从GBIF网站上获取分布点数据(基于rgbif包)
##需要改四个地方的物种名 ----
setwd("D:/SDMs/Thesis/Data/")
library(mapdata)
library(sp)
library(sf)
library(raster)
library(ggplot2)
library(CoordinateCleaner)
library(rgbif)
library(dplyr)
library(beepr)# 这个包是用来进行声音提示的
#第一步:使用''rgbif'' 包从 GBIF 中搜索与目标物种小叶栎(Quercus chenii)相关的数据------
key <- name_suggest(q="Quercus chenii", rank='species')$data$key[1]
key #物种识别号:2879172
occ_count(taxonKey=key, georeferenced=TRUE,
basisOfRecord = "PRESERVED_SPECIMEN"
)
# 266,目前,GBIF 数据库中有 266 条基于博物馆的记录是具有经纬度信息的。
# 需要注意的是,因为数据库一直在更新,这个数字可能会变动,
# 因此使用本流程过程中可能会有部分数值有一些小范围的变动。
data <- occ_search(scientificName = "Quercus chenii", limit= 10000,
hasCoordinate = TRUE, basisOfRecord= "PRESERVED_SPECIMEN"
# ,country = "CN" # 将物种限定在中国境内
) #这一步需要 1-2 分钟完成,具体看网速和物种分布数量。
beepr::beep(8) # 上面代码运行结束会有声音提示
str(data$data)
data$meta # 查看有多少个物种分布点
names(data$data) #检查数据的列名,一共有162列数据
head(data$data[,c("countryCode","country")]) # 例如:检查 54 列和 57 列的列名
datasel<-data$data %>%
dplyr::select(species, decimalLongitude,
decimalLatitude,country,coordinateUncertaintyInMeters) %>%
rename(Species = species, Long = decimalLongitude, Lat =
decimalLatitude, Uncertain=coordinateUncertaintyInMeters) %>%
mutate (Species = as.factor(Species), country=
as.factor(country))
datasel # 现在有了一个包含266 行、5 列数据:Species、Long、Lat、country、Uncertain
datasel %>% group_by(country) %>% summarize(total= n()) %>%
arrange(desc(total)) # GBIF 检索的数据集中每个国家的记录数量
# write.csv(datasel, "Quercus acutissima Carruth.csv")
#第二步:使用 CoordinateCleaner 包删除有问题的分布记录--------------------------
map('worldHires', col=1:10)
# plot(world_map, col="tan2",bg="lightblue", border= NA) #图示地图。
points(datasel$Long, datasel$Lat, bg="red", col='white', pch=21, cex=1.0)
#将从 GBIF 下载的数据映射到地图上 (图 9)。
datasel %>% filter(Long==0, Lat==0) # 统计经纬度等于零的记录
problem_records<- clean_coordinates(datasel, lon = "Long", lat = "Lat",
species = "Species",
tests = c("capitals", "centroids",
"equal","institutions",
"zeros", "seas"))
# 检测记录点坐标是否围绕首都、国家的中心
# 是否落入海洋,为零,或在饲养动物的博物馆(机构)周围。
summary(problem_records) # 2个有问题的分布点被发现
datasel_problematic <- datasel[which(problem_records$.summary==
"FALSE"),] #只选择有问题的点。
# 查看有问题分布点的位置
map('worldHires', col=1:10)
points(datasel_problematic$Long, datasel_problematic$Lat, bg="red",
col='white', pch=22, cex=0.8) # 将有问题的点映射到地图上
# 从数据集中排除所有有问题的记录
datasel_clean <- datasel[which(problem_records$.summary== "TRUE"),]
str(datasel_clean) # 现在还剩下264个正常分布点数据
# 删除坐标中具有较高不确定性的记录。
# 看看有多少个体有这样的信息,并把值的范围图像化
summary(is.na(datasel_clean$Uncertain))
hist(datasel_clean$Uncertain/1000, breaks = 10,
main="Coordinate uncertainy in km", xlab= "Uncertainty in km")
# 注意,GBIF 提供的信息是以米为单位
# 但是为了评估不确定性,使用千米为单位更加方便。 (图 13)。
# 这个命令中用于选择“信息不确定”的点,或者是数值<10 km 以下的点。
datasel_cleanF<-datasel_clean %>% filter(is.na(Uncertain)| Uncertain/1000 < 10)
#最后一步是删除重复的记录。
datasel_cleanF <- datasel_cleanF %>% distinct(Long, Lat,.keep_all=T)
str(datasel_cleanF) # 图 14。
#以下代码可以图示这个最终数据集。
map('worldHires', col=1:10)
points(datasel_cleanF$Long, datasel_cleanF$Lat, bg="blue", col='white',
pch=21, cex=0.5) # 图 15
# 最后记得将数据保存下来
write.csv(datasel_cleanF, "Quercus chenii_gbif.csv")
如果你想要从CVH网站上获得分布点数据,就需要自己去查看里面的标本记录,然后手动记录下来,这个过程可能比较繁琐,目前没有相关的包能够获取上面的记录。
3.2.分布点稀疏处理(基于spThin函数)
为了防止数据空间自相关性对模拟结果的影响,通常需要对分布点数据进行稀疏处理,如10km范围内只保留一个分布点数据。目前常用的方法是在ArcGIS中安装SDMtoolbox (v.2.5) 小插件,并使用其中所提供的 “Spatially Rarefy Occurence Data” 功能进行处理。但是采用这样的方式有时候会报让人摸不着头脑的错误,所以我提供了一个使用R语言进行稀疏处理的方法,如下:
## thin----------------------------------------------
setwd("C:\\Users\\YP\\Desktop\\SDMData")
# library(wallace)
library(spThin) # v.0.2.0
library(readxl)
library(mapdata)
# run_wallace() ##运行run后可直接进入到浏览器界面
QuercusL <- read.csv("Quercus chenii_gbif.csv", sep = ",")
head(QuercusL)
summary(QuercusL)
# remove rows with duplicate coordinates
occs.dups <- duplicated(QuercusL[c('Long', 'Lat')])
occs <- QuercusL[!occs.dups,]
# remove NAs
occs <- occs[complete.cases(occs$Long, occs$Lat), ]
summary(occs)
output <- spThin::thin(occs, 'Lat', 'Long', 'Species',
thin.par = 10, reps = 100, # 稀疏距离参数thin.par设置为10km,并重复运行100次
locs.thinned.list.return = TRUE,
write.files = FALSE,
verbose = TRUE)
plotThin(output)
# 基于最优迭代构建最优输出
maxThin <- which(sapply(output, nrow) == max(sapply(output, nrow)))
# if there's more than one max, pick the first one
maxThin <- output[[ifelse(length(maxThin) > 1, maxThin[1], maxThin)]]
# subset occs to match only thinned occs
occs <- occs[as.numeric(rownames(maxThin)),]
write.csv(occs,"./Quercus chenii_gbif_occs.csv")
# plot
map('worldHires', 'China')
points(occs$Long, occs$Lat, bg="red", col='white', pch=21, cex=1.0)
分布点空间位置:
3.3.如何获取环境变量数据(基于getData函数)
目前,常用的大尺度空间数据主要由WorldClim提供,我之前写了一篇关于如何下载该数据的方法,感兴趣的可以看看:Worldclim(v1.4、v2.1)数据集使用介绍
在这里我可以提供另一种方法:使用R语言下载数据。代码如下:
##使用"Raster"包从 WorldClim 下载环境数据--------------------------------
# 目前该函数只提供了下载2.5, 5, 和 10 min的空间分辨率下的数据,且来源于worldclim1.4(CMIP5)
bioclim<-getData('worldclim', var="bio", res=5) #把 19 个5min分辨率的气候数据图层下载到本地
Future_bioclim2<-getData('CMIP5', var='tmin', res=10, rcp=85, model='AC', year=70) # 获取未来气候数据
names(bioclim) #查看存储在"bioclim"对象中的不同气候变量层的名称
我可能不太建议使用R代码来下载数据,因为这可能受到网速的限制,导致下载突然中断,所以直接到官网里面下载,这样也很方便。
当然还可以使用CHELSA数据集,这是由瑞士森林雪和景观研究所(WSL)提供的。现在也有很多研究使用该数据:
3.4.环境变量多重共线性分析
在构建SDMs时,会存在不同环境变量对分布分布点的贡献程度相似,即可看做这几个变量对物种分布的作用相同,因此为了降低模型拟合的不确定性,就需要从这几个变量中选择一个贡献度最高,最关键的环境变量,因此就需要使用多重共线性分析。常用的方法为Pearson相关性系数以及方差膨胀因子(VIF)。其中Pearson相关性系数和VIF计算方法如下:
library(raster)
library(pheatmap)
library(usdm)
setwd("C:/Users/YP/Desktop/SDMData")
quercus <- read.csv("Quercus chenii_gbif_occs.csv", sep = ",")
summary(quercus)
current_folder1 <- "C:/Users/YP/Desktop/Data/"
# 气候因子1
current_bio <- stack(paste0(current_folder1,
c("wc2.1_30s_bio_1","wc2.1_30s_bio_2",
"wc2.1_30s_bio_3","wc2.1_30s_bio_4",
"wc2.1_30s_bio_5","wc2.1_30s_bio_6",
"wc2.1_30s_bio_7","wc2.1_30s_bio_8",
"wc2.1_30s_bio_9","wc2.1_30s_bio_10",
"wc2.1_30s_bio_11","wc2.1_30s_bio_12",
"wc2.1_30s_bio_13","wc2.1_30s_bio_14",
"wc2.1_30s_bio_15","wc2.1_30s_bio_16",
"wc2.1_30s_bio_17","wc2.1_30s_bio_18",
"wc2.1_30s_bio_19"
),
".tif")
)
# 修改变量名
new_names1 <- paste0("bio", 1:19)
names(current_bio) <- new_names1
# 地形因子2
current_topo <- stack(paste0(current_folder1,
c("DEM_2023", "Aspect", "Slope"),
".tif")
)
new_names2 <- c("Elevation", "Aspect", "Slope")
names(current_topo) <- new_names2
# 土壤因子3
current_soil <- stack(paste0(current_folder1,
c("BD1","PH1", "TK1",
"TN1","TP1","soil_1km1"),
".tif")
)
new_names3 <- c("BD","pH","TK","TN","TP","SoilType")
names(current_soil) <- new_names3
# 将以上不同变量数据合并为一个叠加图层
(current <- stack(current_bio, current_topo, current_soil))
# 进行坐标转换,将分布点数据转换为环境数据空间投影格式,即WGS1984
longlat <- as.data.frame(quercus[, c(4:5)]) # Long Lat
coordinates(longlat) <- c("Long", "Lat") # 这里只能是Long在前,Lat在后,不然下一行代码会报错!!!
proj4string(longlat) <- CRS("+proj=longlat +datum=WGS84")
UTM_Proj <- CRS(st_crs(current)$proj4string)
xy <- spTransform(longlat, UTM_Proj)
# 检查坐标系是否一致
if (st_crs(current)$proj4string == st_crs(xy)$proj4string) {
print("坐标系转换结果:坐标系相同")
} else {
print("坐标系转换结果:坐标系不同")
}
# 从所有的记录中提取环境因子对应值
data <- raster::extract(current, xy)
data_na <- na.omit(data)
# 构建pearson相关分析(pearson>0.80)
cor_matrix <- cor(data_na, method = "pearson") # R在带stats包
cor_matrix
# dev.off() # 聚类热图绘制报错时可以使用,这可以清除之前可能干扰当前绘图的图形设置
pheatmap(cor_matrix,cluster_cols = F,display_numbers = T,main = "环境因子相关性分析")
# 计算方差膨胀因子VIF
vifstep(data_na) # usdm包:vifstep()
多重共线性诊断VIF: 一般文章的材料方法中对VIF的限定不同,有的限定在5以内,有的限定在10以内,也有的限定在20以内,这个根据自己的判断来定,没有十分明确的数值,数值越小越好,一般建议保留在VIF<10的变量。
文章:
-
1. Analysis of the distribution pattern of Chinese Ziziphus jujuba under climate change based on optimized biomod2 and MaxEnt models. Ecological Indicators:
-
2. Predicting the invasive trend of exotic plants in China based on the ensemble model under climate change: A case for three invasive plants of Asteraceae. Science of The Total Environment:
-
3. Projected degradation of quercus habitats in southern China under future global warming scenarios. Forest Ecol. Manag.:
-
4. Incorporating eco-evolutionary information into species distribution models provides comprehensive predictions of species range shifts under climate change. Science of The Total Environment:
Pearson相关性热图:
在得到Pearson相关性系数之后,如果我们要去除Pearson>0.80的变量,那么我们应该选择这两个变量中好解释的变量(不要为难自己,给自己挖坑),或者是你认为更加重要的变量(你就是这个领域的专家,你的知识就是最正确的)。
4.Biomod2模型构建(v.3.5.1)
目前最新的版本应该是v.4.2-1,最新版可能使用起来更加流畅和友好,但是我比较习惯使用v.3.5.1版本构建模型,这个是已经成熟使用过的版本,可能报错少一些,运行更加稳定。