Biomod2 （上）：物种分布模型预备知识总结

今天学习了一下物种分布模型方面的知识，现在进行一个总结和回顾。

1.栅格数据处理

library(raster)
library(tidyverse)
matrix(1:9, nrow = 3, ncol = 3)# 生成一个栅格矩阵
raster(matrix(1:9, 3))
raster(matrix(1:9, 3)) %>% plot()# 生成一个简单的栅格图片
r<-raster(ncol = 30, nrow = 30)
values(r) <- 1:ncell(r) # ncell(r):表示r的栅格个数，共900个
plot(r)

生成一个简单的栅格图像：Raster

自定义一个栅格图像：30*30

1.1 读取一个栅格图片

# 读取栅格图片
raster(system.file("external/rlogo.grd", package = "raster")) %>% plot() #管道函数

这是R语言自带的图片：

1.2 计算数据间的相关系数

# 计算相关系数
cor(runif(100,0,10),runif(100,0,10)) # 提供两个长度一直的向量

0.0046693

1.3 生成多波段的栅格图像

r1 <- raster(matrix(runif(900), 30, 30))
r2 <- raster(matrix(runif(900), 30, 30))
r3 <- raster(matrix(runif(900), 30, 30))
r4 <- raster(matrix(runif(900), 30, 30))
r5 <- raster(matrix(runif(900), 30, 30))
# 将五个栅格图层堆叠在一起，形成多波段影像
stack(r1, r2, r3, r4, r5) 
stack(r1, r2, r3, r4, r5) %>% plot()

结果展示：

1.4 修改变量名称

# 更改变量名称
df <- stack(r1, r2, r3, r4, r5) %>%as.data.frame(xy = T) %>%as_tibble() %>%setNames(c("lon", "lat", "x1", "x2", "x3", "x4", "y")) # 重命名
head(df)

1.4.1 计算多个变量之间的相关性

## 变量间相关性 法1
cor(df[, c("x1","x2","x3","x4","y")])## 变量间相关性，法2
#devtools::install_github("danlwarren/ENMTools") # 注意:大概需要安装46个packages
# 安装的时候会出现一些问题，不太好安装
library(ENMTools) 
raster.cor.matrix(stack(r1, r2, r3, r4, r5))

2. 矢量数据处理

2.1 提取矢量数据

install.packages("mapchina")
library(mapchina)
cd_sf <- mapchina::china %>%dplyr::filter(Name_Perfecture == "成都市") %>%group_by(Name_Province) %>%summarise(geometry = sf::st_union(geometry)) %>%ungroup()
colnames(cd_sf) # see all variable names
plot(cd_sf)

2.2 数据掩膜

# 使用本地数据DEM进行掩膜
china_dem <- raster("D:/Datasets/w001001.adf")
bj_dem <- china_dem %>%crop(bj_sf) %>%mask(bj_sf)
plot(bj_dem)

2.2 栅格计算

类似于一个二分类问题：

