使用 Polygon 裁剪光栅文件并使用相同的文件名写入输出

Posted 2023-03-23

【中文标题】使用 Polygon 裁剪光栅文件并使用相同的文件名写入输出

【英文标题】：Crop raster files with Polygon and writing the output with the same filename

【发布时间】：2019-05-04 16:03:36

【问题描述】：

我有很多光栅文件（卫星图像，所有可用 geotiff .tif 扩展名）。有些文件作为单个文件拆分为所有波段，有些文件有多个波段。由于这在我的硬盘驱动器上占用了大量空间，我想用我感兴趣的区域裁剪每个文件，我将其作为 shapefile 多边形。

我接近我自己的解决方案，并使用以下代码将裁剪后的图像作为新的 .tif 文件：

library(raster)
rasterfiles = list.files(path=getwd(), pattern = "*.TIF", full.names=TRUE)
s = stack(rasterfiles)
shp = readOGR("Area.shp")
rasterfiles_crop = crop(s, extent(shp))

output = writeRaster(rc, 'out.tif', format="GTiff", overwrite=TRUE, bylayer = TRUE)

使用此代码，我收到文件名 out_1.tif、out_2.tif 等...

不幸的是，生成的文件只有 1 个波段，因此 R 仅在涉及多波段 TIF 图像时识别第一个波段。

我想保留所有波段和原始文件名，并在新的末尾添加“_crop”。有人可以在这里帮助我如何更改代码吗？

谢谢

【参考方案1】：

你可以循环编写它们

library(raster)
library(rgdal)

shp <- readOGR("Area.shp")     

infiles <- list.files(path=getwd(), 
                      pattern="*.TIF", 
                      full.names=TRUE)

outfiles <- file.path(YourOutputPath, 
                      paste0(basename(tools::file_path_sans_ext(infiles)),
                      "_crop.tif")
                      )

for (i in seq_along(infiles)) 
  r <- crop(raster(infiles[i]), shp)
  writeRaster(r, filename=outfiles[i])

【讨论】：

谢谢，这段代码完美运行！ readOGR 是 rgdal 包中的一个函数，也必须加载。

我意识到，只有当原始文件在 TIFF 图像中有 1 个波段时，此方法才有效。或者换句话说：otherbadns 不在结果文件中。我不能使用类似 stack(list.files(path.....) 任何提示如何保留所有波段？

哦，你是对的！但是使用堆栈而不是光栅它应该可以工作。如果在循环中使用堆栈，会出现错误吗？

【参考方案2】：

我现在找到了一个解决方案，以下代码列出了文件夹中的所有 TIF 文件，并且多波段 tif 在裁剪过程后保留其波段：

library(raster)
library(rgdal)

setwd("input-folder")

 ## polygon with crop-extend ##
shp <- readOGR("area.shp")

 ## load tif files ##
infiles = list.files(path=getwd(),
                     pattern="*.tif$|*.TIF$")

 ## Filenames with desired suffix and output place ##
outfiles = file.path("D:/Downloads/BDA/Output",
                     paste0(basename(tools::file_path_sans_ext(infiles)),
                     ".tif"))

 ## crop and output settings (compression and datatype)
 for (i in seq_along(infiles)) 
      r = crop(stack(infiles[i]), shp)
      writeRaster(r, filename=outfiles[i],
                  bylayer=FALSE,
                  format="GTiff",
                  datatype="INT1U",
                  options="COMPRESS=ZIP",
                  overwrite=TRUE)
 

感谢 Richard 提供的精彩循环代码！

关于数据类型：如果 R 可以检查输入文件的数据类型并自动为裁剪输出选择相同的数据类型，那就太好了。现在我必须手动指定数据类型。否则，即使输入文件只有 8 位无符号 (INT1U) 或 16 位有符号 (INT2S)，输出文件也会为 float32 (FLT4S)。

 datatype= same.as.input.file