WEEK5

栅格快速采样

library(sf)

library(terra)

library(pbapply)

# 设置terra包使用的并行处理器核心数

terraOptions(parallel=4)

# 选择文件夹

tif_folder <- choose.dir(default = "", caption = "Select the folder containing the TIF files")

# 读取点矢量数据集

points <- st_read("D:/GISstudy/Thesis/samplepoints.shp")

# 如果没有FID字段，可以添加一个

if (!"FID" %in% names(points)) {

points$FID <- seq_len(nrow(points))

}

print(points)

# 列出所有的TIF文件

tif_files <- list.files(tif_folder, pattern = "\\.tif$", full.names = TRUE)

print(tif_files)

# 使用pblapply来遍历所有TIF文件并且展示进度条

samples_list <- pblapply(tif_files, function(tif_file) {

# 读取栅格文件

r <- rast(tif_file)

# 检查栅格是否有CRS，如果没有，则赋予点矢量数据集的CRS

if (is.na(crs(r))) {

crs(r) <- crs(points)

}

# 采样

samples <- extract(r, points)

return(samples[[1]])

}) # extract函数将自动使用terraOptions 设置的并行数

# 从tif_files中提取文件名，用作samples_df的列名

col_names <- sapply(tif_files, function(x) tools::file_path_sans_ext(basename(x)))

# 创建一个空的数据框以存储采样结果

samples_df <- data.frame(matrix(ncol = length(tif_files), nrow = nrow(points)))

# 将采样列表转换成数据框，并为其设置列名

samples_df <- setNames(as.data.frame(do.call(cbind, samples_list)), col_names)

# 为数据框设置行名为FID

rownames(samples_df) <- points$FID

# 写入CSV文件到TIF文件所在的文件夹

csv_file_path <- file.path(tif_folder, "sampled_data.csv")

write.csv(samples_df, csv_file_path, row.names = TRUE)