Close polygons using linestrings and selection of polygons created by area

user_2977164 0
Viewed 132

I'd like to close all the polygons using linestrings and the selection of polygons by area. In my example:

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

Now, I'd like to create new contours just only if the polygons are greater than 900 m^2, but in the first step I try to make the conversion on "LINESTRING" to "MULTIPOLYGON" (close polygons) using directly st_cast() and doesn't work. My final goal is to close the polygons, make the selections of polygons if it area is > 900m^2? Please, any help with it? Thanks in advance

sf r
2 Answers
# 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")
#

plot1

user_29771640

Assuming your example demonstrates the workflow you have, you could also directly polygonize the stars raster:

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

Created on 2022-01-13 by the reprex package (v2.0.1)

user_10820040
Related