循环播放剪辑 tiff 文件并导出到新文件夹

Question

我正在使用Bivand 的书来尝试弄清楚如何将 R 用于空间数据，但遇到了困难。

我有一个文件夹，其中包含从 1950 年到 2010 年的年度 tiff 文件（名称：tas_mean_C_cru_TS31_01_1950，列 4762，行 2557，源类型连续，像素类型浮点，像素深度 32 位，压缩 LZW）。

我想用多边形图层对它们进行剪辑，然后将剪辑后的栅格放入与原始栅格同名的新文件夹中。我有下面的代码，我试图从书中和在线的示例中进行操作，但它不起作用。我的问题出现在循环部分，最后出现错误。

我的两个问题是： 1）我的循环语句有什么问题？我正在尝试遍历 input.path 中定义的文件夹中的文件。 2）如果我退后一步，忽略剪辑，我怎样才能将文件复制并放置在同名的新文件夹中？

谢谢珍

---

#load libraries
library (sp)
library (rgdal)
library (raster)
library (maps)
library (mapproj)

#SET THE WORKING DIRECTORY PATH
input.path  = "F:/Dropbox/GIS/SNAP/Temp/tas50_09/"
#SET THE OUTPUT DIRECTORY
out.file<-"F:/temp/"
#FIND THE SIZE OF THE MATRIX -- WILL NEED THESE TO PRE-ALLOCATE SUBSEQUENT MATRICES
nrow <- 2557
ncol <- 4762
# initial rasters are....
# xmn=-2173225.11814296,xmx= 1498276.88185705, ymn = 409671.150470569, 
# ymx = 2381118.15047057, crs = '+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs')
# polygon to clip rasters
poly<-readOGR("F:/Dropbox/GIS/AKCommunities","AleutianComm50mile")
#reads the path to the files with a tif extention
file.names <- dir(input.path, pattern =".tif")
# then a loop for reading each existing file of type ".tif" as table
for(j in 1:length(input.path)){
  v <- extract(input.path, poly, weights=TRUE, cellnumbers=TRUE)
#Write the data to a geotiff
out <- raster(v ,xmn=-1446037.8545397,xmx= -686037.854539692, ymn = 356262.979377343, ymx = 607262.979377344, crs = '+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs')
writeRaster(v,filename=paste(out.file,'file.names ',sep=''),format = 'GTiff',options='COMPRESS=LZW',datatype='FLT4S',overwrite=T)
}
# I get the following error
Error in (function (classes, fdef, mtable)  : 
unable to find an inherited method for function ‘extract’ for signature ‘"character", "SpatialPolygonsDataFrame"’

Answer 1

以下编辑可能会帮助您发现逻辑中的问题。

当您需要迭代时，一致的命名方案是一件好事：

ipath <- "F:/Dropbox/GIS/SNAP/Temp/tas50_09/"
opath <- "F:/temp/"

ppath <- "F:/Dropbox/GIS/AKCommunities"
pname <- "AleutianComm50mile"

避免使用可能有命名空间冲突的对象名称，包括“poly”、“nrow”和“ncol”。

ppoly <- readOGR(ppath,pname)

请注意，您不会在任何地方使用这些值：

m <- 2557
n <- 4762

从循环中提取复杂的文字值可以提高可读性：

p1 <- '+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154'
p2 <- '+x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'
projargs <- paste(p1,p2)

## output raster limits
xmn <- -1446037.8545397
xmx <- -686037.854539692
ymn <- 356262.979377343
ymx <- 607262.979377344

使用 $ 锚定模式匹配字符串结尾并避免无意的不匹配。

files <- list.files(ipath, pattern='[.]tif$')
stopifnot(length(files)>0)

现在问题区域更容易阅读。

1. 您是否需要 extract() 中的权重和单元格编号似乎值得怀疑。
2. 如果 v 不是矩阵对象，您可能希望在 raster() 中使用 nrow=m 和 ncol=n 参数：

无论如何，现在应该很近了。

for (f in files) {

    ifile <- file.path(ipath,f)
    ofile <- file.path(opath,f)

    v  <- extract(ifile, ppoly, weights=TRUE, cellnumbers=TRUE)

    rv <- raster(v, crs=CRS(projargs), xmn=xmn, xmx=xmx, ymn=ymn, ymx=ymx)

    writeRaster(rv, filename = ofile, format='GTiff',
                options='COMPRESS=LZW', datatype='FLT4S', overwrite=TRUE)
}

Answer 2

我意识到这已经快两年了，但我想我可以将我的解决方案提供给遇到这个问题的任何人。我很高兴我创建了我的第一个 for 循环。我也有同样的错误信息。使用@goangit 的命名约定，我自己创建了一些东西。只是略有不同。大会非常有条理，所以我同意@goangit 并建议这成为未来的一种习惯。它是最短的小东西。不管怎样，下面是一个通用光栅裁剪循环的脚本。

#Set your working directory
setwd("C:/Here_it_is/")
ipath <- "C:/Here_it_is/"   
#I know it repeats the above, but I noticed for me, it needed the wd set
opath <- "C:/Where_its_going/"

#Identify your polygon boundary
ppath <- "C:/GIS/" #directory where the .shp is 
pname <- "boundary" #filename of the actual .shp without the extension.

#Read your polygon boundary
ppoly <- readOGR(ppath,pname) #if this doesn't work theyn maybe you do not have the package loaded? 

#Identify spatial extent
e <- extent(ppoly)

#Projection type and datum, respectively
p1 <- "+proj=UTM"
p2 <- "+datum=WGS84"
projmap <- paste(p1,p2)

#Identify the list of images you wish to perform the loop
files <- list.files(ipath, pattern=".tif$", full.names=FALSE)
stopifnot(length(files)>0)

#add output directory
outfiles <- pasteO(opath, files) #paste0 forgoes any separator between the two vectors "opath" and "files". if you need to add a slash somewhere then you need to change it to the paste function or change your vectors above.

#Change extension
extension(outfiles) <- "tif" #if you want to make it .tif 

for (f in 1:length(files)){
r <- raster(files[f])
rc <- crop(r, e)
rm <- mask(rc, ppoly)
rw <- writeRaster(rm,outfiles[f],overwrite=TRUE)
}

#这应该保存你指定outfiles向量的位置。