R/s2_translate.R

In pobsteta/theia2r: An R toolbox to find, download and preprocess Sentinel-2 data

#' @title Convert from SAFE format
#' @description The function build a virtual raster from a Sentinel2 SAFE product,
#'  eventually translating it in another spatial format.
#'  For now, only L1C and L2a with long name (< 2016/12/06) are recognised.
#'  Output vrt is at 10m resolution.
#' @param infile Full path of the input SAFE folder (alternatively,
#'  full path of the xml file of the product with metadata).
#' @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 directory of `infile`.
#' @param subdirs (optional) Logical: if TRUE, differet output products are
#'  placed in separated `outdir` subdirectories; if FALSE, they are placed in
#'  `outdir` directory; if NA (default), subdirectories are created only if
#'  `prod_type` has length > 1.
#' @param tmpdir (optional) Path where intermediate files (VRT) will be created.
#'  Default is a temporary directory.
#' @param rmtmp (optional) Logical: should temporary files be removed?
#'  (Default: TRUE).
#'  This parameter takes effect only if the output files are not VRT
#'  (in this case temporary files cannot be deleted, because rasters of source
#'  bands are included within them).
#' @param prod_type (optional) Vector of types to be produced as outputs
#'  (see [safe_shortname] for the list of accepted values). Default is
#'  reflectance ('TOA' for level 1C, 'BOA' for level 2A).
#' @param tiles (optional) Character vector with the desired output tile IDs 
#'  (id specified IDs are not present in the input SAFE product, they are not
#'  produced). Default (NA) is to process all the found tiles.
#' @param res (optional) Spatial resolution (one between '10m', '20m' or '60m');
#'  default is '10m'. Notice that, choosing '10m' or '20m', bands with lower
#'  resolution will be rescaled to `res`. Band 08 is used with `res = '10m'`,
#'  band 08A with `res = '20m'` and `res = '60m'`.
#' @param format (optional) Format of the output file (in a
#'  format recognised by GDAL). Default value is 'VRT' (Virtual Raster).
#' @param compress (optional) In the case a GTiff format is
#'  chosen, the compression indicated with this parameter is used.
#' @param vrt_rel_paths (optional) Logical: if TRUE (default on Linux), 
#'  the paths present in the VRT output file are relative to the VRT position;
#'  if FALSE (default on Windows), they are absolute. 
#'  This takes effect only with `format = 'VRT'`.
#' @param utmzone (optional) UTM zone of output products (default:
#'  the first one retrieved from input granules). Note that this function
#'  does not perform reprojections: if no granules refer to the specified
#'  UTM zone, no output is created.
#' @param overwrite Logical value: should existing output files be
#'  overwritten? (default: FALSE)
#' @return A vector with the names of the created output files
#'   (just created or already existing).
#' @author Luigi Ranghetti, phD (2017) \email{ranghetti.l@@irea.cnr.it}
#' @note License: GPL 3.0
#' @importFrom jsonlite fromJSON
#' @export
#' @examples \dontrun{
#' s2_l1c_example <- file.path(
#'   '/existing/path',
#'   'S2A_MSIL1C_20170603T101031_N0205_R022_T32TQQ_20170603T101026.SAFE')
#' s2_l1c_example <- file.path(
#'   '/existing/path',
#'   'S2A_MSIL2A_20170603T101031_N0205_R022_T32TQQ_20170603T101026.SAFE')
#'
#' # Create a single TOA GeoTIFF in the same directory
#' s2_translate(s2_l1c_example, format='GTiff')
#'
#' # Create a single BOA VRT with a custom name
#' s2_translate(s2_l2a_example, '/new/path/example_sentinel2_sr.vrt',
#'   vrt_rel_paths=TRUE)
#'
#' # Create three products (ENVI) in the same directory at 60m resolution
#' s2_translate(s2_example, format='ENVI', prod_type=c('BOA','TCI','SCL'),
#'   res='60m', subdirs=TRUE)
#'}

