使用多边形裁剪栅格数据，保留沿边界的单元。

Question

我想裁剪这个栅格，使其仅包含意大利内的区域。但是输出结果似乎丢失了沿边界的一些单元格。请参见下图中标记为红色的区域：

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

以下是我的脚本：

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)

- Claudia

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)

- jbaums

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)

- fdetsch