在 R 中光栅化重叠的 sf 对象并指定用于像素的值

问题描述 投票：0回答：1

我有一个 R 中包含多个多边形的 sf 数据框（每行都是一个正方形，位置略有不同，具有与其关联的不同值）。多边形之间也存在一些间隙。我想使用 star::st_rasterize 栅格化 sf 数据帧（使用第二个 sf 数据帧作为新栅格的范围），但遇到了 2 个问题。

  1. 当栅格像素中有多个重叠的 sf 对象时，我想指定我希望栅格像素从 sf 对象中获取尽可能低的值。

  2. 在光栅化过程中，我想根据最近的 sf 对象的值使用插值来填充 sf 对象之间的间隙。换句话说，我不希望栅格在第二个 sf df 范围内具有 NA 值。

这里有一些虚拟的 sf 数据框：

a）将成为栅格范围的多边形

### load libraries
library(tidyverse)
library(sf)
library(stars)

### set vertices
polygon_vertices <- matrix(c(-1489008, 1369848,
                             -1488191, 1370128,
                             -1488132, 1369347,
                             -1488985, 1369243,
                             -1489008, 1369848),
                           ncol = 2,
                           byrow = TRUE) 

### make it a polygon
polygon_extent <- st_polygon(list(polygon_vertices))

### make it an sf object and set the crs
polygon_extent <- polygon_extent %>%
  st_sfc(crs = 6350)

b) 具有重叠 sf 对象的 sf 数据框

### set the vertices for the rows of the sf df
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) 

### combine all of them into an sf df
sf_boxes <- st_sf(pix_value = c(1:18),  # Unique identifier for each polygon
  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)

所以看起来像这样：

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

然后我可以对其进行栅格化，但我无法弄清楚如何设置代码，以便 1）具有来自 sf_boxes 的重叠 sf 对象的每个栅格像素获得分配给它的最低值，2）范围内的缺失值的多边形范围已填充。

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

有任何想法吗？谢谢！

r gis spatial-interpolation rasterize
1个回答
0
投票

您可以使用 

terra::rasterize

来做到这一点

你的多边形

wkt <- c('POLYGON ((-1488355 1369737, -1488355 1370072, -1488237 1370112, -1488190 1370112, -1488162 1369737, -1488355 1369737))', 'POLYGON ((-1488355 1369737, -1488355 1370072, -1488237 1370112, -1488190 1370112, -1488162 1369737, -1488355 1369737))', 'POLYGON ((-1488253 1369484, -1488253 1369859, -1488171 1369859, -1488143 1369484, -1488253 1369484))', 'POLYGON ((-1488292 1369795, -1488292 1370093, -1488191 1370128, -1488166 1369795, -1488292 1369795))', 'POLYGON ((-1488260 1369634, -1488154 1369634, -1488132 1369347, -1488260 1369332, -1488260 1369634))', 'POLYGON ((-1488255 1369734, -1488630 1369734, -1488630 1369977, -1488255 1370106, -1488255 1369734))', 'POLYGON ((-1488240 1369419, -1488615 1369419, -1488615 1369794, -1488240 1369794, -1488240 1369419))', 'POLYGON ((-1488212 1370120, -1488191 1370128, -1488190 1370120, -1488212 1370120))', 'POLYGON ((-1488970 1369519, -1488995 1369519, -1489008 1369848, -1488970 1369861, -1488970 1369519))', 'POLYGON ((-1488970 1369519, -1488995 1369519, -1489008 1369848, -1488970 1369861, -1488970 1369519))', 'POLYGON ((-1488518 1369823, -1488518 1370016, -1488191 1370128, -1488168 1369823, -1488518 1369823))', 'POLYGON ((-1488452 1369986, -1488452 1370039, -1488191 1370128, -1488180 1369986, -1488452 1369986))', 'POLYGON ((-1488705 1369618, -1488999 1369618, -1489008 1369848, -1488705 1369952, -1488705 1369618))', 'POLYGON ((-1488804 1369607, -1488999 1369607, -1489008 1369848, -1488804 1369918, -1488804 1369607))', 'POLYGON ((-1488921 1369704, -1489002 1369704, -1489008 1369848, -1488921 1369878, -1488921 1369704))', 'POLYGON ((-1488299 1369901, -1488674 1369901, -1488674 1369962, -1488299 1370091, -1488299 1369901))', 'POLYGON ((-1488904 1369779, -1489005 1369779, -1489008 1369848, -1488904 1369884, -1488904 1369779))', 'POLYGON ((-1488525 1369926, -1488780 1369926, -1488525 1370013, -1488525 1369926))')

library(terra)
v <- vect(wkt)
values(v) <- data.frame(ID=1:nrow(v))

解决方案，使用 

terra::rasterize
和参数 
fun
和 
background

r <- rast(v, res=2)
x <- rasterize(v, r, "ID", fun="min", background=0)
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.