根据经纬度（从nc格式环境数据文件中）提取环境因子

---

前言

本文主要利用nc格式环境数据文件和物种经纬度分布文件，根据经纬度（从nc格式环境数据文件中）提取环境因子

---

一、准备所需文件

1. nc格式环境数据文件
   本文所用的环境数据文件来自Bio-oracle

2. 物种经纬度分布文件
   

二、代码分享

根据经纬度数据从环境数据文件中提取对应的环境因子

# 加载/下载必要的库
install.packages("ncdf4")library(raster)
library(ncdf4)
# 设置参数
nc_dir <- "D:/oceandata/Bio-ORACLE_download/"          # NetCDF文件目录
output_file <- "C:/Users/www/Desktop/pop_gps_all_env.txt" # 输出文件
gps_file <- "C:/Users/www/Desktop/pop.gps.txt"  #物种经纬度分布文件# 读取GPS数据
pop_gps <- read.table(gps_file, header = TRUE)# 获取所有.nc文件路径
nc_files <- list.files(path = nc_dir, pattern = "\\.nc$", full.names = TRUE)# 创建进度条
pb <- txtProgressBar(min = 0, max = length(nc_files), style = 3)# 循环处理每个nc文件
for(i in seq_along(nc_files)){tryCatch({# 读取当前nc文件current_raster <- raster(nc_files[i])# 坐标系检查与转换if (!grepl("+proj=longlat", crs(current_raster))) {current_raster <- projectRaster(current_raster, crs = "+proj=longlat +datum=WGS84")}# 从文件名提取变量名var_name <- tools::file_path_sans_ext(basename(nc_files[i]))# 提取环境数据（带双线性插值）env_values <- extract(current_raster, pop_gps[, c("Longitude", "Latitude")],method = "bilinear")# 检查NA值并尝试用最近邻插值填充na_indices <- which(is.na(env_values))if (length(na_indices) > 0) {env_values[na_indices] <- extract(current_raster, pop_gps[na_indices, c("Longitude", "Latitude")],method = "simple")}# 添加到数据框pop_gps[[var_name]] <- env_values# 更新进度条setTxtProgressBar(pb, i)}, error = function(e){message(sprintf("\n文件 %s 处理失败: %s", nc_files[i], e$message))})
}# 关闭进度条
close(pb)# 保存结果
write.table(pop_gps, output_file, row.names = FALSE, sep = "\t")# 结果验证
cat("\n处理完成！共添加", length(nc_files), "个环境变量\n")
cat("最终数据维度：", dim(pop_gps), "\n")
cat("NA值统计：\n")
print(colSums(is.na(pop_gps)))

输出文件示例：pop_gps_all_env.txt

总结

– 2020-8-6