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

【问题标题】：R crop raster using polygon keeping cells along the borderR使用多边形裁剪栅格沿边界保持单元格
【发布时间】：2017-05-17 11:14:04
【问题描述】：

我正在尝试在意大利上空裁剪 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)

【问题讨论】：

标签： r polygon crop mask raster

【解决方案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)

【讨论】：

【解决方案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)

【讨论】：