# 加载所需包
library(raster)
library(ecospat)
library(ggplot2)
library(dplyr)
library(ade4)
# ---------- 数据准备 ----------
setwd("C:/Users/86157/Desktop/run") # 设置工作目录
# 读取气候数据（6个生物气候变量）
climate_files <- list.files(pattern = "bio.*\\.asc$")
climate_stack <- stack(climate_files)
# 读取6个物种的分布数据
species_list <- c("BT", "ABT", "CBBT", "MGBT", "XKBT", "XYYM")
species_data <- lapply(species_list, function(x){
df <- read.csv(paste0(x, ".csv"))
if(all(c("Longitude", "Latitude") %in% colnames(df))){
df %>%
dplyr::select(x = Longitude, y = Latitude) %>%
na.omit()
} else {
stop(paste("列名不匹配，请检查", x, "的CSV文件"))
}
})
# ---------- 关键修改部分 ----------
# 主成分分析（增强数据清洗）
all_points <- do.call(rbind, species_data)
# 提取环境变量并去除NA值
env_data <- extract(climate_stack, all_points)
valid_rows <- complete.cases(env_data)
env_clean <- env_data[valid_rows, ]
all_points_clean <- all_points[valid_rows, ]
# 执行PCA分析（添加数据标准化）
pca <- dudi.pca(env_clean, scannf = FALSE, nf = 2, center = TRUE, scale = TRUE)
# 创建环境背景网格（修正参数设置）
grid <- ecospat.grid.clim.dyn(
glob = pca$li, # 使用PCA得分作为全局背景
glob1 = pca$li, # 与glob相同数据集
sp = pca$li[0,], # 空数据框创建背景网格
R = 100,
th.env = 0.05
)
# ---------- 生态位模型构建（增强数据验证）----------
niche_models <- lapply(species_data, function(sp){
tryCatch({
# 数据清洗与验证
sp_env <- extract(climate_stack, sp)
valid_sp <- complete.cases(sp_env)
sp_clean <- sp[valid_sp, ]
if(nrow(sp_clean) < 5) {
warning("有效记录不足5条，跳过该物种")
return(NULL)
}
# PCA得分投影验证
sp_scores <- suprow(pca, sp_env[valid_sp, ])$li
if(ncol(sp_scores) != 2) stop("PCA维度不匹配")
# 构建生态位模型
model <- ecospat.grid.clim.dyn(
glob = pca$li, # 使用原始PCA空间
glob1 = pca$li,
sp = sp_scores,
R = 100,
th.env = 0.05
)
# 模型验证
if(is.null(model$z)) stop("密度估计失败")
return(model)
}, error = function(e) {
message("模型构建失败: ", conditionMessage(e))
return(NULL)
})
})
# ---------- 生态位重叠分析（增强空值处理）----------
overlap_matrix <- matrix(NA, nrow=6, ncol=6,
dimnames=list(species_list, species_list))
for(i in 1:6){
for(j in 1:6){
if(!is.null(niche_models[[i]]) && !is.null(niche_models[[j]])){
tryCatch({
overlap_matrix[i,j] <- ecospat.niche.overlap(
niche_models[[i]],
niche_models[[j]],
method = "D"
)$D
}, error = function(e) {
message(sprintf("物种%s与%s重叠计算失败: %s",
species_list[i], species_list[j],
conditionMessage(e)))
})
}
}
}
# ---------- 可视化与输出（增强容错）----------
# 仅绘制有效模型
valid_models <- sapply(niche_models, function(x) !is.null(x))
if(sum(valid_models) > 0){
plot_colors <- rainbow(6)[valid_models]
# 创建基础绘图窗口
ecospat.plot.niche(niche_models[[which(valid_models)[1]]],
col=plot_colors[1],
title="生态位分布",
name.axis1="PC1 (主要环境梯度)",
name.axis2="PC2 (次要环境梯度)")
# 添加其他物种分布
if(sum(valid_models) > 1){
for(i in 2:sum(valid_models)){
ecospat.plot.niche(niche_models[[which(valid_models)[i]]],
col=plot_colors[i],
add=TRUE)
}
}
# 添加图例
legend("topright",
legend=species_list[valid_models],
fill=plot_colors,
cex=0.8,
title="物种代码")
} else {
warning("没有有效模型可供绘制")
}
# 结果输出
print("生态位重叠矩阵:")
print(round(overlap_matrix, 3))
write.csv(overlap_matrix, "Niche_Overlap_Matrix.csv")
# 保存工作空间
save.image("Niche_Analysis_Backup.RData")这串代码产生了错误于.local(obj, ...):
cannot derive coordinates from non-numeric matrix，修改并生成完整代码
本文探讨了如何通过df%>%mutate函数应用正则表达式，有效去除DataFrame中所有列的特定模式，提升数据处理效率。
扫一扫
6404
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
