Makes thresholded sf multipolygons from user classifications for floating forests with a sample data set

threshold_classifications_sf.R

library(sf)
library(fasterize)
library(spex)
library(dplyr)
library(purrr)
library(ggplot2)
library(lwgeom)

test_set <- devtools::source_gist("c95701bd444cda3e342838fd9a090fb3",
                                  filename = "test_set.R")$value

#make into a raster at the 10m scale to that it doesn't kill spex
test_raster <- fasterize(test_set %>% mutate(value = 1), 
                         raster(test_set, res = 10), 
                         field = "value", fun = "sum")

#plot the raster
plot(test_raster)

#polygonize by threshold
make_threshold_raster <- function(arast, threshold){
  arast <-   arast >= threshold
  raster::values(arast)[!raster::values(arast)] <- NA
 # raster::values(arast)[raster::values(arast)] <- threshold
  arast
}

out_polys <- map(1:13, ~make_threshold_raster(test_raster, .x)) %>%
  map(polygonize) 

out_union <-  map(out_polys, st_union)

op_mp <- out_union %>%
  do.call(c, .) %>%
  st_sf(threshold = 1:13, geometry = .) %>%
  st_make_valid() 

#plot the result
ggplot(op_mp) +
  geom_sf() +
  facet_wrap(~threshold)