如何基于Shapefile中的单个多边形裁剪.tif卫星影像？

核心问题定位 #

你叠加失败的核心原因是卫星影像与地块矢量的投影坐标系（CRS）不一致，这是空间数据叠加操作失败的最常见诱因。以下方案基于当前R空间数据处理的主流工具terra（替代已逐步弃用的raster/rgdal包，性能更优、API更统一）实现。

步骤1：数据读取与投影统一 #

# 安装并加载必要工具包
install.packages(c("terra", "ggplot2"))
library(terra)
library(ggplot2)

# 读取卫星影像（转为terra标准栅格格式SpatRaster）
sat_img <- rast("你的卫星影像文件路径/2022_079000_449000_RGB_hrl.tif")

# 读取地块Shapefile（转为terra标准矢量格式SpatVector）
parcels <- vect("你的Shapefile文件夹路径/Essellanden_percelen.shp")

# 检查并统一投影
cat("影像CRS：", crs(sat_img), "\n")
cat("地块CRS：", crs(parcels), "\n")

if (!crs(sat_img) == crs(parcels)) {
  parcels <- project(parcels, crs(sat_img))
}

步骤2：叠加可视化（验证对齐效果） #

方法1：原生快速绘图（适配需求场景） #

# 绘制卫星影像（RGB对应影像3个波段）
plot(sat_img, rgb=1:3)
# 叠加红色地块边界
lines(parcels, col="red", lwd=1)

执行后即可得到你需要的叠加效果（与示例图一致）。

方法2：ggplot2精细化绘图（适合导出高清图） #

# 转换为ggplot兼容格式
img_df <- as.data.frame(sat_img, xy=TRUE)
colnames(img_df) <- c("x", "y", "R", "G", "B")

# 绘制叠加图
ggplot() +
  geom_raster(data=img_df, aes(x=x, y=y, fill=rgb(R/255, G/255, B/255))) +
  geom_sf(data=st_as_sf(parcels), color="red", fill=NA, lwd=0.8) +
  scale_fill_identity() +
  coord_sf(crs=crs(sat_img)) +
  theme_void()

步骤3：批量裁剪每个地块并导出为JPG/PNG #

# 创建存储导出文件的文件夹
dir.create("parcel_images", showWarnings=FALSE)

# 循环处理每个地块
for (i in 1:nrow(parcels)) {
  # 提取单个地块
  single_parcel <- parcels[i, ]
  
  # 裁剪+掩膜：保留地块内的影像区域
  clipped_img <- crop(sat_img, single_parcel)
  masked_img <- mask(clipped_img, single_parcel)
  
  # 导出为JPG（可替换为地块唯一ID作为文件名，需根据你的矢量属性列调整）
  output_path <- paste0("parcel_images/parcel_", i, ".jpg")
  writeRaster(masked_img, output_path, format="JPEG", overwrite=TRUE)
  
  # 如需导出PNG，修改format参数即可
  # writeRaster(masked_img, paste0("parcel_images/parcel_", i, ".png"), format="PNG", overwrite=TRUE)
}

注意事项 #

• 若需统一导出图像的尺寸（适配图像识别的标准化需求），可对裁剪后的影像用terra::resample重采样到固定分辨率，或用ggsave指定输出宽度/高度。
• 确保地块矢量包含唯一标识列（如ID），可将其作为文件名前缀，方便后续关联识别结果。

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