r - 如何在 R 中批量裁剪栅格图层并更改投影

标签 r gis raster map-projections

我正在处理空间数据,为分析做好准备 - 我在我的研究区域的所需范围内有一个 DEM,尽管我在国家范围(美国)还有大约 39 个其他图层。有没有办法一次性将所有这 39 个图层裁剪到与 DEM 相同的程度?

此外,我将在不同的投影中将输出与其他图层叠加。是否可以调整输出层的投影和像素大小?

我正在尝试尽可能多地使用免费软件来进行数据操作...

最佳答案

我遇到了上面的问题,但在 R 中编写了一个函数来批量执行所有这些操作 - 见下文。我有美国大陆规模的 39 个气候数据图层(来自 PRISM 气候数据组; http://www.prism.oregonstate.edu/ ),并且希望将它们裁剪到加利福尼亚州南部 DEM 的范围,重新投影它们,并将它们导出以便于导入并与 SAGA GIS 中的其他图层一起使用。下面是代码,以及如何运行它的示例,将工作目录设置为包含要裁剪的图层的文件夹,并且仅包含这些图层。

在处理过程中,所有数据都存储在内存中,因此对于庞大的数据集,它可能会因为内存不足而挂起......这可能是一个值得改进的地方。

此外,R 论坛上的回复也提供了一种更短、更优雅的方法:http://permalink.gmane.org/gmane.comp.lang.r.geo/18320

我希望有人觉得它有用!

#########################################
#BatchCrop Function ###
#by Mike Treglia, <a href="https://stackoverflow.com/cdn-cgi/l/email-protection" class="__cf_email__" data-cfemail="f29f868097959e9b93b2959f939b9edc919d9f" rel="noreferrer noopener nofollow">[email protected]</a> ###
###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 - 如何在 R 中批量裁剪栅格图层并更改投影,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/17624977/

相关文章:

r - 通过使用多种颜色填充渐变来平铺密度图

可以处理重复的不规则时间序列的滚动窗口函数

r - 使用带有简单特征 (SF) 函数的 apply()

python - 如何使用 python 查找 ESRI shapefile 中特定纬度/经度的信息?

java - 免费GIS服务器获取地理数据

r - 如何在 R 中使用 for 循环保存具有不同名称的文件?

r - R中使用marrangeGrob和ggsave的页码位置

gis - 将栅格数据导入 NetLogo 会生成一列所有补丁变量 = 0

r - 在 R 中将路径/路线图编写为 GeoTiff

r - 将单个栅格设置为 NA,其中栅格堆栈的值为 NA