基于最大重叠，将多边形转换为栅格图像（使用 R 语言的 terra 或 stars 包）

Question

我有一个关于最大重叠多边形光栅化的问题，即将具有最高面积重叠的多边形的值分配给光栅单元。

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

问题在于terra包的函数（和类似stars的）将单元格值从包含单元格中心的多边形分配。如果一个光栅单元格包含多个多边形，则我更愿意选择光栅单元格中覆盖面积最大的多边形（土壤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个单元格），情况并非如此。 随着单元格与多边形相比越来越大，问题似乎会变得更加严重。因此，我可以从高分辨率栅格图开始，然后使用聚合函数（例如mode）重新采样到所需的（较低）分辨率。但是，也许有人有更有效的想法吗？

感谢您的帮助！

- paulsw

Answer 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那样优雅）

- Robert Hijmans

---

非常感谢你，Robert！你的第一个解决方案符合我的初始想法（高分辨率光栅化和聚合），对于更大的目标光栅来说可能是更好的方法。我也会尝试你的第二个想法。 - paulsw

非常感谢您分享使用 terra 的所有不同可能性。这真的非常启发人。干杯。 - lovalery

Answer 2

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

将SpatRaster r转换为SpatVector，然后再转换为sf对象，以利用sf::st_join()函数的largest=TRUE参数。其余代码则是简单地使用terra::rasterize()函数将sf对象再次转换为SpatVector，然后转换为SpatRaster。

因此，请查看下面详细介绍该过程的示例。 示例代码

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")

• 输出结果

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

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 (v2.0.1) 于2022-02-10创建的。}

---

使用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包（版本为2.0.1）于2022-02-10创建的。}

- lovalery

---

这个方案完美地运作，非常优雅！抱歉我的代码和图形不匹配，我稍微调整了一下分辨率，然后复制了错误的图形。 - paulsw

@paulsw，非常感谢您的反馈。很高兴我能帮到您。祝您工作顺利。干杯！ - lovalery

Answer 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)

- dbaston