使用线串关闭多边形并选择按区域创建的多边形答案

【问题标题】：Close polygons using linestrings and selection of polygons created by area使用线串关闭多边形并选择按区域创建的多边形
【发布时间】：2021-08-18 01:19:50
【问题描述】：

我想使用线串和按区域选择多边形来关闭所有多边形。在我的例子中：

# Packages
library(stars)
library(sf)
library(ggplot2)

#Vectorizing a raster object to an sf object
tif = system.file("tif/L7_ETMs.tif", package = "stars")
# Let's stay with only the first band, indicated in the final dimension
x = read_stars(tif)[, 1:50, 1:50, 1]
x = round(x/5)
# Calculate the min and max of the raster values
purrr::map(x, min)
# 10
purrr::map(x, max)
# 28
# Change values lower than 15 for 1 or 2
x[[1]] = ifelse(x[[1]]<15,1,2)

# Take a look at the image
plot(x)

# Obtain the contours
l =  st_contour(x, 
                # Remove NA
                na.rm = T,
                # Obtain contour lines instead of polygons
                contour_lines = TRUE,
                # raster values at which to draw contour levels
                breaks = 1.5)

# Plot the contours
raster_line <- st_as_sf(l, coords = c("x", "y"), crs = 32725) %>%
  st_cast("LINESTRING")
ggplot(raster_line) +
  geom_sf(aes(color =  L7_ETMs.tif), show.legend = "line")  +
  theme_bw() +
  theme(
   panel.grid.major = element_blank(),
   panel.grid.minor = element_blank(),
   legend.position = "none")
#

现在，我只想在多边形大于 900 m^2 的情况下创建新轮廓，但在第一步中，我尝试使用“LINESTRING”转换为“MULTIPOLYGON”（闭合多边形）直接st_cast() 并且不起作用。我的最终目标是关闭多边形，如果面积 > 900m^2，则选择多边形？请问，有什么帮助吗？提前致谢

【问题讨论】：

标签： r sf

【解决方案1】：

# Convert to polygon and multipolygons
p <- st_cast(raster_line,"POLYGON")
mult.p <- st_cast(p,"MULTIPOLYGON")

# Calculate polygons area
mult.p$area <- st_area(mult.p)

# Selection of polygons on interest
mult.p.target <-  mult.p %>% 
        dplyr::filter(as.numeric(mult.p$area) > 900) 
#

# Plot the new contours
raster_poly <- st_as_sf(mult.p.target, coords = c("x", "y"), crs = 32725) %>%
  st_cast("MULTIPOLYGON")
ggplot(raster_poly) +
  geom_sf(aes(color =  L7_ETMs.tif))  +
  theme_bw() +
  theme(
   panel.grid.major = element_blank(),
   panel.grid.minor = element_blank(),
   legend.position = "none")
#

【讨论】：

【解决方案2】：

假设您的示例演示了您所拥有的工作流程，您还可以直接对星星栅格进行多边形化：

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(sf)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

#Vectorizing a raster object to an sf object
tif = system.file("tif/L7_ETMs.tif", package = "stars")
# Let's stay with only the first band, indicated in the final dimension
x = read_stars(tif)[, 1:50, 1:50, 1]
x = round(x/5)

x %>% 
  cut(., breaks = c(14,Inf)) %>% 
  st_as_sf(., as_points = FALSE, merge = TRUE) %>% 
  mutate(area = st_area(.)) %>% 
  filter(area > units::set_units(900, "m^2", mode = "standard")) %>% 
  dplyr::select(-area) %>% 
  plot()

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

【讨论】：