#' @title Compute maps of spectral indices
#' @description Create maps of a set of spectral indices. Since
#' `gdal_calc.py` is used to perform computations, output files
#' are physical rasters (no output VRT is allowed).
#' @param infiles A vector of input filenames. Input files are paths
#' of BOA (or TOA) products already converted from SAFE format to a
#' format managed by GDAL (use [s2_translate] to do it);
#' their names must be in the sen2r naming convention
#' ([safe_shortname]).
#' @param indices Character vector with the names of the required
#' indices. Values should be included in names corresponding to the
#' Abbreviations of the following indices:
#' [IDB](https://www.indexdatabase.de/db/is.php?sensor_id=96).
#' @param outdir (optional) Full name of the output directory where
#' the files should be created (default: current directory).
#' `outdir` can bot be an existing or non-existing directory (in the
#' second case, its parent directory must exists).
#' If it is a relative path, it is expanded from the common parent
#' directory of `infiles`.
#' @param parameters (optional) Values of index parameters. This variable
#' must be a named list, in which each element is a list of parameters,
#' i.e.:
#' `parameters = list("SAVI" = list("a" = 0.5))`
#' Values can be both numeric values or band names (e.g. "band_1").
#' If not specified, parameters are set to default values.
#' @param source (optional) Vector with the products from which computing
#' the indices. It can be `"BOA"`, `"TOA"` or both (default). If both values
#' are provided, indices are computed from the available products (`"TOA"`
#' if TOA is available, BOA if BOA is available); in the case both are
#' available, two files are produced (they can be distinguished from the
#' level component - S2x1C or S2x2A - in the filename).
#' @param format (optional) Format of the output file (in a
#' format recognised by GDAL). Default is the same format of input images
#' (or "GTiff" in case of VRT input images).
#' @param subdirs (optional) Logical: if TRUE, different indices are
#' placed in separated `outfile` subdirectories; if FALSE, they are placed in
#' `outfile` directory; if NA (default), subdirectories are created only if
#' more than a single spectral index is required.
#' @param tmpdir (optional) Path where intermediate files (GTiff) will be
#' created in case `format` is `"VRT"`.
#' @param compress (optional) In the case a GTiff format is
#' present, the compression indicated with this parameter is used.
#' @param bigtiff (optional) Logical: if TRUE, the creation of a BigTIFF is
#' forced (default is FALSE).
#' This option is used only in the case a GTiff format was chosen.
#' @param dataType (optional) Numeric datatype of the output rasters.
#' if "Float32" or "Float64" is chosen, numeric values are not rescaled;
#' if `"Int16"` (default) or `"UInt16"`, values are multiplicated by `scaleFactor` argument;
#' if `"Byte"`, values are shifted by 100, multiplicated by 100 and truncated
#' at 200 (so that range -1 to 1 is coerced to 0-200), and nodata value
#' is assigned to 255.
#' @param scaleFactor (optional) Scale factor for output values when an integer
#' datatype is chosen (default values are 10000 for `"Int16"` and `"UInt16"`,
#' 1E9 for `"Int32"` and `"UInt32"`). Notice that, using `"UInt16"` and
#' `"UInt32"` types, negative values will be truncated to 0.
#' @param proc_mode (optional) Character: if `"gdal_calc"`,
#' `gdal_calc` routines are used to compute indices;
#' if `"raster"` (default) or `"stars"`, R functions are instead used
#' (using respectively `raster` or `stars` routines).
#' **Notes**:
#' 1. default value (`"raster"`) is the only fully supported mode.
#' `"gdal_calc"` can be used only if a runtime GDAL environment can be properly
#' configured (no assistance is provided in case of GDAL-related problems).
#' `"raster"` mode is experimental.
#' 2. There is a difference in which the two modes manage values out
#' of ranges (e.g. -32768 to 32767 in Int16 and 0 to 255 in Byte):
#' `"raster"` and `"stars"` modes set these values to NA,
#' `"gdal_calc"` mode clip them to the minimum/maximum values;
#' @param parallel (optional) Logical: if TRUE, the function is run using parallel
#' processing, to speed-up the computation for large rasters.
#' The number of cores is automatically determined; specifying it is also
#' possible (e.g. `parallel = 4`).
#' If FALSE (default), single core processing is used.
#' Multiprocess masking computation is always performed in singlecore mode
#' @param overwrite Logical value: should existing output files be
#' overwritten? (default: FALSE)
#' @param .log_message (optional) Internal parameter
#' (it is used when the function is called by `sen2r()`).
#' @param .log_output (optional) Internal parameter
#' (it is used when the function is called by `sen2r()`).
#' @return A vector with the names of the created products.
#' @export
#' @importFrom foreach foreach "%do%" "%dopar%"
#' @importFrom doParallel registerDoParallel
#' @importFrom parallel makeCluster stopCluster detectCores
#' @importFrom jsonlite fromJSON
#' @import data.table
#' @importFrom raster brick
#' @author Luigi Ranghetti, phD (2020)
#' @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
#' @examples
#' # Define file names
#' ex_in <- system.file(
#' "extdata/out/S2A2A_20190723_022_Barbellino_BOA_10.tif",
#' package = "sen2r"
#' )
#'
#' # Run function
#' ex_out <- s2_calcindices(
#' infiles = ex_in,
#' indices = "EVI",
#' outdir = tempdir(),
#' dataType = "Float32"
#' )
#' ex_out
#'
#' # Show output
#' oldpar <- par(mfrow = c(1,2), mar = rep(0,4))
#' image(stars::read_stars(ex_in), rgb = 4:2, maxColorValue = 3500, useRaster = TRUE)
#' par(mar = rep(2/3,4))
#' image(stars::read_stars(ex_out), useRaster = TRUE)
#' par(oldpar)
s2_calcindices <- function(
infiles,
indices,
outdir = ".",
parameters = NULL,
source = c("TOA","BOA"),
format = NA,
subdirs = NA,
tmpdir = NA,
compress = "DEFLATE",
bigtiff=FALSE,
dataType = "Int16",
scaleFactor = NA,
proc_mode = "raster",
parallel = FALSE,
overwrite = FALSE,
.log_message = NA,
.log_output = NA
) {
# to avoid NOTE on check
prod_type <- . <- i <- NULL
# Check proc_mode and GDAL external dependency
if (!proc_mode %in% c("gdal_calc", "raster", "stars")) {
print_message(
type = "warning",
"proc_mode = \"",proc_mode,"\" is not recognised; ",
"switching to \"raster\"."
)
proc_mode <- "raster"
}
if (proc_mode == "gdal_calc" && is.null(load_binpaths()$gdal_calc)) {
tryCatch(
check_gdal(abort = TRUE),
error = function(e) {
print_message(
type = "warning",
"External GDAL binaries are required with 'proc_mode = \"gdal_calc\"'; ",
"please configure them using function check_gdal() ",
"or through a GUI with check_sen2r_deps(). ",
"Now switching to proc_mode = \"raster\"."
)
proc_mode <- "raster"
}
)
}
# create outdir if not existing (and dirname(outdir) exists)
suppressWarnings(outdir <- expand_path(outdir, parent=comsub(infiles,"/"), silent=TRUE))
if (!dir.exists(dirname(outdir))) {
print_message(
type = "error",
"The parent folder of 'outdir' (",outdir,") does not exist; ",
"please create it."
)
}
dir.create(outdir, recursive=FALSE, showWarnings=FALSE)
# Compute n_cores
n_cores <- if (is.numeric(parallel)) {
min(as.integer(parallel), length(infiles))
} else if (parallel==FALSE) {
1
} else {
min(parallel::detectCores()-1, length(infiles), 11) # use at most 11 cores
}
if (n_cores<=1) {
`%DO%` <- `%do%`
parallel <- FALSE
n_cores <- 1
} else {
`%DO%` <- `%dopar%`
}
# generate indices.json if missing and read it
create_indices_db()
indices_db <- list_indices(
c("n_index","name","longname","s2_formula","a","b","c","d"),
all = TRUE
)
# check that the required indices exists
if (!all(indices %in% indices_db$name)) {
print_message(
type="error",
if (!any(indices %in% indices_db$name)) {"The "} else {"Some of the "},
"requested index names (\"",
paste(indices[!indices %in% indices_db$name], collapse="\", \""),
"\") are not recognisable; please use accepted ",
"values. To list accepted index names, type ",
"'sort(list_indices(\"name\", all=TRUE))'.")
}
if (!all(indices %in% indices_db$name)) {
print_message(
type="warning",
"Some of the specified index names (",
paste(indices[!indices %in% indices_db$name],collapse=", "),
") are not recognisable and will be skipped.")
indices <- indices[indices %in% indices_db$name]
}
# extract needed indices_db
indices_info <- indices_db[match(indices,indices_db$name),]
# check output format
gdal_formats <- fromJSON(system.file("extdata/settings/gdal_formats.json",package="sen2r"))$drivers
if (!is.na(format)) {
sel_driver <- gdal_formats[gdal_formats$name == format,]
if (nrow(sel_driver)==0) {
print_message(
type="error",
"Format \"",format,"\" is not recognised; ",
"please use one of the formats supported by your GDAL installation."#\n\n",
# "To list them, use the following command:\n",
# "gdalUtils::gdalinfo(formats=TRUE)\n\n",
# "To search for a specific format, use:\n",
# "gdalinfo(formats=TRUE)[grep(\"yourformat\", gdalinfo(formats=TRUE))]"
)
}
}
# assign scaleFactor and nodata values
if (grepl("Int",dataType) & is.na(scaleFactor)) {
scaleFactor <- ifelse(grepl("Int32",dataType),1E9,1E4)
}
sel_nodata <- switch(
dataType,
Int16=-2^15, UInt16=2^16-1, Int32=-2^31, UInt32=2^32-1,
Float32=-9999, Float64=-9999, Byte=255
)
# Get files metadata
infiles_meta <- sen2r_getElements(infiles)
infiles <- infiles[infiles_meta$prod_type %in% source]
infiles_meta <- infiles_meta[prod_type %in% source,]
# create subdirs (if requested)
if (is.na(subdirs)) {
subdirs <- ifelse(length(indices)>1, TRUE, FALSE)
}
if (subdirs) {
sapply(file.path(outdir,indices), dir.create, showWarnings=FALSE)
}
# read TOA/BOA image
if (n_cores > 1) { # nocov start
cl <- makeCluster(
n_cores,
type = if (Sys.info()["sysname"] == "Windows") {"PSOCK"} else {"FORK"}
)
registerDoParallel(cl)
print_message(
type = "message",
date = TRUE,
"Starting parallel computation of indices..."
)
} # nocov end
outfiles <- foreach(
i = seq_along(infiles),
.packages = c("raster","stars","sen2r"),
.combine=c,
.errorhandling="remove"
) %DO% {
# redirect to log files
if (n_cores > 1) { # nocov start
if (!is.na(.log_output)) {
sink(.log_output, split = TRUE, type = "output", append = TRUE)
}
if (!is.na(.log_message)) {
logfile_message = file(.log_message, open = "a")
sink(logfile_message, type="message")
}
} # nocov end
sel_infile <- infiles[i]
sel_infile_meta <- c(infiles_meta[i,])
sel_format <- suppressWarnings(
ifelse(!is.na(format), format, raster_metadata(sel_infile, "outformat", format = "list")[[1]]$outformat)
)
sel_out_ext <- gdal_formats[gdal_formats$name==sel_format,"ext"][1]
# check bands to use
if (sel_infile_meta$prod_type=="TOA") {
gdal_bands <- data.frame("letter"=LETTERS[1:12],"number"=1:12,"band"=paste0("band_",1:12), stringsAsFactors = FALSE)
} else if (sel_infile_meta$prod_type=="BOA") {
gdal_bands <- data.frame("letter"=LETTERS[1:11],"number"=1:11,"band"=paste0("band_",c(1:9,11:12)), stringsAsFactors = FALSE)
} else {
print_message(type="error", "Internal error (this should not happen).")
}
# compute single indices
# (this cycle is not parallelised)
sel_outfiles <- character(0)
for (j in seq_along(indices)) {
print_message(
type = "message",
date = TRUE,
paste0(
"Computing index ", indices[j],
" on date ",sel_infile_meta$sensing_date,"..."
)
)
# extract parameters
sel_parameters <- parameters[[indices[j]]]
# define output filename
sel_outfile <- paste0(
"S2",sel_infile_meta$mission,sel_infile_meta$level,"_",
strftime(sel_infile_meta$sensing_date,"%Y%m%d"),"_",
sel_infile_meta$id_orbit,"_",
sel_infile_meta$extent_name,"_",
indices_info[j,"name"],"_",
gsub("m$","",sel_infile_meta$res),".",
sel_out_ext)
# define subdir
out_subdir <- ifelse(subdirs, file.path(outdir,indices[j]), outdir)
# if output already exists and overwrite==FALSE, do not proceed
if (!file.exists(file.path(out_subdir,sel_outfile)) | overwrite==TRUE) {
# if indices[j] requires different parameters basing on sensor,
# retrieve the correct information (replace all except index name)
# (for now, this is done "by hand" since only CR indices requires it)
if (indices[j] %in% c("CRred","BDred","CRred2")) {
sel_index_info <- indices_db[indices_db$name==paste0(indices[j],"-2",sel_infile_meta$mission),]
sel_index_info$name <- indices[j]
} else {
sel_index_info <- indices_info[j,]
}
# change index formula to be used with bands
for (sel_par in c("a","b","c","d")) {
assign(sel_par, if (is.null(sel_parameters[[sel_par]])) {sel_index_info[,sel_par]} else {sel_parameters[[sel_par]]})
}
sel_formula <- sel_index_info[,"s2_formula"]
for (sel_par in c("a","b","c","d")) {
sel_formula <- gsub(paste0("([^0-9a-zA-Z])par\\_",sel_par,"([^0-9a-zA-Z])"),
paste0("\\1",get(sel_par),"\\2"),
sel_formula)
}
# Edit formula depending on proc_mode
if (proc_mode == "gdal_calc") {
for (sel_band in rev(seq_len(nrow(gdal_bands)))) {
sel_formula <- gsub(
gdal_bands[sel_band,"band"],
paste0("(",gdal_bands[sel_band,"letter"],".astype(float)/10000)"),
sel_formula
)
}
} else if (proc_mode == "raster") {
for (sel_band in rev(seq_len(nrow(gdal_bands)))) {
sel_formula <- gsub(
gdal_bands[sel_band,"band"],
paste0("(v[,",gdal_bands[sel_band,"number"],"]/10000)"),
sel_formula
)
}
} else if (proc_mode == "stars") {
for (sel_band in rev(seq_len(nrow(gdal_bands)))) {
sel_formula <- gsub(
gdal_bands[sel_band,"band"],
paste0("(v[",gdal_bands[sel_band,"number"],"]/10000)"),
sel_formula
)
}
}
# Clip within the datatype range
if (grepl("Int", dataType)) {
sel_formula <- paste0(
"clip(",
scaleFactor,"*(",sel_formula,"),",
switch(dataType, Int16=-2^15+2, UInt16=0, Int32=-2^31+4, UInt32=0),",",
switch(dataType, Int16=2^15-1, UInt16=2^16-2, Int32=2^31-3, UInt32=2^32-4),")"
)
}
if (dataType == "Byte") {
sel_formula <- paste0("clip(100+100*(",sel_formula,"),0,200)")
}
# if sel_format is VRT, create GTiff as intermediate source files
# (cannot create directly .tif files without breaking _req / _exi names)
if (sel_format == "VRT") {
# define and create tmpdir
if (is.na(tmpdir)) {
tmpdir <- file.path(out_subdir, ".vrt")
}
dir.create(tmpdir, recursive=FALSE, showWarnings=FALSE)
sel_format0 <- "GTiff"
sel_out_ext0 <- gdal_formats[gdal_formats$name==sel_format0,"ext"][1]
out_subdir0 <- tmpdir
sel_outfile0 <- gsub(paste0(sel_out_ext,"$"), sel_out_ext0, sel_outfile)
} else {
sel_format0 <- sel_format
out_subdir0 <- out_subdir
sel_outfile0 <- sel_outfile
}
# Launch the processing
if (proc_mode == "gdal_calc") {
gdalUtil(
"calc",
source = rep(sel_infile, nrow(gdal_bands)),
destination = file.path(out_subdir0,sel_outfile0),
formula = sel_formula,
options = c(
c(apply(gdal_bands,1,function(l){
c(paste0("--",l["letter"],"_band"), which(gdal_bands$letter==l["letter"]))
})),
"--type", dataType,
"--NoDataValue", sel_nodata,
"--format", sel_format0,
if (overwrite==TRUE) {"--overwrite"},
if (sel_format0 == "GTiff") {c(
"--co", paste0("COMPRESS=",toupper(compress)),
"--co", "TILED=YES"
)},
if (sel_format0=="GTiff" & bigtiff==TRUE) {c("--co", "BIGTIFF=YES")}
),
quiet = TRUE
)
} else if (proc_mode == "raster") {
calcindex_raster(
raster::brick(sel_infile),
sel_formula,
out_file = file.path(out_subdir0,sel_outfile0),
NAflag = sel_nodata,
sel_format = sel_format0,
compress = compress,
datatype = dataType,
overwrite = overwrite
)
} else if (proc_mode == "stars") {
calcindex_stars(
sel_infile,
sel_formula,
out_file = file.path(out_subdir0,sel_outfile0),
NAflag = sel_nodata,
sel_format = sel_format0,
compress = compress,
datatype = dataType,
overwrite = overwrite
)
}
# fix for envi extension (writeRaster use .envi)
if (sel_format0=="ENVI") {fix_envi_format(file.path(out_subdir0,sel_outfile0))}
if (sel_format == "VRT") {
gdalUtil(
"buildvrt",
source = file.path(out_subdir0,sel_outfile0),
destination = file.path(out_subdir,sel_outfile),
quiet = TRUE
)
}
} # end of overwrite IF cycle
sel_outfiles <- c(sel_outfiles, file.path(out_subdir,sel_outfile))
} # end of indices FOR cycle
# stop sinking
if (n_cores > 1) { # nocov start
if (!is.na(.log_output)) {
sink(type = "output")
}
if (!is.na(.log_message)) {
sink(type = "message"); close(logfile_message)
}
} # nocov end
sel_outfiles
} # end cycle on infiles
if (n_cores > 1) { # nocov start
stopCluster(cl)
print_message(
type = "message",
date = TRUE,
"Parallel computation of indices done."
)
} # nocov end
return(outfiles)
}
# Accessory functions to interpret NumPy functions power() and clip()
power <- function(x,y) {x^y}
clip <- function(x,min,max) {(x+min+2*max+abs(x-min)-abs(x+min-2*max+abs(x-min)))/4}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.