s2_translate <- function(infile, outdir = ".", subdirs = NA, tmpdir = NA, rmtmp = TRUE, prod_type = NULL, tiles = NA, res = "10m", format = "VRT", compress = "DEFLATE", 
    vrt_rel_paths = NA, utmzone = "", overwrite = FALSE) {
    
    # Define vrt_rel_paths
    if (is.na(vrt_rel_paths)) {
        vrt_rel_paths <- Sys.info()["sysname"] != "Windows"
    }
    
    # Load GDAL paths
    binpaths <- load_binpaths("gdal")
    
    # check res (and use the resolutions >= specified res)
    if (!res %in% c("10m", "20m", "60m")) {
        print_message(type = "error", "\"res\" value is not recognised (accepted values are '10m', '20m' ", "and '60m').")
    }
    if (res == "10m") {
        res <- c("10m", "20m", "60m")
    } else if (res == "20m") {
        res <- c("20m", "60m")
    }
    
    # Check 'tiles' argument
    if (is.null(tiles)) {
        tiles <- NA
    }
    if (anyNA(tiles)) {
        tiles <- NA
    } else if (all(tiles == "")) {
        tiles <- NA
    }
    
    # check output format
    gdal_formats <- fromJSON(system.file("extdata", "gdal_formats.json", package = "theia2r"))
    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))]")
    }
    
    # Check GDAL installation
    check_gdal(abort = TRUE)
    # Retrieve xml required metadata
    infile_meta <- safe_getMetadata(infile, c("xml_main", "xml_granules", "utm", "level", "tiles", "jp2list"))
    infile_dir = dirname(infile_meta$xml_main)
    
    # define output directory
    suppressWarnings(outdir <- expand_path(outdir, parent = dirname(infile_dir), silent = TRUE))
    dir.create(outdir, recursive = FALSE, showWarnings = FALSE)
    # create subdirs
    if (is.na(subdirs)) {
        subdirs <- ifelse(length(prod_type) > 1, TRUE, FALSE)
    }
    if (subdirs) {
        sapply(file.path(outdir, prod_type), dir.create, showWarnings = FALSE)
    }
    
    # check compression value
    if (format == "GTiff") {
        if (!compress %in% c("JPEG", "LZW", "PACKBITS", "DEFLATE", "CCITTRLE", "CCITTFAX3", "CCITTFAX4", "LZMA", "NONE")) {
            print_message(type = "warning", "'", toupper(compress), "' is not a valid compression value; ", "the default 'DEFLATE' value will be used.")
            compress <- "DEFLATE"
        }
    }
    
    # retrieve UTM zone
    if (utmzone == "") {
        print_message(type = "message", "Using UTM zone ", sel_utmzone <- infile_meta$utm[1], ".")
    } else {
        sel_utmzone <- which(infile_meta$utm == as.integer(utmzone))
        if (length(sel_utmzone) == 0) {
            print_message(type = "warning", "Tiles with UTM zone ", utmzone, " are not present: zone ", sel_utmzone <- infile_meta$utm[1], " will be used.")
        }
    }
    
    # select default product type if missing
    if (is.null(prod_type)) {
        if (infile_meta$level == "1C") {
            prod_type <- "TOA"
        } else if (infile_meta$level == "2A") {
            prod_type <- "BOA"
        }
    }
    
    # define output extension
    out_ext <- sel_driver[1, "ext"]
    
    # create a file / set of files for each prod_type
    out_names <- character(0)  # names of created files
    for (sel_prod in prod_type) {
        try({
            
            if (sel_prod %in% c("BOA", "TOA")) {
                sel_type <- "MSI"
            } else {
                sel_type <- sel_prod
            }
            
            # define NA flag
            sel_na <- switch(sel_prod, BOA = "65535", TOA = "65535", SCL = "0", TCI = NA, NA)
            # define output subdir
            out_subdir <- ifelse(subdirs, file.path(outdir, sel_prod), outdir)
            
            # TODO check that required bands are present
            
            # define and create tmpdir
            if (is.na(tmpdir)) {
                tmpdir <- if (format == "VRT") {
                  rmtmp <- FALSE  # force not to remove intermediate files
                  if (!missing(outdir)) {
                    file.path(outdir, ".vrt")
                  } else {
                    tempfile(pattern = "s2translate_")
                  }
                } else {
                  tempfile(pattern = "s2translate_")
                }
            }
            dir.create(tmpdir, recursive = FALSE, showWarnings = FALSE)
            
            # cycle on granules (with compact names, this runs only once; with old name, one or more)
            for (sel_granule in infile_meta$xml_granules) {
                try({
                  
                  sel_tile <- safe_getMetadata(dirname(sel_granule), "nameinfo")$id_tile
                  
                  # continue only if sel_tile is within desired tiles
                  if (anyNA(tiles) | sel_tile %in% tiles) 
                    {
                      
                      # define output basename
                      out_prefix <- safe_shortname(sel_granule, prod_type = sel_prod, res = res[1], full.name = FALSE, abort = TRUE)
                      # complete output filename
                      out_name <- file.path(out_subdir, paste0(out_prefix, ".", out_ext))
                      
                      # if out_name already exists and overwrite==FALSE, do not proceed
                      if (!file.exists(out_name) | overwrite == TRUE) 
                        {
                          
                          # select required bands from the list and order them by resolution
                          jp2df_selbands <- infile_meta$jp2list[infile_meta$jp2list$type == sel_type & infile_meta$jp2list$tile == sel_tile, ]
                          jp2df_selbands <- jp2df_selbands[with(jp2df_selbands, order(band, res)), ]
                          # remove lower resolutions and keep only the best resolution for each band
                          if (!any(jp2df_selbands$res == "")) {
                            jp2df_selbands <- jp2df_selbands[as.integer(substr(jp2df_selbands$res, 1, 2)) >= as.integer(substr(res[1], 1, 2)), ]
                          } else {
                            # for oldname L1C (which do not have 'res' in granule name) remove B08 or B8A basing on the requested 'res'
                            if (as.integer(substr(res[1], 1, 2)) < 20) {
                              jp2df_selbands <- jp2df_selbands[-grep("B8A.jp2$", jp2df_selbands$layer), ]
                            } else {
                              jp2df_selbands <- jp2df_selbands[-grep("B08.jp2$", jp2df_selbands$layer), ]
                            }
                          }
                          jp2df_selbands <- jp2df_selbands[!duplicated(with(jp2df_selbands, paste(band, tile))), ]
                          # extract vector of paths
                          jp2_selbands <- file.path(infile_dir, jp2df_selbands[, "relpath"])
                          
                          # create final vrt with all the bands (of select final raster with a single band)
                          if (length(jp2_selbands) > 1) {
                            final_vrt_name <- ifelse(format == "VRT", out_name, paste0(tmpdir, "/", out_prefix, ".vrt"))
                            system(paste0(binpaths$gdalbuildvrt, " -separate ", "-resolution highest ", "\"", final_vrt_name, "\" ", paste(paste0("\"", jp2_selbands, 
                              "\""), collapse = " ")), intern = Sys.info()["sysname"] == "Windows")
                            if (vrt_rel_paths == TRUE) {
                              gdal_abs2rel(final_vrt_name)
                            }
                          } else {
                            final_vrt_name <- jp2_selbands
                          }
                          
                          # create output file (or copy vrt file)
                          if (format != "VRT" | length(jp2_selbands) == 1) {
                            system(paste0(binpaths$gdal_translate, " -of ", format, " ", if (format == "GTiff") {
                              paste0("-co COMPRESS=", toupper(compress), " ")
                            }, if (!is.na(sel_na)) {
                              paste0("-a_nodata ", sel_na, " ")
                            }, "\"", final_vrt_name, "\" ", "\"", out_name, "\""), intern = Sys.info()["sysname"] == "Windows")
                            if (format == "VRT" & vrt_rel_paths == TRUE) {
                              gdal_abs2rel(out_name)
                            }
                          }
                          
                        }  # end of 'overwite' IF cycle
                      
                      out_names <- c(out_names, out_name)
                      
                    }  # end of 'sel_tile %in% tiles' IF cycle
                  
                })
            }  # end of sel_granule cycle
            
        })
    }  # end of prod_type cycle
    
    # Remove temporary files
    if (rmtmp == TRUE) {
        unlink(tmpdir, recursive = TRUE)
    }
    
    print_message(type = "message", length(out_names), " output files were correctly created.")
    return(out_names)
    
}