How to preserve holes when combining contours into polygons via R?

Part of the workflow for my project requires creating contours from a raster. I am able to generate the contours using raster::rasterToContour and can convert to a polygon using the sf package in R, but am having some trouble getting the desired outcome.

Here is a scaled-down example of what I am currently doing:

library(raster)
library(sf)

raster <- raster(matrix(c(rep(0,5),
                          0,1,1,1,0,
                          0,1,0,1,0,
                          0,1,1,1,0,
                          rep(0,5)), nrow = 5, ncol = 5))
plot(raster)

contour <- rasterToContour(raster, level = 0.5)
plot(contour, add = TRUE)

sf <- st_as_sf(contour) %>% st_polygonize()
plot(sf, col = "gray50", border = "black")

The initial contour plot seems to work appropriately for what I need. However, I can't seem to figure out how to avoid the central region becoming filled, as shown here. I would like the central portion to become a hole in the larger polygon so that it looks like this. I'm sure I can achieve this by pulling the polygons into GIS software, but it would be helpful to do so within R.

Solution

You could try st_cast("POLYGON") instead of st_polygonize():

library(raster)
library(stars)
library(sf)

m <- matrix(c(rep(0,5),
              0,1,1,1,0,
              0,1,0,1,0,
              0,1,1,1,0,
              rep(0,5)), nrow = 5, ncol = 5)

contour_1 <- 
  raster(m) |>
  rasterToContour(level = 0.5) |>
  st_as_sf() |> 
  st_cast("POLYGON")
contour_1
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 0.2 ymin: 0.2 xmax: 0.8 ymax: 0.8
#> CRS:           NA
#>     level                       geometry
#> C_1   0.5 POLYGON ((0.3 0.2, 0.2 0.3,...
plot(contour_1)

Or perhaps stars and stars::st_contour():

contour_2 <- 
  st_as_stars(m) |> 
  st_contour()
contour_2
#> Simple feature collection with 2 features and 3 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 0 ymin: 0 xmax: 5 ymax: 5
#> CRS:           NA
#>           A1  Min Max                       geometry
#> 1 [-0.5,0.5) -0.5 0.5 MULTIPOLYGON (((4.5 5, 5 5,...
#> 2  [0.5,1.5)  0.5 1.5 MULTIPOLYGON (((4.000001 3....
plot(contour_2[,"A1"], key.pos = NULL)

^{Created on 2024-02-09 with reprex v2.1.0}