根據經緯度(從nc格式環境數據文件中)提取環境因子
文章目錄
- 前言
- 一、準備所需文件
- 二、代碼分享
- 總結
前言
本文主要利用nc格式環境數據文件和物種經緯度分布文件,根據經緯度(從nc格式環境數據文件中)提取環境因子
一、準備所需文件
-
nc格式環境數據文件
本文所用的環境數據文件來自Bio-oracle -
物種經緯度分布文件
二、代碼分享
根據經緯度數據從環境數據文件中提取對應的環境因子
# 加載/下載必要的庫
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