如何批量裁剪R中的栅格图层并更改投影 [英] How can you crop raster layers in R in a batch and change projection
问题描述
我正在处理空间数据，以准备进行分析-尽管在我的国家范围内(美国)还有39个其他层次，但我在研究区域的期望范围内拥有DEM.有没有一种方法可以将这39个图层全部裁剪到与DEM相同的程度?
I was working with spatial data to get ready for analyses - I have a DEM at the desired extent of my study area, though I have ~39 other layers at the national scale (US). Is there a way to crop all of these 39 layers to the same extent as the DEM at once?
此外，我将在其他投影中将输出与其他图层叠加.是否可以调整输出层的投影和像素大小?
Also, I will be overlaying the output with other layers in a different projection. Is it possible to adjust the projection and pixel size of the output layers?
我正在尝试尽可能多地使用免费软件进行数据操作...
I am trying to use freeware as much as possible for my data manipulation...
推荐答案
我遇到了上述问题，但是已经在R中编写了一个函数来批量完成所有这些工作-参见下文.我拥有39个美国大陆范围的气候数据层(来自PRISM气候数据组； http://www.并希望将其裁剪到DEM的程度，在加利福尼亚州南部进行重新投影，然后将其导出以方便导入并与SAGA GIS中的其他图层一起使用.下面是代码，并提供了如何运行该代码的示例，将工作目录设置为包含要裁剪的图层且只有这些图层的文件夹.
I had the problem above, but have written a function in R to do all of this in a batch - see below. I had 39 climate data layers at the scale of the continental U.S. (from PRISM Climate Data group; http://www.prism.oregonstate.edu/), and wanted to clip them to the extent of a DEM, in southern California, reproject them, and export them for easy import and use with other layers in SAGA GIS. Below is the code, with an example of how it would be run, setting the working directory to the folder that has the layers that you want to crop, and only those layers.
在处理过程中，所有数据都存储在内存中，因此对于庞大的数据集，由于缺少内存而可能会被挂起……这可能会有所改善.
During the processing, all data are stored in memory, so with huge datasets, it might get hung up because of lack of memory... that would probably be something that would be good to improve.
此外，R论坛上的回复也提供了一种更简短，更优雅的方法: http://permalink.gmane.org/gmane.comp.lang.r.geo/18320
Also, a response on the R Forum provided a shorter, more elegant way to do it too: http://permalink.gmane.org/gmane.comp.lang.r.geo/18320
我希望有人觉得它有用！
I hope somebody finds it useful!
#########################################
#BatchCrop Function ###
#by Mike Treglia, mtreglia@gmail.com ###
###Tested in R Version 3.0.0 (64-bit), using 'raster' version 2.1-48 and 'rgdal' version 0.8-10
########################################
#This function crops .asc raster files in working directory to extent of another layer (referred to here as 'reference' layer), converts to desired projection, and saves as new .asc files in the working directory. It is important that the original raster files and the reference layer are all in the same projection, though different pixel sizes are OK. The function can easily be modified to use other raster formats as well
#Note, Requires package 'raster'
#Function Arguments:
#'Reference' refers to name of the layer with the desired extent; 'OutName' represents the intended prefix for output files; 'OutPrj' represents the desired output projection; and 'OutRes' represents the desired Output Resolution
BatchCrop<-function(Reference,OutName,OutPrj,OutRes){
filenames <- list.files(pattern="*.asc", full.names=TRUE) #Extract list of file names from working directory
library(raster) #Calls 'raster' library
#Function 'f1' imports data listed in 'filenames' and assigns projection
f1<-function(x,z) {
y <- raster(x)
projection(y) <- CRS(z)
return(y)
}
import <- lapply(filenames,f1,projection(Reference))
cropped <- lapply(import,crop,Reference) #Crop imported layers to reference layer, argument 'x'
#Function 'f2' changes projectection of cropped layers
f2<-function(x,y) {
x<-projectRaster(x, crs=OutPrj, res=OutRes)
return(x)
}
output <- lapply(cropped,f2,OutPrj)
#Use a 'for' loop to iterate writeRaster function for all cropped layers
for(i in (1:max(length(filenames)))){ #
writeRaster(output[[i]],paste(deparse(substitute(OutName)), i), format='ascii')
}
}
#############################################
###Example Code using function 'BatchCrop'###
#############################################
#Data layers to be cropped downloaded from: http://www.prism.oregonstate.edu/products/matrix.phtml?vartype=tmax&view=maps [testing was done using 1981-2010 monthly and annual normals; can use any .asc layer within the bounds of the PRISM data, with projection as lat/long and GRS80]
#Set Working Directory where data to be cropped are stored
setwd("D:/GIS/PRISM/1981-2010/TMin")
#Import Reference Layer
reference<-raster("D:/GIS/California/Hab Suitability Files/10m_DEM/10m DEM asc/DEM_10m.asc")
#Set Projection for Reference Layer
projection(reference) <- CRS("+proj=longlat +ellps=GRS80")
#Run Function [desired projection is UTM, zone 11, WGS84; desired output resolution is 800m]
BatchCrop(Reference=reference,OutName=TMinCrop,OutPrj="+proj=utm +zone=11 +datum=WGS84",OutRes=800)
这篇关于如何批量裁剪R中的栅格图层并更改投影的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
