根据最大重叠栅格化多边形（使用 R 包 terra 或 stars）答案

【问题标题】：Rasterize polygons based on maximum overlap (using R packages terra or stars)根据最大重叠栅格化多边形（使用 R 包 terra 或 stars）
【发布时间】：2022-06-28 23:10:07
【问题描述】：

我有一个关于通过最大重叠对多边形进行栅格化的问题，即分配与栅格单元重叠面积最大的多边形的值。

现实世界的练习是在 R 中栅格化土壤 ID 的多边形，以便生成分辨率相对较低的土壤特性图作为模型输入。

问题在于 terra 包的 rasterize() 函数（以及类似的星星的 st_rasterize()）从包含单元中点的多边形中分配单元值。如果一个栅格单元包含多个多边形，我宁愿选择多边形的值（soil-ID），它在栅格单元中具有最高的面积覆盖。

这是一个使用 terra 可视化我的问题的小型独立示例。

library(terra)

f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
r <- rast(v, ncols = 3, nrow = 3)
rcc <- vect(xyFromCell(r, cell = 1:ncell(r)))

x <- rasterize(v, r, field = "NAME_2")
plot(x)
lines(r, col = "light gray")
lines(v)
points(rcc)

大多数情况下，包含单元中心的多边形似乎也具有最高的面积份额。但是，在某些情况下（顶行，第 3 个单元格），情况并非如此。与多边形相比，单元格越大，问题似乎越严重。因此，我可以从高分辨率栅格开始，然后使用聚合函数（例如模式）重新采样到所需的（较低）分辨率。但是，也许有人有更有效的想法？

感谢您的帮助！

【问题讨论】：

标签： r sf terra

【解决方案1】：

你可以这样做：

#Example data
library(terra)    
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
r <- rast(v, ncols = 3, nrow = 3)

n <- 10
r <- disagg(r, n)
r <- rasterize(v, r, "ID_2")
x <- aggregate(r, n, "modal")

plot(x)
lines(x)
lines(v, lwd=2)
text(v, col="red", halo=T)
text(x, col="blue", halo=T)

另一种方式，可能效率较低（特别是如果您有很多 ID）：

z <- lapply(1:nrow(v), \(i) rasterize(v[i,], r, cover=TRUE))
z <- which.max(rast(z))

但是如果你想要非常高的精度，你可以用 exactextractr::coverage_fraction 替换 rasterize

我想效率甚至更低：

r <- rast(v, ncols = 3, nrow = 3)
values(r) <- 1:ncell(r)
# get weights
e <- extract(r, v, weights=TRUE)
e <- as.matrix(e)
head(e)
#    ID lyr.1 weight
#[1,]  1     1   0.38
#[2,]  1     2   0.49
#[3,]  2     2   0.06
#[4,]  2     4   0.05
#[5,]  2     5   0.52
#[6,]  2     6   0.06

# find cell with max weight (you can use dplyr or data.table intead) 
x <- sapply(unique(e[,2]), function(i) { 
    d <- e[e[,2] == i, ,drop=FALSE]
    d[which.max(d[,3]), 2:1]
})

# remove values 
r <- rast(r)
# assign ID to cells
r[x[1,]] <- x[2,]

您可以使用多边形交集来实现相同的效果，但这不能很好地扩展到大型栅格

r <- rast(v, ncols = 3, nrow = 3)
values(r) <- 1:9
v$ID <- 1:nrow(v)
i <- intersect(v[,"ID"], as.polygons(r))
i$area <- expanse(i)
i <- data.frame(i) 
x <- sapply(split(i, i[,2]), 
    \(x) { x[which.max(x[,3]), 2:1] |> unlist()}
)
r <- rast(r)
r[x[1,]] <- x[2,]

（或许没有lovalery提出的st_join那么优雅）

【讨论】：

非常感谢，罗伯特！您的第一个解决方案提出了我的初始想法（以高分辨率和聚合进行栅格化），并且可能是较大目标栅格的更好方法。我也会尝试你的第二个想法。
非常感谢您仅使用terra 分享所有不同的可能性。这真的很有启发性。干杯。

【解决方案2】：

请使用terra 和sf 库找到一种可能的解决方案。

这个想法是将SpatRaster r 转换为SpatVector，然后转换为sf 对象，以便使用sf::st_join() 函数使用largest = TRUE 参数。剩下的代码包括简单地将sf 对象转换回SpatVector，然后使用terra::rasterize() 函数转换为SpatRaster。

所以，请在下面找到详细说明该过程的代表。

Reprex

代码

library(terra)
library(sf)

# Your data
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
r <- rast(v, ncols = 3, nrow = 3)
rcc <- vect(xyFromCell(r, cell = 1:ncell(r)))

# Convert the 'SpatRaster' 'r' into a 'SpatVector (i.e. 'r_poly')
r_poly <- terra::as.polygons(r)