res <- stack(r1, r2) %>%calc(x = ., fun = function(x) {ifelse(x[1] > 0.5 & x[2] > 0.5, 0, 1) # 对数据进行二分类处理,# 同时满足返回0，否则返回1})
res
plot(res)

使用case_when函数能解决上面的问题，更方便

# 使用case_when函数
stack(r1, r2) %>%as.data.frame(xy = T) %>% # 将数据转换为数据框格式mutate(res = case_when(layer.1 > 0.5 & layer.2 > 0.5 ~ 0,TRUE ~ 1)) %>%dplyr::select(1, 2, last_col()) %>%as_tibble() %>%plot()getValues(res) %>% head(10)

2.3 拓展插件的使用

更方便查看和统计数据

# 拓展插件
#install.packages("flextable")
library(flextable)
tibble(model = LETTERS[1:10], mean = runif(10, 0.4, 1)) %>%mutate(sd = runif(10, 0.01, 0.1)) %>%mutate(weight = mean / sum(mean)) %>%regulartable() %>%flextable::align(align = "center", part = "all") %>%theme_booktabs()

3. 图表绘制

3.1 遥感影像绘制

library(terra)
library(tidyterra)
library(ggspatial)
ggplot()+geom_spatraster(data = rast(bj_dem)) +scale_fill_whitebox_c(palette = "muted",na.value = "white") +annotation_scale() +annotation_north_arrow(which_north = "grid")+theme()

3.2 柱状图分析图绘制

AUC、TSS展示

# AUC、TSS绘制
tibble(model1 = LETTERS[1:10], AUC = runif(10, 0.6, 1), TSS = runif(10, 0.6, 1)) %>% ggplot()+geom_col(aes(model1,AUC))##---------------------------------------------------------------------
tibble(model = LETTERS[1:10], AUC = runif(10, 0.6, 1), TSS = runif(10, 0.6, 1)) %>%pivot_longer(-model,names_to = "vars",values_to = "value") %>% ggplot()+geom_col(aes(model,value,fill = vars),position = 'dodge')+ggsci::scale_fill_lancet()+labs(x = "建模算法",y = "AUC/TSS",fill = "")+theme_light()

3.3 在图表中生成标准线（平均值或期望值）

# 生成选择线
install.packages("extrafont")
install.packages("showtext")
library(extrafont)
library(showtext)
showtext_auto(enable = TRUE)
font_add("Times", regular ="./font/Times New Roman.ttf") # 字体选择
font_add("KaiTi", regular ="./font/KaiTi.ttf")
font_add("wqy", regular = "wqy-microhei.ttc")set.seed(123)
tibble(model = LETTERS[1:10], AUC = runif(10, 0.6, 1), TSS = runif(10, 0.6, 1)) %>%pivot_longer(-model,names_to = "vars",values_to = "value") %>% ggplot()+geom_col(aes(model,value,fill = vars),position = 'dodge')+geom_hline(yintercept = 0.6,linetype = 2,size = 0.5)+ # 选择线的位置，即y的大小geom_hline(yintercept = 0.75,linetype = 1.5,size = 1)+ggsci::scale_fill_lancet()+labs(x = "建模算法",y = "AUC/TSS",fill = "")+theme_light()+theme(axis.text.x = element_text(size = 13, family = "KaiTi",angle = 45),axis.title.x = element_text(size = 15,family = "KaiTi"),legend.text = element_text(size = 11,family = "Times"),axis.text.y = element_text(size = 11,family = "Times"),axis.title.y = element_text(size = 13,family = "Times"),legend.position = "bottom")

3.4 箱线图绘制

这是一个非常好用的package：DT，能够更方便的显示数据表

tibble(vars = paste0("vars",rep(1:5,50)),model = rep(LETTERS[1:5],each = 50),weight = runif(250,0.5,1),run = rep(rep(1:5,each = 5),10)) %>% DT::datatable()  # 调用DT中的datatable函数显示数据

tibble(vars = paste0("vars",rep(1:5, 50)), # ,1到5重复50次，并与vars粘贴model = rep(LETTERS[1:5],each = 50),weight = runif(250,0.5,1),run = rep(rep(1:5,each = 5),10)) %>% ggplot()+geom_boxplot(aes(vars,weight,fill = model),alpha = 0.8)+theme_bw()+theme(panel.grid = element_blank())+theme(text=element_text(size=12,  family="serif"))

3.4.1 图像优化

# 图像优化
tibble(vars = paste0("vars",rep(1:5,50)),model = rep(LETTERS[1:5],each = 50),weight = runif(250,0.5,1),run = rep(rep(1:5,each = 5),10)) %>% ggplot()+geom_boxplot(aes(vars,weight,fill = model),alpha = 0.8)+ggsci::scale_fill_lancet()+labs(x = "预测变量",y = "贡献度",fill = "")+ggthemes::theme_economist_white()+theme_bw()+theme(axis.text.x = element_text(size = 13, family = "RMN",angle = 45),axis.title.x = element_text(size = 15,family = "KT"),legend.text = element_text(size = 11,family = "RMN"),axis.text.y = element_text(size = 11,family = "RMN"),axis.title.y = element_text(size = 13,family = "KT"),legend.position = "bottom")

4. 线性回归模型

# 线性回归模型
lm_df <- data.frame(x = iris$Sepal.Length,y = iris$Sepal.Width)
lm_model <- lm(data = lm_df,y ~ x)
broom::tidy(lm_model)