R中栅格化重叠sf对象：指定像素值与填充空缺区域

我有一个包含多个带关联数值多边形的sf数据框，多边形位置略有差异且存在重叠，需要基于另一sf多边形的范围，使用stars::st_rasterize完成栅格化，需解决以下两个问题：

• 当栅格像素内存在多个重叠多边形时，取这些对象中的最小数值；
• 通过最近邻插值填充多边形间的空隙，确保目标范围内的栅格无NA值。

示例数据代码 #

a) 作为栅格范围的多边形 #

library(tidyverse)
library(sf)
library(stars)

# 设置顶点
polygon_vertices <- matrix(c(-1489008, 1369848,
                             -1488191, 1370128,
                             -1488132, 1369347,
                             -1488985, 1369243,
                             -1489008, 1369848),
                           ncol = 2,
                           byrow = TRUE) 

# 创建多边形并转为sf对象
polygon_extent <- st_polygon(list(polygon_vertices)) %>%
  st_sfc(crs = 6350)

b) 包含重叠多边形的sf数据框 #

# 定义各多边形顶点
p1_vertices <- matrix(c(-1488355, 1369737,  
                        -1488355, 1370072,  
                        -1488237, 1370112, 
                        -1488190, 1370112,
                        -1488162, 1369737,
                        -1488355, 1369737),
                      ncol = 2, byrow = TRUE)
p2_vertices <- matrix(c(-1488355, 1369737,  
                        -1488355, 1370072,  
                        -1488237, 1370112, 
                        -1488190, 1370112,
                        -1488162, 1369737,
                        -1488355, 1369737),
                      ncol = 2, byrow = TRUE)
p3_vertices <- matrix(c(-1488253, 1369484,  
                        -1488253, 1369859,  
                        -1488171, 1369859, 
                        -1488143, 1369484,
                        -1488253, 1369484),
                      ncol = 2, byrow = TRUE) 
p4_vertices <- matrix(c(-1488292, 1369795,  
                        -1488292, 1370093,  
                        -1488191, 1370128, 
                        -1488166, 1369795,
                        -1488292, 1369795),
                      ncol = 2, byrow = TRUE) 
p5_vertices <- matrix(c(-1488260, 1369634,  
                        -1488154, 1369634,  
                        -1488132, 1369347, 
                        -1488260, 1369332,
                        -1488260, 1369634),
                      ncol = 2, byrow = TRUE) 
p6_vertices <- matrix(c(-1488255, 1369734,  
                        -1488630, 1369734,  
                        -1488630, 1369977, 
                        -1488255, 1370106,
                        -1488255, 1369734),
                      ncol = 2, byrow = TRUE) 
p7_vertices <- matrix(c(-1488240, 1369419,  
                        -1488615, 1369419,  
                        -1488615, 1369794, 
                        -1488240, 1369794,
                        -1488240, 1369419),
                      ncol = 2, byrow = TRUE) 
p8_vertices <- matrix(c(-1488212, 1370120,  
                        -1488191, 1370128,  
                        -1488190, 1370120, 
                        -1488212, 1370120),
                      ncol = 2, byrow = TRUE)
p9_vertices <- matrix(c(-1488970, 1369519,  
                        -1488995, 1369519,  
                        -1489008, 1369848, 
                        -1488970, 1369861,
                        -1488970, 1369519),
                      ncol = 2, byrow = TRUE) 
p10_vertices <- matrix(c(-1488970, 1369519,  
                        -1488995, 1369519,  
                        -1489008, 1369848, 
                        -1488970, 1369861,
                        -1488970, 1369519),
                      ncol = 2, byrow = TRUE) 
p11_vertices <- matrix(c(-1488518, 1369823,  
                         -1488518, 1370016,  
                         -1488191, 1370128, 
                         -1488168, 1369823,
                         -1488518, 1369823),
                       ncol = 2, byrow = TRUE)
p12_vertices <- matrix(c(-1488452, 1369986,  
                         -1488452, 1370039,  
                         -1488191, 1370128, 
                         -1488180, 1369986,
                         -1488452, 1369986),
                       ncol = 2, byrow = TRUE)
