library(raster)
library(ade4)
library(ecospat)
library(dplyr)
library(ggplot2)
增强型气候数据加载器（新增投影校验逻辑）
load_climate <- function() {
climate_files <- list.files(pattern = “bio.*\.asc$”)
if(length(climate_files) == 0) stop(“未找到气候文件，请检查：\n1.文件扩展名是否为.asc\n2.文件名是否包含’bio’前缀\n3.工作目录设置正确”)
climate_stack <- stack(lapply(climate_files, function(f){
r <- raster(f)
# 自动修复常见投影问题
if(is.na(proj4string®)) {
message(paste(f, “缺少坐标系，自动设置为WGS84”))
proj4string® <- “+proj=longlat +datum=WGS84”
}
return®
}))
return(climate_stack)
}
强化物种数据加载器（新增坐标转换校验）
load_species <- function(species_list) {
lapply(species_list, function(sp){
csv_file <- paste0(sp, “.csv”)
if(!file.exists(csv_file)) stop(paste(“文件不存在:”, csv_file))
df <- tryCatch({
read.csv(csv_file) %>%
dplyr::select(Longitude, Latitude) %>%
mutate(
Longitude = as.numeric(as.character(Longitude)),
Latitude = as.numeric(as.character(Latitude))
) %>%
na.omit()
}, error = function(e) stop(paste(sp, "数据读取失败:", e$message)))
# 坐标范围二次验证
valid_coords <- df %>%
filter(between(Longitude, -180, 180),
between(Latitude, -90, 90))
if(nrow(valid_coords) < 5) stop(paste(sp, "有效记录不足5条"))
return(as.matrix(valid_coords))
})
}
主分析流程（新增物种数据追踪机制）
main_analysis <- function() {
数据加载
climate_stack <- load_climate()
species_list <- c(“BT”, “ABT”, “CBBT”, “MGBT”, “XKBT”, “XYYM”)
带物种标签的数据加载
species_data <- tryCatch({
data_list <- load_species(species_list)
names(data_list) <- species_list
data_list
}, error = function(e) stop(paste(“数据加载失败:”, e$message)))
构建追踪索引（关键修复点）
all_points <- do.call(rbind, lapply(seq_along(species_data), function(i){
cbind(species_data[[i]], Species = names(species_data)[i])
}))
环境数据提取
env_data <- raster::extract(climate_stack, all_points[,1:2])
valid_idx <- complete.cases(env_data)
if(sum(valid_idx) < 10) stop(“有效环境数据点不足”)
env_clean <- env_data[valid_idx, ]
all_points_clean <- all_points[valid_idx, ]
稳健PCA分析
pca <- dudi.pca(env_clean, scannf = FALSE, nf = 2, center = TRUE, scale = TRUE)
pca_scores <- as.data.frame(pca$li)
colnames(pca_scores) <- c(“PC1”, “PC2”)
按物种分割PCA分数（关键修改）
species_pca <- lapply(species_list, function(sp){
idx <- all_points_clean[,“Species”] == sp
as.matrix(pca_scores[idx, ]) # 确保转换为数值矩阵
})
names(species_pca) <- species_list
环境网格构建（修复坐标类型问题）
niche_models <- lapply(species_list, function(sp){
ecospat.grid.clim.dyn(
glob = as.matrix(pca_scores), # 全局背景
glob1 = as.matrix(pca_scores), # 研究区域背景
sp = species_pca[[sp]], # 物种分布
R = 100,
th.env = 0.05
)
})
names(niche_models) <- species_list
return(list(pca = pca, models = niche_models))
}
生态位模型构建函数
build_models <- function(species_data, pca) {
lapply(seq_along(species_data), function(i){
sp <- species_data[[i]]
tryCatch({
sp_env <- raster::extract(climate_stack, sp)
valid_sp <- complete.cases(sp_env)
sp_clean <- sp[valid_sp, ]
# 增强型数据转换
sp_scores <- tryCatch(
{
scores <- suprow(pca, sp_env[valid_sp, ])$li
as.matrix(scores)
},
error = function(e) stop("PCA投影失败，请检查环境变量一致性")
)
if(nrow(sp_scores) < 5) stop("有效PCA投影点不足")
ecospat.grid.clim.dyn(
glob = pca_scores,
glob1 = pca_scores,
sp = sp_scores,
R = 100,
th.env = 0.05
)
}, error = function(e){
message(sprintf("%s建模失败: %s", species_list[i], e$message))
return(NULL)
})
})
}
niche_models <- build_models(species_data, pca)这段代码产生了错误: 找不到对象’species_data’，修改并生成完整代码
扫一扫
739
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