# Convert 'r_poly' into a 'sf' object (i.e. 'r_poly_sf')
r_poly_sf <- sf::st_as_sf(r_poly)

# Convert 'v' into a 'sf' object (i.e. 'v_sf')
v_sf <- sf::st_as_sf(v)

# Left join r_poly_sf with v_sf based on the largest overlap
results_sf <- sf::st_join(r_poly_sf, v_sf, largest = TRUE)

# Convert 'results_sf' into a SpatVector (i.e. 'results_vect')
results_vect <- terra::vect(results_sf)

# Rasterize 'results_vect' to get a 'SpatRaster' (i.e. 'results')  
results <- terra::rasterize(results_vect, r, field = "NAME_2")

输出

注意：请注意右上角的单元格是NA，因为来自r 的多边形没有与v 重叠（如果需要，您仍然可以使用@987654341 为不重叠的单元格设置值terra::rasterize() 函数中的 @ 参数）。

results
#> class       : SpatRaster 
#> dimensions  : 3, 3, 1  (nrow, ncol, nlyr)
#> resolution  : 0.2613707, 0.2446047  (x, y)
#> extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : memory 
#> name        :   NAME_2 
#> min value   : Capellen 
#> max value   :   Remich

terra::values(results, dataframe=TRUE)
#>       NAME_2
#> 1   Clervaux
#> 2   Clervaux
#> 3       <NA>
#> 4    Redange
#> 5     Mersch
#> 6 Echternach
#> 7   Capellen
#> 8 Luxembourg
#> 9     Remich

可视化

plot(results)
lines(r, col = "light gray")
lines(v)
points(rcc)

^{由reprex package 创建于 2022-02-10 (v2.0.1)}

与上面完全相同，但使用rast(v, ncols = 5, nrow =5) 可以获得可以与您在问题中提供的数字进行比较的结果。

library(terra)
library(sf)

# Your data
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
r <- rast(v, ncols = 5, nrow = 5)
rcc <- vect(xyFromCell(r, cell = 1:ncell(r)))

# Convert the 'SpatRaster' 'r' into a 'SpatVector (i.e. 'r_poly')
r_poly <- terra::as.polygons(r)

# Convert 'r_poly' into a 'sf' object (i.e. 'r_poly_sf')
r_poly_sf <- sf::st_as_sf(r_poly)

# Convert 'v' into a 'sf' object (i.e. 'v_sf')
v_sf <- sf::st_as_sf(v)

# Left join r_poly_sf with v_sf based on the largest overlap
results_sf <- sf::st_join(r_poly_sf, v_sf, largest = TRUE)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

# Convert 'results_sf' into a SpatVector (i.e. 'results_vect')
results_vect <- terra::vect(results_sf)

# Rasterize 'results_vect' to get a 'SpatRaster' (i.e. 'results')  
results <- terra::rasterize(results_vect, r, field = "NAME_2")
results
#> class       : SpatRaster 
#> dimensions  : 5, 5, 1  (nrow, ncol, nlyr)
#> resolution  : 0.1568224, 0.1467628  (x, y)
#> extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : memory 
#> name        :   NAME_2 
#> min value   : Capellen 
#> max value   :    Wiltz

terra::values(results, dataframe=TRUE)
#>              NAME_2
#> 1          Clervaux
#> 2          Clervaux
#> 3          Clervaux
#> 4              <NA>
#> 5              <NA>
#> 6             Wiltz
#> 7             Wiltz
#> 8           Vianden
#> 9           Vianden
#> 10             <NA>
#> 11          Redange
#> 12          Redange
#> 13           Mersch
#> 14       Echternach
#> 15       Echternach
#> 16         Capellen
#> 17         Capellen
#> 18       Luxembourg
#> 19     Grevenmacher
#> 20     Grevenmacher
#> 21 Esch-sur-Alzette
#> 22 Esch-sur-Alzette
#> 23 Esch-sur-Alzette
#> 24           Remich
#> 25           Remich

plot(results)
lines(r, col = "light gray")
lines(v)
points(rcc)

^{由reprex package (v2.0.1) 于 2022-02-10 创建}

【讨论】：

这很完美，是一个非常优雅的解决方案！对不起，我的代码和图形不匹配，我尝试了一下分辨率，然后复制了错误的图形。
@paulsw，非常感谢您的反馈。很高兴我能帮上忙。我祝你工作顺利。干杯。

【解决方案3】：

使用exactextractr >= 0.9（目前在 GitHub 上），这可以通过

library(terra)
library(exactextractr)
library(sf)

v <- st_read(system.file("ex/lux.shp", package="terra"))
r <- rast(ext(v), ncols = 3, nrow = 3)
r <- rasterize_polygons(v, r)

plot(r)
plot(st_geometry(v), add=TRUE)

【讨论】：