p13_vertices <- matrix(c(-1488705, 1369618,  
                         -1488999, 1369618,  
                         -1489008, 1369848, 
                         -1488705, 1369952,
                         -1488705, 1369618),
                       ncol = 2, byrow = TRUE) 
p14_vertices <- matrix(c(-1488804, 1369607,  
                         -1488999, 1369607,  
                         -1489008, 1369848, 
                         -1488804, 1369918,
                         -1488804, 1369607),
                       ncol = 2, byrow = TRUE)
p15_vertices <- matrix(c(-1488921, 1369704,  
                         -1489002, 1369704,  
                         -1489008, 1369848, 
                         -1488921, 1369878,
                         -1488921, 1369704),
                       ncol = 2, byrow = TRUE) 
p16_vertices <- matrix(c(-1488299, 1369901,  
                         -1488674, 1369901,  
                         -1488674, 1369962, 
                         -1488299, 1370091,
                         -1488299, 1369901),
                       ncol = 2, byrow = TRUE)  
p17_vertices <- matrix(c(-1488904, 1369779,  
                         -1489005, 1369779,  
                         -1489008, 1369848, 
                         -1488904, 1369884,
                         -1488904, 1369779),
                       ncol = 2, byrow = TRUE) 
p18_vertices <- matrix(c(-1488525, 1369926,  
                         -1488780, 1369926,  
                         -1488525, 1370013, 
                         -1488525, 1369926),
                       ncol = 2, byrow = TRUE) 

# 组合为sf数据框
sf_boxes <- st_sf(pix_value = c(1:18),
  geometry = st_sfc(st_polygon(list(p1_vertices)),
                    st_polygon(list(p2_vertices)),
                    st_polygon(list(p3_vertices)),
                    st_polygon(list(p4_vertices)),
                    st_polygon(list(p5_vertices)),
                    st_polygon(list(p6_vertices)),
                    st_polygon(list(p7_vertices)),
                    st_polygon(list(p8_vertices)),
                    st_polygon(list(p9_vertices)),
                    st_polygon(list(p10_vertices)),
                    st_polygon(list(p11_vertices)),
                    st_polygon(list(p12_vertices)),
                    st_polygon(list(p13_vertices)),
                    st_polygon(list(p14_vertices)),
                    st_polygon(list(p15_vertices)),
                    st_polygon(list(p16_vertices)),
                    st_polygon(list(p17_vertices)),
                    st_polygon(list(p18_vertices))),
  crs = 6350)

c) 数据可视化代码 #

ggplot() +
  geom_sf(data = polygon_extent, color = "gray") +
  geom_sf(data = sf_boxes, aes(color = as.character(pix_value),
                                  fill = as.character(pix_value)),
          alpha = 0.2) +
  NULL

d) 当前栅格化代码（未满足需求） #

dummy_raster <- st_rasterize(sf_boxes,
                          st_as_stars(st_bbox(polygon_extent),
                                      nx = 399,
                                      ny = 460))
plot(dummy_raster)

解决方案代码 #

步骤1：创建目标栅格模板 #

先定义与目标范围、分辨率一致的空栅格，确保栅格化的输出符合要求：

target_raster <- st_as_stars(st_bbox(polygon_extent), nx = 399, ny = 460)

步骤2：栅格化并处理重叠区域取最小值 #

使用st_rasterize的aggregation参数指定min函数，让重叠像素取最小的pix_value：

# 仅传入数值列进行栅格化，提升效率
rasterized_min <- st_rasterize(sf_boxes["pix_value"], target_raster, aggregation = min)

步骤3：最近邻插值填充空隙 #

用stars::st_interpolate_na的最近邻方法填充所有NA值，确保目标范围内无空隙：

final_raster <- st_interpolate_na(rasterized_min, method = "nearest")

可选：裁剪到目标多边形范围 #

如果需要严格限制在polygon_extent内（避免bbox边缘的插值超出目标多边形），可以添加裁剪步骤：

final_raster_cropped <- st_crop(final_raster, polygon_extent)

查看最终结果 #

plot(final_raster_cropped)

内容的提问来源于stack exchange，提问作者kseagull