当前位置：首页 > news >正文

【环境数据处理】-基于R批量对环境数据克里金插值提高数据精度

news 2025/11/14 5:46:49

# 加载必要的包
if (!require("sp")) install.packages("sp")
if (!require("gstat")) install.packages("gstat")
if (!require("raster")) install.packages("raster")
if (!require("automap")) install.packages("automap")library(sp)
library(gstat)
library(raster)
library(automap)# 设置工作目录和文件路径
input_dir <- "F:/西安科硕信息/2.环境数据已处理/CurrentALL/4"
output_dir <- "F:/西安科硕信息/2.环境数据已处理/CurrentALL/3"# 创建输出目录（如果不存在）
if (!dir.exists(output_dir)) {dir.create(output_dir, recursive = TRUE)
}# 获取所有ASC文件
asc_files <- list.files(path = input_dir, pattern = "\\.asc$", full.names = TRUE)# 检查是否找到文件
if (length(asc_files) == 0) {stop("在指定目录中未找到ASC文件")
}cat("找到", length(asc_files), "个ASC文件\n")# 处理每个ASC文件
for (i in 1:length(asc_files)) {asc_file <- asc_files[i]filename <- basename(asc_file)output_name <- sub("\\.asc$", "_kriged.tif", filename)output_path <- file.path(output_dir, output_name)cat("正在处理文件:", filename, "\n")tryCatch({# 1. 读取ASC栅格数据raster_data <- raster(asc_file)cat("  原始像元大小:", res(raster_data), "\n")cat("  原始行列数:", nrow(raster_data), "x", ncol(raster_data), "\n")# 2. 将栅格转换为空间点数据框points_data <- as(raster_data, "SpatialPointsDataFrame")names(points_data) <- "value"# 移除NA值points_data <- points_data[!is.na(points_data$value), ]# 3. 自动克里金插值cat("  正在进行克里金插值...\n")# 定义新的更高分辨率网格（像元大小缩小一半）new_res <- res(raster_data) / 2new_extent <- extent(raster_data)# 创建新的网格grd <- expand.grid(x = seq(from = new_extent@xmin + new_res[1]/2, to = new_extent@xmax - new_res[1]/2, by = new_res[1]),y = seq(from = new_extent@ymin + new_res[2]/2, to = new_extent@ymax - new_res[2]/2, by = new_res[2]))coordinates(grd) <- ~x + ygridded(grd) <- TRUE# 执行自动克里金插值kriging_result <- autoKrige(value ~ 1, points_data, grd)# 4. 将结果转换为栅格kriged_raster <- raster(kriging_result$krige_output["var1.pred"])cat("  插值后像元大小:", res(kriged_raster), "\n")cat("  插值后行列数:", nrow(kriged_raster), "x", ncol(kriged_raster), "\n")# 5. 设置相同的CRScrs(kriged_raster) <- crs(raster_data)# 6. 导出为TIFF文件writeRaster(kriged_raster, output_path, format = "GTiff", datatype = "FLT4S", overwrite = TRUE)cat("  成功导出:", output_name, "\n\n")}, error = function(e) {cat("  处理文件时出错:", filename, "-", e$message, "\n\n")})
}cat("所有文件处理完成！\n")
cat("输出文件保存在:", output_dir, "\n")# 可选：显示处理结果的统计信息
if (length(asc_files) > 0) {cat("\n处理结果统计:\n")for (i in 1:length(asc_files)) {asc_file <- asc_files[i]filename <- basename(asc_file)output_name <- sub("\\.asc$", "_kriged.tif", filename)output_path <- file.path(output_dir, output_name)if (file.exists(output_path)) {result_raster <- raster(output_path)cat(filename, "-> 像元大小:", res(result_raster), "行列数:", nrow(result_raster), "x", ncol(result_raster), "\n")}}
}