在 R 中创建非矩形栅格

编程问答 2022-05-13

问题描述

在我的分析中，我必须遍历栅格的单元格 ID（即单元格编号），这是非可选的，并且在我无法熟练掌握的不同软件设置 (BUGScode) 中完成。我的栅格必须完全适合 shapefile（我的研究区域），否则我将循环遍历我没有任何值或兴趣的单元格编号，然后......你知道会发生什么;）。

可重现的例子

### get a shapefile
GER<-getData("GADM",country="DEU",level=1,download=T,path=getwd())

### make a raster
raster<-raster(ext=extent(GER),vals=0,crs=crs(GER))

plot(raster);plot(GER,add=T)

### Now make the raster fit exactly to the German borders
mask.raster<-mask(raster,GER)

plot(mask.raster) ### looks right,but WAIT
plot(mask.raster,colNA="black") ### the border cells were only set to NA,which is problematic if your raster contains NA in the first place,i.e. if you used rasterize to count points in the raster cells

trim.raster<-trim(raster,values=NA)

plot(trim.raster,colNA="black"); plot(GER,add=T) ### trim only trims complete rows/columns from the raster

crop.raster<-crop(raster,GER)

plot(crop.raster,colNA="black");plot(GER,add=T) ### crop only crops to extent

如何使我的栅格与我的研究区域的轮廓完全吻合，从而不成为矩形？

这在各种 GIS 软件中定期完成，但我目前仅限于 R。

解决方法

正如我最近发现的那样，R 中的 raster 意味着某种投影矩阵，因此只能是矩形（或在特殊情况下为对角线），但实际上从未适合形状。

至于当您只是简单地屏蔽光栅文件时会发生什么。 如果您循环说 cell in 1:ncell(raster)，您可能会遇到以下问题：

  1. 处理时间更长，因为您必须遍历可能更多的单元格。
  2. 合并总和等时的错误消息。
  3. 就我而言，有一个问题，如果您使用 rasterize 计算栅格中的某些点，并且计数为 0（这被解释为 NA ... 无计可数）同一类单元格在您的学习区域内被跳过。通过将“计数”为空 = NA 的栅格值设置为 0，然后屏蔽栅格，可以轻松避免这种情况，但这是许多解决方法的一种情况...
,

根据大多数定义以及支持它们的文件格式，栅格是矩形的。您断言在各种 GIS 软件中都不是这种情况，但我认为这不是真的，而且您没有提供任何证据。也许您正在考虑（类似光栅的）矩形多边形？

另外，我不知道会发生什么。大多数有能力的 R 软件只会跳过 NA 的栅格单元。

rrr-rasterspatial