Close polygons using linestrings and selection of polygons created by area
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
Solution 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")
#
Solution 2:
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)