WEEK3

消除异常值

library(raster)

library(pbapply)

# 设置工作文件夹路径

work_dir <- "D:/Rlab/thesis/sr900hklcz.tif_mw1"

setwd(work_dir)

# 列出文件夹中所有的 .tif 文件

tif_files <- list.files(pattern = "\\.tif$", full.names = TRUE) # 添加full.names=TRUE以获取完整路径

print(tif_files)

# 创建一个新的子文件夹来存放处理后的文件，如果它还不存在的话

new_dir <- file.path(work_dir, "knocked_out")

if (!dir.exists(new_dir)) {

dir.create(new_dir)

}

# 定义一个函数来敲除特定的值

knockout_value <- function(raster_file, value_to_replace = c(-999, 999), replacement = NA) {

r <- raster(raster_file)

r[r %in% value_to_replace] <- replacement

return(r)

}

# 应用定义好的函数knockout_value到每个.tif文件，并显示进度条

rasters_knocked_out <- pblapply(tif_files, knockout_value)

# 设置投影坐标系

projection <- "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

# 为每个栅格数据定义投影信息

rasters_with_projection <- pblapply(rasters_knocked_out, function(r) {

projection(r) <- projection

return(r)

})

# 保存修改后的栅格文件到新创建的子文件夹，并显示进度条

pblapply(seq_along(rasters_with_projection), function(i) {

outFile <- file.path(new_dir, basename(tif_files[i]))

projected_raster <- projectRaster(rasters_with_projection[[i]], crs = projection)

writeRaster(projected_raster, filename = outFile, format = "GTiff", overwrite = TRUE)

})