R裁剪栅格使用多边形保持沿边界的单元格

Question

我正在尝试在意大利裁剪 this raster，但输出结果似乎遗漏了边界沿线的一些像元。请参见下图中以红色突出显示的区域：

如何保留所有越界的单元格？

下面是我的脚本：

library(raster)

# Load data
x <- raster("x.nc")
IT <- getData(name = "GADM", country = "Italy", level = 0)

# Mask and crop
x_masked <- mask(x, IT)
x_masked_cropped <- crop(x_masked, IT)

# Plot
plot(x_masked_cropped)
plot(IT, add = T)

Answer 1

这是一种方法。我们使用 gdalUtils::gdal_rasterize 创建一个二进制掩码栅格，使用 at=TRUE 确保值 1 被烧入意大利多边形接触的所有单元格。 gdal_rasterize 指的是磁盘上的文件，因此首先将 IT 写入 OGR 支持的文件。

library(gdalUtils)
library(rgdal)
x_crop <- crop(x, IT)
writeOGR(IT, tempdir(), f <- basename(tempfile()), 'ESRI Shapefile')
gdal_rasterize(sprintf('%s/%s.shp', tempdir(), f), 
               f2 <- tempfile(fileext='.tif'), at=T,
               tr=res(x_crop), te=c(bbox(x_crop)), burn=1, 
               init=0, a_nodata=0, ot='Byte')

plot(x_crop*raster(f2)) # multiply the raster by 1 or NA
plot(IT, add=TRUE)

Answer 2

您可以识别所有与意大利相交的栅格像元并将剩余的（即不相交的像素）设置为 NA。确保通过 cellFromPolygon(..., weights = TRUE) 检索具有各自权重的单元格 - 否则，只会返回中心位于意大利境内的单元格（另请参见 ?raster::extract）。

## identify cells covering italy and set all remaining pixels to NA
cls <- cellFromPolygon(x, IT, weights = TRUE)[[1]][, "cell"]
x[][-cls] <- NA

plot(trim(x))
plot(IT, add = TRUE)

R裁剪栅格使用多边形保持沿边界的单元格

R crop raster using polygon keeping cells along the border

r

mask

crop

polygon

raster