raster_stack_to_unique_polygons.R

# method: raster stack to uniquely attributed polygons
# using GRASS 7.4 to avoid rasterToPolygons
library(sp)
library(sf)
library(raster)
library(rgrass7)
library(dsmartr) # just for the demo dataset 
# https://github.com/obrl-soil/dsmartr/blob/master/data/heronvale_covariates.rda
data("heronvale_covariates")

# the categorical rasters in heronvale_covariates are layers 3, 5, and 8 
# (geomorphon most-common-forms, MrVBF, QLD regional ecosystems)

hv_cat <- heronvale_covariates[[c(3,5,8)]]

# cast data to matrix, then get unique rows
# these are the unique combinations of attributes in use
ucombos <- unique(hv_cat) # MARGIN defaults to cell vector, convenient

# 338 unique combinations of these three rasters' values (from ~20K cells)
# cast to list
ucombos <- split(ucombos, row(ucombos))

# make a raster with values that match the index of ucombos (1-338) 
# NB this includes areas with any or all NA values
# even with clusterR this will take a while
beginCluster(n =  max(1, (parallel::detectCores() - 1)))
assign('ucombos', ucombos, envir = parent.frame())
cat_r <- clusterR(x = hv_cat,
                  fun = calc,
                  args = list(fun = function(cell) {
                    Position(function(x) identical(as.vector(cell), x), ucombos, nomatch = NA)
                  }),
                  export = 'ucombos')
endCluster()

### Polygonise 
# set up a temp location
# https://obrl-soil.github.io/r-osgeo4w-windows/ for setup stuff
initGRASS(gisBase = Sys.getenv('GISBASE'), tempdir(), override = T, mapset = 'PERMANENT')

# set working area and projection 
# (point at one of the on-disk sources for your stack)
execGRASS("g.proj", flags = c('t', 'c'), 
          parameters = list(georef = 'C:/data/HVC_ALL.tif'))

# read raster over to GRASS
writeRAST(as(cat_r, 'SpatialGridDataFrame'), 'cat_r')

# polygonise it (experiment with flag 's'!)
# https://grass.osgeo.org/grass74/manuals/r.to.vect.html
execGRASS('r.to.vect', input = 'cat_r', output = 'cat_v', type = 'area', flags = c('v', 'overwrite'))

# return vector layer to GRASS
cat_v <- readVECT('cat_v')

# cast to sf object
cat_sf <- st_as_sf(cat_v)

# 'cat' column can be matched to ucombos list pos, but I'm cheating a bit below b/c lazy
cat_sf$GM_MCF_pattern <- sapply(ucombos, function(x) unlist(x)[1], USE.NAMES = FALSE)
cat_sf$NETD_MRVBF <- sapply(ucombos, function(x) unlist(x)[2], USE.NAMES = FALSE)
cat_sf$PCRE_DBVG1M <- sapply(ucombos, function(x) unlist(x)[3], USE.NAMES = FALSE)

#  tidy up
cat_sf <- dplyr::select(cat_sf, -cat, -label) 
# ditch all-NA polygon(s) (ocean in this case)
keep <- apply(st_set_geometry(cat_sf, NULL), MARGIN = 1, FUN = function(x) all(is.na(x)))
cat_sf <- cat_sf[keep == FALSE, ]

plot(cat_sf['GM_MCF_pattern'], border = NA)
# compare with 
plot(hv_cat[[1]])

# can also cast to single polygons
cat_sfs <- st_cast(cat_sf, 'POLYGON')

#...?
# profit!