#' @title Buffer cloud masks
#' @description Internal function (used by [s2_mask]) which smooths
#' and buffers a 0-1 mask image in order to reduce the roughness of the mask
#' obtained from SCL classification (which is done pixel by pixel).
#' See details.
#' @param inmask The path of the input 0-1 mask (where 0 represents the area
#' to be masked, 1 the clean surface).
#' @param tmpdir (optional) Path where intermediate files (VRT) will be created.
#' Default is a temporary directory.
#' @param radius (optional) Numerical (positive): the size (in the unit of
#' `inmask`, typically metres) to be used as radius for the smoothing
#' (the higher it is, the more smooth the output mask will result).
#' @param buffer (optional) Numerical (positive or negative): the size of the
#' buffer (in the unit of `inmask`, typically metres) to be applied to the
#' masked area after smoothing it (positive to enlarge, negative to reduce).
#' @param namask (optional) The path of an input 0-1 mask where 0 represents
#' the area of the original file with NA values (which should not be
#' smoothed / buffered).
#' Default (NULL) means that no NA values are present.
#' @param bigtiff (optional) Logical: if TRUE, the creation of a BigTIFF is
#' forced (default is FALSE).
#' @return The path of the smoothed mask.
#' @importFrom methods is
#' @author Luigi Ranghetti, phD (2019)
#' @references L. Ranghetti, M. Boschetti, F. Nutini, L. Busetto (2020).
#' "sen2r": An R toolbox for automatically downloading and preprocessing
#' Sentinel-2 satellite data. _Computers & Geosciences_, 139, 104473.
#' \doi{10.1016/j.cageo.2020.104473}, URL: \url{https://sen2r.ranghetti.info/}.
#' @note License: GPL 3.0
#' @keywords internal
smooth_mask <- function(
inmask,
tmpdir = tempdir(),
radius = 250,
buffer = 250,
namask = NULL,
bigtiff = FALSE
) {
# if inmask is a raster use the path (it should not happen)
inmask_path0 <- if (is(inmask, "character")) {
inmask
} else {
inmask@file@name
}
# if inmask is a raster use the path (it should not happen)
namask_path <- if (is.null(namask) | is(namask, "character")) {
namask
} else {
namask@file@name
}
# convert radius and buffer from metres to number of pixels
# (as required by gdal_fillnodata)
inmask_res <- mean(raster_metadata(inmask_path0, "res", format = "list")[[1]]$res)
radius_npx <- abs(radius / inmask_res)
buffer_npx <- buffer / inmask_res
# check GDAL external dependency
if (is.null(load_binpaths()$gdal_calc)) {
tryCatch(
check_gdal(abort = TRUE),
error = function(e) {
print_message(
type = "error",
"External GDAL binaries are required to smooth masks; ",
"please configure them using function check_gdal() ",
"or through a GUI with check_sen2r_deps()."
)
}
)
}
# 1. set inmask=1 (clear sky) and namask=0 (nodata) to NA
inmask_path1 <- file.path(tmpdir,basename(tempfile(pattern = "mask_", fileext = ".tif")))
if (!is.null(namask)) {
gdalUtil(
"calc",
source = c(inmask_path0, namask_path),
destination = inmask_path1,
formula = "A+1-B",
options = c(
"--type", "Byte", "--NoDataValue", "1", "--format", "GTiff",
"--co", "COMPRESS=LZW", if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
} else {
gdalUtil(
"translate",
source = inmask_path0,
destination = inmask_path1,
options = c(
"-of", "GTiff", "-co", "COMPRESS=LZW", "-a_nodata", "1",
if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
}
if (radius_npx != 0) {
# 2. first positive buffer (1/2 radius)
inmask_path2 <- gsub("\\.tif$","_2.tif",inmask_path1)
gdalUtil(
"fillnodata",
source = inmask_path1,
destination = inmask_path2,
options = c(
"-md", radius_npx*3/4,
"-co", "COMPRESS=LZW", "-si", "0", "-of", "GTiff"
),
quiet = TRUE
)
# 3. invert the image: set inmask=0 (clouds) and namask=0 (nodata) to NA
inmask_path3 <- gsub("\\.tif$","_3.tif",inmask_path1)
if (!is.null(namask)) {
inmask_path2b <- gsub("\\.tif$","_2b.tif",inmask_path1)
gdalUtil(
"translate",
source = inmask_path2,
destination = inmask_path2b,
options = c(
"-of", "GTiff", "-co", "COMPRESS=LZW", "-a_nodata", "none",
if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
gdalUtil(
"calc",
source = c(inmask_path2b, namask_path),
destination = inmask_path3,
formula = "A*B",
options = c(
"--type", "Byte", "--NoDataValue", "1", "--format", "GTiff",
"--co", "COMPRESS=LZW", if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
} else {
gdalUtil(
"translate",
source = inmask_path2,
destination = inmask_path3,
options = c(
"-of", "GTiff", "-co", "COMPRESS=LZW", "-a_nodata", "0",
if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
}
# 2. second negative buffer (3/4 radius + 5/4 radius)
inmask_path4 <- gsub("\\.tif$","_4.tif",inmask_path1)
gdalUtil(
"fillnodata",
source = inmask_path3,
destination = inmask_path4,
options = c(
"-md", radius_npx*2,
"-co", "COMPRESS=LZW", "-si", "0", "-of", "GTiff"
),
quiet = TRUE
)
# 5. invert the image: set inmask=1 (clear sky) and namask=0 (nodata) to NA
inmask_path5 <- gsub("\\.tif$","_5.tif",inmask_path1)
if (!is.null(namask)) {
inmask_path4b <- gsub("\\.tif$","_4b.tif",inmask_path1)
gdalUtil(
"translate",
source = inmask_path4,
destination = inmask_path4b,
options = c(
"-of", "GTiff", "-co", "COMPRESS=LZW", "-a_nodata", "none",
if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
gdalUtil(
"calc",
source = c(inmask_path4b, namask_path),
destination = inmask_path5,
formula = "A*B+1-B",
options = c(
"--type", "Byte", "--NoDataValue", "1", "--format", "GTiff",
"--co", "COMPRESS=LZW", if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
} else {
gdalUtil(
"translate",
source = inmask_path4,
destination = inmask_path5,
options = c(
"-of", "GTiff", "-co", "COMPRESS=LZW", "-a_nodata", "1",
if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
}
} else {
inmask_path5 <- inmask_path1
} # end of IF radius > 0 cycle
# 6. third positive buffer (5/4 radius to complete smooth, buffer to buffer if < 0, to 3/2 buffer if >0)
inmask_path6 <- gsub("\\.tif$","_6.tif",inmask_path1)
gdalUtil(
"fillnodata",
source = inmask_path5,
destination = inmask_path6,
options = c(
"-md", radius_npx*5/4+ifelse(buffer_npx>0,buffer_npx*3/2,buffer_npx),
"-co", "COMPRESS=LZW", "-si", "0", "-of", "GTiff"
),
quiet = TRUE
)
# 7-8 if buffer_npx > 0
if (buffer_npx > 0) {
# 7. invert the image: set inmask=0 (clouds) and namask=0 (nodata) to NA
inmask_path7 <- gsub("\\.tif$","_7.tif",inmask_path1)
if (!is.null(namask)) {
inmask_path6b <- gsub("\\.tif$","_6b.tif",inmask_path1)
gdalUtil(
"translate",
source = inmask_path6,
destination = inmask_path6b,
options = c(
"-of", "GTiff", "-co", "COMPRESS=LZW", "-a_nodata", "none",
if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
gdalUtil(
"calc",
source = c(inmask_path6b, namask_path),
destination = inmask_path7,
formula = "A*B",
options = c(
"--type", "Byte", "--NoDataValue", "1", "--format", "GTiff",
"--co", "COMPRESS=LZW", if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
} else {
gdalUtil(
"translate",
source = inmask_path6,
destination = inmask_path7,
options = c(
"-of", "GTiff", "-co", "COMPRESS=LZW", "-a_nodata", "0",
if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
}
# 8. fourth negative buffer (1/2)
inmask_path8 <- gsub("\\.tif$","_8.tif",inmask_path1)
gdalUtil(
"fillnodata",
source = inmask_path7,
destination = inmask_path8,
options = c(
"-md", buffer_npx/2,
"-co", "COMPRESS=LZW", "-si", "0", "-of", "GTiff"
),
quiet = TRUE
)
}
# 9. remove nodata labels
inmask_path9 <- gsub("\\.tif$","_9.tif",inmask_path1)
gdalUtil(
"translate",
source = ifelse(buffer_npx>0, inmask_path8, inmask_path6),
destination = inmask_path9,
options = c(
"-of", "GTiff", "-co", "COMPRESS=LZW", "-a_nodata", "none",
if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
if (!is.null(namask)) {
inmask_path9b <- gsub("\\.tif$","_9b.tif",inmask_path1)
gdalUtil(
"calc",
source = c(inmask_path9, namask_path),
destination = inmask_path9b,
formula = "A*B",
options = c(
"--type", "Byte", "--NoDataValue", "1", "--format", "GTiff",
"--co", "COMPRESS=LZW", if (bigtiff==TRUE) {c("-co", "BIGTIFF=YES")}
),
quiet = TRUE
)
inmask_path9b
} else {
inmask_path9
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.