I am trying to figure out the best approach to convert an raster or stars object into its boundary returned as an sf object. I think I have a pretty decent handle on how to do it with raster (shown below):

#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
#> Loading required package: abind
#> Loading required package: sp

## path
test <- system.file("external/test.grd", package = "raster")

## Using raster
x <- raster(test)
r <- x > -Inf
pp <- rasterToPolygons(r, dissolve = TRUE)
#> Loading required namespace: rgeos
boundary_x <- st_as_sf(pp)

Here using stars i am forced to rely on st_union for that final aggregation step. For this example, it is pretty quick but for larger examples I first have to create many little polygons (using st_as_sf.stars then use st_union) which takes a long time. Is there a more efficient way of doing this?

## Using stars
y <- read_stars(test)
boundary_y <- st_as_sf(y, merge = TRUE) %>%

  • are you hoping to keep the exact boundary, with the ragged edge, or do you want a simplified geometry that contains everything (i.e.which might look "smoother" but doesn't match exactly the original grid, with a slightly larger area)? – stragu Dec 10 '20 at 7:43
  • @stragu simplified would totally be fine – boshek Dec 10 '20 at 19:04
  • interestingly, running a benchmark of your two alternatives with 100 iterations, the raster method is significantly slower than the stars method (median of 1232 ms vs 435 ms), so not sure where your issue comes from... I'm using sf 0.9-6, stars 0.4-3 and raster 3.4-5 (same geo libraries as you). – stragu Dec 11 '20 at 4:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.