通过在R中填充NA,使光栅达到相同的程度



我有几个具有不同几何体/轮廓的裁剪光栅。具体来说,同一块田地几年的空间产量图,但程度各不相同——测量结果并不总是整个田地的整体,但在某些年份只是其中的一部分。我想计算这些地图的平均值,并将它们组合成一个平均值栅格。然而,这意味着,并不是说5层/光栅中的每个像素都有一个值。我可以接受这些缺失的值是NA,所以最终的平均值只能通过假设字段部分的3个光栅来计算,其中映射不重叠。

我想用"extend{raster}"扩展光栅,用NA值填充非重叠部分:

y <- extend(y, shape, value=NA)#Shape是一个矩形,它分解了所有的成品图光栅

这对所有光栅都很有效。但他们仍然没有相同的程度。即使我通过setExtent()extent() <- extent()将范围调整为矩形文件的范围,甚至调整为其他扩展光栅之一,我仍然得到:

compareMaster(x)中的错误:不同的数量或列

当我想堆叠它们并使用calc(y, fun=mean,...)时。原始光栅范围差异太大,无法重新采样。但它们确实有相同的分辨率和CRS。

有人知道如何解决这个问题吗?

如果您想实现calc(y, fun=mean,...)这样的操作,您可以获得光栅的最小公共范围,并在将它们堆叠在一起并应用该操作之前将它们全部裁剪到该范围。

假设你有三个光栅:

# Generate 3 dummy rasters with different extents
r1 <- raster( crs="+proj=utm +zone=31")
extent(r1) <- extent(0, 100, 0, 500)
res(r1) <- c(5, 5)
values(r1) <- sample(10, ncell(r1), replace=TRUE)
r2 <- raster( crs="+proj=utm +zone=31")
extent(r2) <- extent(10, 120, -10, 400) 
res(r2) <- c(5, 5)
values(r2) <- runif(ncell(r2), 1, 10)
r3 <- raster( crs="+proj=utm +zone=31")
extent(r3) <- extent(50, 150, 30, 200) 
res(r3) <- c(5, 5)
values(r3) <- runif(ncell(r3), 1, 10)

第一种方法:

# Summing your rasters will only work where they are not NA
r123 = r1+r2+r3 # r123 has the minimal common extent
r1 = crop(r1, r123) # crop to that minimal extent
r2 = crop(r2, r123)
r3 = crop(r3, r123)
s123 = stack(r1, r2, r3)
s123.mean = calc(s123, fun=mean)

另一个:

# Manually compute the minimal extent
xmin <- max(bbox(r1)[1,1], bbox(r2)[1,1], bbox(r3)[1,1])
xmax <- min(bbox(r1)[1,2], bbox(r2)[1,2], bbox(r3)[1,2])  
ymin <- max(bbox(r1)[2,1], bbox(r2)[2,1], bbox(r3)[2,1])  
ymax <- min(bbox(r1)[2,2], bbox(r2)[2,2], bbox(r3)[2,2])  
newextent=c(xmin, xmax, ymin, ymax)
r1 = crop(r1, newextent)
r2 = crop(r2, newextent)
r3 = crop(r3, newextent)
s123 = stack(r1, r2, r3)
s123.mean = calc(s123, fun=mean)

我不太确定你到底想要什么,但如果你真的想保持所有光栅的完整性(没有crop操作),你也可以计算最大的范围(与上面的第二种方法相同,只需将minmax反转)和extend所有你的光栅在堆叠它们之前都是那个(同样的事情,你最终得到的只是更大的光栅,填充了NA……不确定这真的是想要的):

xmin <- min(bbox(r1)[1,1], bbox(r2)[1,1], bbox(r3)[1,1])
xmax <- max(bbox(r1)[1,2], bbox(r2)[1,2], bbox(r3)[1,2])  
ymin <- min(bbox(r1)[2,1], bbox(r2)[2,1], bbox(r3)[2,1])  
ymax <- max(bbox(r1)[2,2], bbox(r2)[2,2], bbox(r3)[2,2])  
newextent=c(xmin, xmax, ymin, ymax)
r1 = extend(r1, newextent)
r2 = extend(r2, newextent)
r3 = extend(r3, newextent)
s123 = stack(r1, r2, r3)
s123.mean = calc(s123, fun=mean)

这有帮助吗?

注意

你可能不想玩setExtent()extent() <- extent(),因为你可能会以光栅的错误地理坐标结束(即,范围会被修改,但内容不会被修改,所以你实际上可能会翻译光栅,而不是针对它,因为它们之间的重叠在物理上毫无意义)。

我也遇到了同样的问题,要破坏一些光栅,这些光栅包含bioclim数据。结果发现,我确实正确裁剪了,但下载的分辨率不同。查看每个光栅的列数。只需在R中键入光栅(文件)的名称即可完成。它们可能大小相同。我希望我能帮忙。

最新更新