Biomod1
- 1.栅格数据处理
- 1.1 读取一个栅格图片
- 1.2 计算数据间的相关系数
- 1.3 生成多波段的栅格图像
- 1.4 修改变量名称
- 1.4.1 计算多个变量之间的相关性
- 2. 矢量数据处理
- 2.1 提取矢量数据
- 2.2 数据掩膜
- 2.2 栅格计算
- 2.3 拓展插件的使用
- 3. 图表绘制
- 3.1 遥感影像绘制
- 3.2 柱状图分析图绘制
- 3.3 在图表中生成标准线（平均值或期望值）
- 3.4 箱线图绘制
- 3.4.1 图像优化
- 4. 线性回归模型
今天学习了一下物种分布模型方面的知识，现在进行一个总结和回顾。
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)
