如何将R中ggplot生成的克里金插值图导出为指定区域TIFF栅格

已在R中完成普通克里金插值，代码如下：

#ordinary kriging
mydata.krig <- krige(Clay~1, train_data, newdata=mask, vgm)

names(mydata.krig)
names(mydata.krig)[1] <- "Clay.pred"
names(mydata.krig)

min(mydata.krig$Clay.pred, na.rm=T); max(mydata.krig$Clay.pred, na.rm=T) 

ggplot() +
geom_stars(data = mydata.krig["Clay.pred"]) +
scale_fill_gradient(low = "yellow", high = "dark blue", limits = c(15.5486,61.41131)) +
geom_sf(data = thessaly_smus, color = "black", fill = NA, size = 1) +
labs(title = "%Clay Predicted Values - Ordinary Kriging, OK")

已通过ggplot生成预测值地图，需要将插值结果导出为仅保留thessaly_smus多边形区域的TIFF栅格文件，可在GIS软件中打开。

解决方案 #

直接对克里金插值生成的stars对象进行裁剪与导出，是最精准且GIS友好的方式，步骤如下：

1. 加载依赖包 #

确保已安装并加载核心包：

# 未安装则先运行：install.packages(c("stars", "sf"))
library(stars)
library(sf)

2. 裁剪插值结果到目标多边形 #

先裁剪到多边形的外接范围，再筛选出多边形内部的栅格值，将外部值设为NA：

# 第一步：裁剪到多边形外接矩形
clay_cropped <- st_crop(mydata.krig["Clay.pred"], thessaly_smus)

# 第二步：仅保留多边形内部的栅格值，外部设为NA
clay_cropped <- st_set_values(
  clay_cropped,
  values = ifelse(
    st_intersects(st_as_sf(clay_cropped), thessaly_smus, sparse = FALSE),
    clay_cropped$Clay.pred,
    NA
  )
)

3. 导出为GIS兼容TIFF #

使用write_stars()导出，设置GIS可识别的无数据值与优化参数：

write_stars(
  clay_cropped,
  dsn = "clay_ok_predicted_cropped.tiff", # 输出文件路径
  driver = "GTiff",
  options = c("COMPRESS=LZW", "TILED=YES"), # 可选：压缩文件大小，提升GIS加载速度
  NA_value = NA # 标记无数据值，GIS软件可识别
)

内容的提问来源于stack exchange，提问作者Dimitris K