计算 R 中栅格单元内的多边形

Question

我正在尝试做一些我认为很简单的事情，即计算每个栅格单元中的多边形数量。我在这里尝试两种解决方案：R：如何计算网格的每个单元格中有多少个点？但它们都不起作用，而且我对使用的所有不同包感到疯狂。我想要一个

terra

解决方案，但不确定如何实现这一点，目前我对任何软件包持开放态度。 到目前为止，这是一些虚拟数据和我的代码：

# Load packages
pacman::p_load(terra, raster, rnaturalearth, rts, sf)

# Load dummy polygons
cntrs <- ne_countries(country = c("united kingdom", "argentina", "singapore"), type = "countries")

# Reproject
cntrs <- project(vect(cntrs), "+proj=moll +lon_0=0 +x_0=0 +y_0=0")

# Create empty raster grid
r <- terra::rast(xmin = -180, xmax = 180, ymin = -90, ymax = 90,
          res = 10,
          crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

r[] <- 0

### First solution of the link: rasterize
terra::rasterize(cntrs, r, fun = "count")

# gets a warning message and an empty raster
Warning message:
  Failed to compute min/max, no valid pixels found in sampling. (GDAL error 1) 

我不知道为什么光栅是空的。在链接的示例中，它们也有一些不与任何栅格单元重叠的点。

### Second solution: cellFromXY
# change polygons format
cntrs <- as(st_geometry(st_as_sf(cntrs)), "Spatial")
table(cellFromXY(r, cntrs))

# Gets error message:
Error in (function (classes, fdef, mtable)  : 
            unable to find an inherited method for function ‘cellFromXY’ for signature ‘"SpatRaster", "SpatialPolygons"’

即使我的

rts

版本是

rts_1.1-14

，与手册上相同（https://cran.r-project.org/web/packages/rts/rts.pdf），他们在其中写道“对象”可以是

SpatRaster

...

# Tried converting SpatRaster to RasterLayer as in the examples linked
r <- raster::raster(r)
table(cellFromXY(r, cntrs))

# Gets error message:
Error: Not compatible with requested type: [type=S4; target=double].

在示例中，他们使用

SpatialPoints

但我有多边形，所以我尝试使用

SpatialPolygons

，这似乎也不起作用。任何见解表示赞赏。

Answer 1

这是一种让你的代码正常工作的方法。也就是说，用匹配的坐标和 crs 定义栅格，然后对其进行转换。

library(rnaturalearth)
library(terra)
cntrs <- ne_countries(country = c("united kingdom", "argentina", "singapore"), type = "countries")

mollw <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0"
cntrs <- project(vect(cntrs), mollw)

r <- terra::rast(xmin = -180, xmax = 180, ymin = -90, ymax = 90, res = 10, crs = "+proj=longlat")
r <- project(r, mollw)

x <- terra::rasterize(cntrs, r)

rasterize

非常适合计数。但既然你正在处理多边形，我认为你追求的是

rasterGeom

。例如：

f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
r <- rast(v, res=.05)
x <- rasterizeGeom(v, r, "count")