2013-07-12 23 views
1

我與空間數據的準備去分析裁剪柵格圖層中的R - 我在我的研究領域所期望的程度DEM,雖然我在〜39其它層全國範圍(美國)。有沒有辦法將所有這39層與DEM同時進行裁剪?你怎麼能在一個批處理和改變投影

另外,我將在覆蓋不同的投影與其它層的輸出。是否可以調整輸出圖層的投影和像素大小?

我想使用免費軟件儘可能地爲我的數據操作...

+0

嘗試使用package sp或rgdal更改所有圖層的投影。 –

+1

謝謝!我想出了一個解決方案,並希望得到它貼在這裏,以便其他人可以使用它,但因爲我是一個新的用戶我回答之前等待8小時:-)。我結束了使用軟件包柵格(也調用sp和gdal)來完成這個... – mtreg

回答

3

我有上面的問題,但已經寫在R A函數來做到這一切在一個批次 - 見下文。我有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, [email protected] ### 
###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)