根据最大重叠栅格化多边形（使用 R 包 terra 或 stars）
Rasterize polygons based on maximum overlap (using R packages terra or stars)
我有一个关于通过最大重叠对多边形进行栅格化的问题，即分配与栅格单元重叠面积最大的多边形的值。
现实世界的练习是在 R 中栅格化土壤 ID 的多边形，以便生成分辨率相对较低的土壤特性图作为模型输入。
问题是 terra 包的 rasterize() 函数（和类似的星星 st_rasterize()）从包含像元中点的多边形分配像元值。如果栅格单元格包含多个多边形，我宁愿 select 多边形 (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 个单元格），情况并非如此。
与多边形相比，单元格越大，问题似乎越严重。因此，我可以从高分辨率栅格开始，然后使用聚合函数（例如模式）重新采样到所需的（较低的）分辨率。但是，也许有人有更有效的想法？
感谢您的帮助！
请使用 terra 和 sf 库找到一种可能的解决方案。
想法是将 SpatRaster r 转换为 SpatVector，然后转换为 sf 对象，以便利用 sf::st_join()使用 largest = TRUE 参数的函数。其余代码包括简单地将 sf 对象转换回 SpatVector，然后使用 terra::rasterize() 函数将 SpatRaster 转换为 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 重叠（如果需要，您仍然可以设置单元格的值不要在 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 package (v2.0.1)
于 2022-02-10 创建
你可以这样做：
#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 替换光栅化
效率更低，我想：
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优雅）
我有一个关于通过最大重叠对多边形进行栅格化的问题，即分配与栅格单元重叠面积最大的多边形的值。
现实世界的练习是在 R 中栅格化土壤 ID 的多边形，以便生成分辨率相对较低的土壤特性图作为模型输入。
问题是 terra 包的 rasterize() 函数（和类似的星星 st_rasterize()）从包含像元中点的多边形分配像元值。如果栅格单元格包含多个多边形，我宁愿 select 多边形 (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 个单元格），情况并非如此。 与多边形相比，单元格越大，问题似乎越严重。因此，我可以从高分辨率栅格开始，然后使用聚合函数（例如模式）重新采样到所需的（较低的）分辨率。但是，也许有人有更有效的想法？
感谢您的帮助！
请使用 terra 和 sf 库找到一种可能的解决方案。
想法是将 SpatRaster r 转换为 SpatVector，然后转换为 sf 对象，以便利用 sf::st_join()使用 largest = TRUE 参数的函数。其余代码包括简单地将 sf 对象转换回 SpatVector，然后使用 terra::rasterize() 函数将 SpatRaster 转换为 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重叠（如果需要，您仍然可以设置单元格的值不要在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 package (v2.0.1)
于 2022-02-10 创建你可以这样做：
#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 替换光栅化
效率更低，我想：
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优雅）