#' Apply minimum mapping unit to categorical raster
#'
#' @param x Raster object
#' @param mmu Minimum mapping unit as number of pixels
#' @param ncores Number of cores to use if parallel processing is wanted
#' @param filename Optional. filename for writing raster to disk.
#'
#' @return A raster layer with the designated mmu
#' @export
minimum_mapping_unit <- function(x, mmu, ncores, filename = NULL) {
layers <- raster::unstacK(raster::layerize(x, falseNA = TRUE))
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
patches <- foreach::foreach(i=layers, .packages=c("raster")) %dopar% {
patches <- raster::clump(i)
f <- as.data.frame(raster::freq(patches))
exclude_id <- f$value[which(f$count < mmu)]
rcl <- matrix(c(exclude_id, rep(NA, length(exclude_id))), ncol = 2)
reclassify(patches, rcl)
}
parallel::stopCluster(cl)
stacked <- do.call(raster::cover, patches)
masked <- raster::mask(x, stacked)
# apply majority filter
focal_fun <- function(v) {
mode_val <- as.numeric(names((which.max(table(v)))))
if (!purrr::is_empty(mode_val)) {
return(mode_val)
} else {
return(NA)
}
}
additional_nas <- raster::freq(masked, value = NA)- raster::freq(x, value = NA)
wsize <- 3
while(additional_nas > 0){
filled <- raster::focal(masked, matrix(1, nrow=wsize, ncol=wsize), fun=focal_fun, pad=TRUE, NAonly = TRUE)
filled_mask <- raster::mask(filled, x)
additional_nas <- raster::freq(filled_mask, value = NA) - raster::freq(x, value = NA)
wsize <- wsize + 2
}
if(!is.null(filename)) {
raster::writeRaster(filled_mask, filename, overwrite = TRUE)
}
return(filled_mask)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.