Mercurial > repos > ecology > srs_process_data
diff Lib_preprocess_S2.r @ 0:cf69ad260611 draft default tip
planemo upload for repository https://github.com/Marie59/Sentinel_2A/srs_tools commit b32737c1642aa02cc672534e42c5cb4abe0cd3e7
| author | ecology |
|---|---|
| date | Mon, 09 Jan 2023 13:36:02 +0000 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Lib_preprocess_S2.r Mon Jan 09 13:36:02 2023 +0000 @@ -0,0 +1,1335 @@ +# == == == == == == == == == == == == == == == == == == == == == == == == == == == +# preprocS2 +# Lib_preprocess_S2.R +# == == == == == == == == == == == == == == == == == == == == == == == == == == == +# PROGRAMMERS: +# Jean-Baptiste FERET <jb.feret@teledetection.fr> +# Copyright 2021/08 Jean-Baptiste FERET +# == == == == == == == == == == == == == == == == == == == == == == == == == == == +# This Library contains functions to preprocess Sentinel-2 images downloaded from +# different data hubs, such as THEIA, PEPS or SCIHUB +# == == == == == == == == == == == == == == == == == == == == == == == == == == == + +#" This function adjusts information from ENVI header +#" +#" @param dsn character. path where to store the stack +#" @param bands list. should include "bandname", and if possible "wavelength" +#" @param sensor character. Name of the sensor used to acquire the image +#" @param stretch boolean. Set TRUE to get 10% stretching at display for reflectance, mentioned in hdr only +#" +#" @return None +#" @importFrom utils read.table +#" @importFrom raster hdr raster +#" @export +adjust_envi_hdr <- function(dsn, bands, sensor = "Unknown", stretch = FALSE) { + + # Edit hdr file to add metadata + hdr <- read_envi_header(get_hdr_name(dsn)) + hdr$`band names` <- bands$bandname + if (length(bands$wavelength) == length(bands$bandname)) { + hdr$wavelength <- bands$wavelength + }else { + hdr$wavelength <- NULL + } + if (stretch == TRUE) { + hdr$`default stretch` <- "0.000000 1000.000000 linear" + } + hdr$`z plot range` <- NULL + hdr$`data ignore value` <- "-Inf" + hdr$`sensor type` <- sensor + write_envi_header(hdr = hdr, hdrpath = get_hdr_name(dsn)) + + # remove unnecessary files + file2remove <- paste(dsn, ".aux.xml", sep = "") + if (file.exists(file2remove)) file.remove(file2remove) + file2remove <- paste(dsn, ".prj", sep = "") + if (file.exists(file2remove)) file.remove(file2remove) + file2remove <- paste(dsn, ".stx", sep = "") + if (file.exists(file2remove)) file.remove(file2remove) + return(invisible()) +} + +#" This function saves reflectance files +#" +#" @param s2sat character. Sentinel-2 mission ("2A" or "2B") +#" @param tile_s2 character. S2 tile name (2 numbers + 3 letters) +#" @param dateacq_s2 double. date of acquisition +#" +#" @return s2mission character. name of the S2 mission (2A or 2B) +#" @importFrom sen2r safe_getMetadata check_scihub_connection s2_list +#" @export +check_s2mission <- function(s2sat, tile_s2, dateacq_s2) { + + # is mission already defined by user? + if (!is.null(s2sat)) { + if (s2sat == "2A") { + s2mission <- "2A" + }else if (s2sat == "2B") { + s2mission <- "2B" + }else { + message("Could not identify if image from Sentinel-2A or -2B") + message("Defining central wavelength of spectral bands based on S2A") + s2mission <- "2A" + } + }else { + message("Could not identify if image from Sentinel-2A or -2B") + message("Defining central wavelength of spectral bands based on S2A") + s2mission <- "2A" + } + return(s2mission) +} + +#" this function aims at computing directory size +#" @param path character. path for directory +#" @param recursive boolean . set T if recursive +#" +#" @return size_files numeric. size in bytes +#" - image stack +#" - path for individual band files corresponding to the stack +#" - path for vector (reprojected if needed) +#" +#" @importFrom raster raster +#" @importFrom tools file_path_sans_ext file_ext +#" @export +dir_size <- function(path, recursive = TRUE) { + stopifnot(is.character(path)) + files <- list.files(path, full.names = TRUE, recursive = recursive) + vect_size <- sapply(files, function(x) file.size(x)) + size_files <- sum(vect_size) + return(size_files) +} + +#" This function reads S2 data from L2A directories downloaded from +#" various data hubs including THEIA, PEPS & SCIHUB (SAFE format & LaSRC) +#" @param path_dir_s2 character. path for S2 directory +#" @param path_vector character. path for vector file +#" @param s2source character. type of directory format (depends on atmospheric correction: SAFE produced from Sen2Cor) +#" @param resolution numeric. buffer applied to vector file (in meters) +#" @param interpolation character. method for resampling. default = "bilinear" +#" @param fre_sre character. SRE or FRE products from THEIA +#" +#" @return listout list. +#" - image stack +#" - path for individual band files corresponding to the stack +#" - path for vector (reprojected if needed) +#" +#" @importFrom raster raster +#" @importFrom tools file_path_sans_ext file_ext +#" @export +extract_from_s2_l2a <- function(path_dir_s2, path_vector = NULL, s2source = "SAFE", + resolution = 10, interpolation = "bilinear", fre_sre = "FRE") { + # Get list of paths corresponding to S2 bands and depending on S2 directory + s2_bands <- get_s2_bands(path_dir_s2 = path_dir_s2, + s2source = s2source, + resolution = resolution, + fre_sre = fre_sre) + + if (length(s2_bands$s2bands_10m) > 0) { + rastmp <- raster::raster(s2_bands$s2bands_10m[[1]]) + } else if (length(s2_bands$s2bands_20m) > 0) { + rastmp <- raster::raster(s2_bands$s2bands_20m[[1]]) + } + # check if vector and raster share the same projection. if not, re-project vector + if (!is.null(path_vector)) { + raster_proj <- raster::projection(rastmp) + path_vector_reproj <- paste(tools::file_path_sans_ext(path_vector), "_reprojected.shp", sep = "") + path_vector <- reproject_shp(path_vector_init = path_vector, + newprojection = raster_proj, + path_vector_reproj = path_vector_reproj) + } + # Extract data corresponding to the vector footprint (if provided) & resample data if needed + if (length(s2_bands$s2bands_10m) > 0) { + stack_10m <- read_s2bands(s2_bands = s2_bands$s2bands_10m, path_vector = path_vector, + resampling = 1, interpolation = interpolation) + } + if (length(s2_bands$s2bands_20m) > 0) { + if (resolution == 10 && s2source != "LaSRC") { + resampling <- 2 + }else { + resampling <- 1 + } + stack_20m <- read_s2bands(s2_bands = s2_bands$s2bands_20m, path_vector = path_vector, + resampling = resampling, interpolation = interpolation) + } + # get full stack including 10m and 20m spatial resolution + if (length(s2_bands$s2bands_10m) > 0 && length(s2_bands$s2bands_20m) > 0) { + diffxstart <- attributes(stack_10m)$dimensions[[1]]$from - attributes(stack_20m)$dimensions[[1]]$from + diffxstop <- attributes(stack_10m)$dimensions[[1]]$to - attributes(stack_20m)$dimensions[[1]]$to + diffystart <- attributes(stack_10m)$dimensions[[2]]$from - attributes(stack_20m)$dimensions[[2]]$from + diffystop <- attributes(stack_10m)$dimensions[[2]]$to - attributes(stack_20m)$dimensions[[2]]$to + if (!diffxstop == 0) { + # size of 20m > size of 10m --> reduce 20m + # size of 10m > size of 20m --> reduce 10m + if (diffxstop > 0) { + stack_10m <- stack_10m[, 1:(dim(stack_10m)[1] - diffxstop), , ] + }else if (diffxstop < 0) { + stack_20m <- stack_20m[, 1:(dim(stack_20m)[1] + diffxstop), , ] + } + } + if (!diffystop == 0) { + if (diffystop > 0) { + stack_10m <- stack_10m[, , 1:(dim(stack_10m)[2] - diffystop), ] + }else if (diffystop < 0) { + stack_20m <- stack_20m[, , 1:(dim(stack_20m)[2] + diffystop), ] + } + } + if (!diffxstart == 0) { + if (diffxstart > 0) { + stack_20m <- stack_20m[, (1 + diffxstart):dim(stack_20m)[1], , ] + }else if (diffxstart < 0) { + stack_10m <- stack_10m[, (1 - diffxstart):dim(stack_10m)[1], , ] + } + } + if (!diffystart == 0) { + if (diffystart > 0) { + stack_20m <- stack_20m[, , (1 + diffystart):dim(stack_20m)[2], ] + }else if (diffystart < 0) { + stack_10m <- stack_10m[, , (1 - diffystart):dim(stack_10m)[2], ] + } + } + # reorder bands with increasing wavelength + s2bands <- c("B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12", "Cloud") + namebands <- c(names(s2_bands$s2bands_10m), names(s2_bands$s2bands_20m)) + reorder_bands <- match(s2bands, namebands) + namebands <- namebands[reorder_bands] + listfiles <- c(stack_10m$attr, stack_20m$attr)[reorder_bands] + + # adjust size to initial vector footprint without buffer + # --> buffer is needed in order to ensure that extraction following + # footprint of vector matches for images of different spatial resolution + # get bounding box corresponding to footprint of image or image subset + bb_xycoords <- get_bb(path_raster = listfiles[1], + path_vector = path_vector, buffer = 0) + + # prepare reading data for extent defined by bounding box + nxoff <- bb_xycoords$UL$col + nyoff <- bb_xycoords$UL$row + nxsize <- bb_xycoords$UR$col - bb_xycoords$UL$col + 1 + nysize <- bb_xycoords$LR$row - bb_xycoords$UR$row + 1 + nbufxsize <- nxsize + nbufysize <- nysize + s2_stack <- stars::read_stars(listfiles, along = "band", + RasterIO = list(nXOff = nxoff, nYOff = nyoff, + nXSize = nxsize, nYSize = nysize, + nBufXSize = nbufxsize, nBufYSize = nbufysize, + resample = "nearest_neighbour"), proxy = TRUE) + + + names(s2_stack$attr) <- namebands + }else if (length(s2_bands$s2bands_10m) > 0) { + s2_stack <- stack_10m + namebands <- names(s2_bands$s2bands_10m) + names(s2_stack$attr) <- namebands + }else if (length(s2_bands$s2bands_20m) > 0) { + s2_stack <- stack_20m + namebands <- names(s2_bands$s2bands_20m) + names(s2_stack$attr) <- namebands + } + + listout <- list("s2_stack" = s2_stack, "s2_bands" = s2_bands, "path_vector" = path_vector, + "namebands" = namebands) + return(listout) +} + +#" This function gets coordinates of a bounding box defined by a vector (optional) and a raster +#" +#" @param path_raster character. path for raster file +#" @param path_vector character. path for vector file +#" @param buffer numeric. buffer applied to vector file (in meters) +#" +#" @return bb_xycoords list. Coordinates (in pixels) of the upper/lower right/left corners of bounding box +#" @export +get_bb <- function(path_raster, path_vector = NULL, buffer = 0) { + + if (!is.null(path_vector)) { + # get bounding box with a 50m buffer in order to allow for interpolation + bb_xycoords <- get_bb_from_vector(path_raster = path_raster, + path_vector = path_vector, + buffer = buffer) + }else if (is.null(path_vector)) { + bb_xycoords <- get_bb_from_fullimage(path_raster) + } + return(bb_xycoords) +} + +#" This function gets extreme coordinates of a bounding box corresponding to a full image +#" +#" @param path_raster character. path for raster file +#" +#" @return bb_xycoords list. Coordinates (in pixels) of the upper/lower right/left corners of bounding box +#" @importFrom raster raster +#" @export +get_bb_from_fullimage <- function(path_raster) { + # get raster coordinates corresponding to Full image + rasterobj <- raster::raster(path_raster) + bb_xycoords <- list() + bb_xycoords[["UL"]] <- data.frame("row" = 1, "col" = 1) + bb_xycoords[["UR"]] <- data.frame("row" = 1, "col" = dim(rasterobj)[2]) + bb_xycoords[["LL"]] <- data.frame("row" = dim(rasterobj)[1], "col" = 1) + bb_xycoords[["LR"]] <- data.frame("row" = dim(rasterobj)[1], "col" = dim(rasterobj)[2]) + return(bb_xycoords) +} + +#" This gets bounding box corresponding to a vector from a raster (UL, UR, LL, LR corners) +#" +#" @param path_raster character. path for raster file +#" @param path_vector character. path for vector file +#" @param buffer numeric. buffer applied to vector file (in meters) +#" +#" @return bb_xycoords list. Coordinates (in pixels) of the upper/lower right/left corners of bounding box +#" @importFrom sf st_read st_bbox st_crop +#" @importFrom rgeos gbuffer bbox2SP +#" @importFrom sp SpatialPoints bbox +#" @importFrom raster projection extract extent raster +#" @importFrom methods as +#" @export +get_bb_from_vector <- function(path_raster, path_vector, buffer = 0) { + + data_raster <- raster::raster(path_raster) + # extract BB coordinates from vector + bb_vector <- rgeos::gbuffer(spgeom = as(sf::st_read(dsn = path_vector, quiet = TRUE), "Spatial"), + width = buffer, byid = TRUE) + # extract BB coordinates from raster + bb_raster <- rgeos::bbox2SP(bbox = bbox(data_raster)) + # compute intersection + intersect <- rgeos::gIntersection(bb_vector, bb_raster) + bbext <- raster::extent(intersect) + xmin <- bbext[1] + xmax <- bbext[2] + ymin <- bbext[3] + ymax <- bbext[4] + # get coordinates of bounding box corresponding to vector + corners <- list() + corners[["UR"]] <- sp::SpatialPoints(coords = cbind(xmax, ymax)) + corners[["LR"]] <- sp::SpatialPoints(coords = cbind(xmax, ymin)) + corners[["UL"]] <- sp::SpatialPoints(coords = cbind(xmin, ymax)) + corners[["LL"]] <- sp::SpatialPoints(coords = cbind(xmin, ymin)) + raster::projection(corners[["UL"]]) <- raster::projection(corners[["UR"]]) <- + raster::projection(corners[["LL"]]) <- raster::projection(corners[["LR"]]) <- + raster::projection(sf::st_read(dsn = path_vector, quiet = TRUE)) + # get coordinates for corners of bounding box + bb_xycoords <- list() + for (corner in names(corners)) { + ex_df <- as.data.frame(raster::extract(data_raster, corners[[corner]], cellnumbers = TRUE)) + colrow <- ind2sub(data_raster, ex_df$cell) + bb_xycoords[[corner]] <- data.frame("row" = colrow$row, "col" = colrow$col) + } + return(bb_xycoords) +} + +#" get hdr name from image file name, assuming it is BIL format +#" +#" @param impath path of the image +#" +#" @return corresponding hdr +#" @import tools +#" @export +get_hdr_name <- function(impath) { + if (tools::file_ext(impath) == "") { + impathhdr <- paste(impath, ".hdr", sep = "") + }else if (tools::file_ext(impath) == "bil") { + impathhdr <- gsub(".bil", ".hdr", impath) + }else if (tools::file_ext(impath) == "zip") { + impathhdr <- gsub(".zip", ".hdr", impath) + }else { + impathhdr <- paste(tools::file_path_sans_ext(impath), ".hdr", sep = "") + } + + if (!file.exists(impathhdr)) { + message("WARNING : COULD NOT FIND hdr FILE") + print(impathhdr) + message("Process may stop") + } + return(impathhdr) +} + +#" This function returns path for the spectral bands to be used +#" +#" @param path_dir_s2 character. Path for the directory containing S2 data. either L2A .SAFE S2 file or THEIA directory +#" @param s2source character. defines if data comes from SciHub as SAFE directory, from THEIA or from LaSRC +#" @param resolution numeric. spatial resolution of the final image: 10m or 20m +#" @param fre_sre character. SRE or FRE products from THEIA +#" +#" @return listbands list. contains path for spectral bands corresponding to 10m and 20m resolution +#" @export +get_s2_bands <- function(path_dir_s2, s2source = "SAFE", resolution = 10, fre_sre = "FRE") { + + if (s2source == "SAFE" || s2source == "Sen2Cor") { + listbands <- get_s2_bands_from_sen2cor(path_dir_s2 = path_dir_s2, resolution = resolution) + }else if (s2source == "THEIA") { + listbands <- get_s2_bands_from_theia(path_dir_s2 = path_dir_s2, resolution = resolution, + fre_sre = fre_sre) + }else if (s2source == "LaSRC") { + listbands <- get_s2_bands_from_lasrc(path_dir_s2 = path_dir_s2, resolution = resolution) + }else { + message("The data source (Atmospheric correction) for Sentinel-2 image is unknown") + message("Please provide S2 images from one of the following data sources:") + message("- LaSRC (atmospheric correction: LaSRC)") + message("- THEIA (atmospheric correction: MAJA)") + message("- SAFE (atmospheric correction: Sen2Cor)") + s2bands_10m <- s2bands_20m <- granule <- mtdfile <- metadata_msi <- metadata_lasrc <- NULL + listbands <- list("s2bands_10m" = s2bands_10m, "s2bands_20m" = s2bands_20m, "GRANULE" = granule, + "metadata" = mtdfile, "metadata_MSI" = metadata_msi, + "metadata_lasrc" = metadata_lasrc) + } + return(listbands) +} + +#" This function returns path for the spectral bands in SAFE / sen2Cor directory +#" +#" @param path_dir_s2 character. Path for the SAFE directory containing S2 data +#" @param resolution numeric. spatial resolution of the final image: 10m or 20m +#" +#" @return listbands list. contains path for spectral bands corresponding to 10m and 20m resolution, as well name of as granule +#" @export +get_s2_bands_from_sen2cor <- function(path_dir_s2, resolution = 10) { + # build path for all bands + if (resolution == 10) { + b10m <- c("B02", "B03", "B04", "B08") + b20m <- c("B05", "B06", "B07", "B8A", "B11", "B12") + }else { + b10m <- c() + b20m <- c("B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12") + } + # get granule directory & path for corresponding metadata XML file + granule <- list.dirs(list.dirs(path_dir_s2, recursive = FALSE)[grep(pattern = "GRANULE", + x = list.dirs(path_dir_s2, recursive = FALSE))], recursive = FALSE) + mtdfile <- file.path(granule, "MTD_TL.xml") + if (file.exists(file.path(path_dir_s2, "MTD_MSIL2A.xml"))) { + mtd_msi_file <- file.path(path_dir_s2, "MTD_MSIL2A.xml") + } else { + mtd_msi_file <- NULL + } + + # Define path for bands + s2bands_20m_dir <- file.path(granule, "IMG_DATA", "R20m") + s2bands_10m_dir <- file.path(granule, "IMG_DATA", "R10m") + s2bands_10m <- s2bands_20m <- list() + for (band in b20m) { + s2bands_20m[[band]] <- file.path(s2bands_20m_dir, list.files(s2bands_20m_dir, pattern = band)) + } + for (band in b10m) { + s2bands_10m[[band]] <- file.path(s2bands_10m_dir, list.files(s2bands_10m_dir, pattern = band)) + } + # get cloud mask + cloud <- "MSK_CLDPRB_20m" + cloud_20m_dir <- file.path(granule, "QI_DATA") + s2bands_20m[["Cloud"]] <- file.path(cloud_20m_dir, list.files(cloud_20m_dir, pattern = cloud)) + listbands <- list("s2bands_10m" = s2bands_10m, + "s2bands_20m" = s2bands_20m, + "GRANULE" = granule, + "metadata" = mtdfile, + "metadata_MSI" = mtd_msi_file, + "metadata_lasrc" = NULL) + return(listbands) +} + +#" This function returns path for the spectral bands in LaSRC directory +#" +#" @param path_dir_s2 character. Path for the SAFE directory containing S2 data +#" @param resolution numeric. spatial resolution of the final image: 10m or 20m +#" +#" @return listbands list. contains path for spectral bands corresponding to 10m and 20m resolution, as well name of as granule +#" @importFrom stringr str_subset +#" @export +get_s2_bands_from_lasrc <- function(path_dir_s2, resolution = 10) { + + # get granule directory & path for corresponding metadata XML file + granule <- path_dir_s2 + mtdfile <- file.path(granule, "MTD_TL.xml") + if (file.exists(file.path(path_dir_s2, "MTD_MSIL1C.xml"))) { + mtd_msi_file <- file.path(path_dir_s2, "MTD_MSIL1C.xml") + } else { + mtd_msi_file <- NULL + } + + # build path for all bands + b10m <- c("band2", "band3", "band4", "band5", "band6", "band7", "band8", "band8a", "band11", "band12") + b10m_standard <- c("B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12") + # Define path for bands + s2bands_10m <- s2bands_20m <- list() + for (i in 1:seq_along(b10m)) { + s2bands_10m[[b10m_standard[i]]] <- file.path(path_dir_s2, + list.files(path_dir_s2, + pattern = paste(b10m[i], ".tif", sep = ""))) + } + + # get metadata file containing offset + mtd_lasrc <- str_subset(list.files(path_dir_s2, pattern = "S2"), ".xml$") + if (file.exists(file.path(path_dir_s2, mtd_lasrc))) { + metadata_lasrc <- file.path(path_dir_s2, mtd_lasrc) + } else { + metadata_lasrc <- NULL + } + # get cloud mask + cloud <- "CLM" + s2bands_10m[["Cloud"]] <- file.path(path_dir_s2, list.files(path_dir_s2, pattern = cloud)) + listbands <- list("s2bands_10m" = s2bands_10m, + "s2bands_20m" = s2bands_20m, + "GRANULE" = granule, + "metadata" = mtdfile, + "metadata_MSI" = mtd_msi_file, + "metadata_lasrc" = metadata_lasrc) + return(listbands) +} + +#" This function returns path for the spectral bands in THEIA directory +#" +#" @param path_dir_s2 character. Path for the SAFE directory containing S2 data +#" @param resolution numeric. spatial resolution of the final image: 10m or 20m +#" @param fre_sre character. SRE or FRE products from THEIA +#" +#" @return listbands list. contains path for spectral bands corresponding to 10m and 20m resolution, as well name of as granule +#" @export +get_s2_bands_from_theia <- function(path_dir_s2, resolution = 10, fre_sre = "FRE") { + + # build path for all bands + if (resolution == 10) { + b10m <- c("B02", "B03", "B04", "B08") + b20m <- c("B05", "B06", "B07", "B8A", "B11", "B12") + } else { + b10m <- c() + b20m <- c("B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12") + } + + # get path_tile_s2 directory & path for corresponding metadata XML file + path_tile_s2 <- list.dirs(path_dir_s2, recursive = FALSE) + files_tile_s2 <- list.files(path_tile_s2, recursive = FALSE) + mtdfile <- file.path(path_tile_s2, files_tile_s2[grep(pattern = "MTD_ALL.xml", x = files_tile_s2)]) + + # Define path for bands + s2bands_10m_dir <- s2bands_20m_dir <- path_tile_s2 + s2bands_10m <- s2bands_20m <- list() + for (band in b20m) { + band_20m_pattern <- paste0(gsub("0", "", band), ".tif") # for THEAI band 2 is "B2" ("B02" for SAFE) + list_files_20m <- list.files(s2bands_20m_dir, pattern = band_20m_pattern) + s2bands_20m[[band]] <- file.path(s2bands_20m_dir, list_files_20m)[grep(pattern = fre_sre, + x = file.path(s2bands_20m_dir, list_files_20m))] + } + for (band in b10m) { + band_10m_pattern <- paste0(gsub("0", "", band), ".tif") # for THEAI band 2 is "B2" ("B02" for SAFE) + list_files_10m <- list.files(s2bands_10m_dir, pattern = band_10m_pattern) + s2bands_10m[[band]] <- file.path(s2bands_10m_dir, list_files_10m)[grep(pattern = fre_sre, + x = file.path(s2bands_10m_dir, list_files_10m))] + } + + # get cloud mask 10m + cloud_10m <- "CLM_R1" + cloud_10m_dir <- file.path(path_tile_s2, "MASKS") + s2bands_10m[["Cloud"]] <- file.path(cloud_10m_dir, list.files(cloud_10m_dir, pattern = cloud_10m)) + + # get cloud mask 20m + cloud_20m <- "CLM_R2" + cloud_20m_dir <- file.path(path_tile_s2, "MASKS") + s2bands_20m[["Cloud"]] <- file.path(cloud_20m_dir, list.files(cloud_20m_dir, pattern = cloud_20m)) + + # return list bands + listbands <- list("s2bands_10m" = s2bands_10m, + "s2bands_20m" = s2bands_20m, + "path_tile_s2" = path_tile_s2, + "metadata" = mtdfile) + return(listbands) +} + +#" This function check S2 data level: +#" - L2A: already atmospherically corrected +#" - L1C: requires atmospheric corrections with sen2cor +#" +#" @param prodname character. original name for the S2 image +#" +#" @return s2level character. S2 level: L1C or L2A +#" @export +get_s2_level <- function(prodname) { + prodname <- basename(prodname) + if (length(grep(pattern = "L1C_", x = prodname)) == 1) { + s2level <- "L1C" + } else if (length(grep(pattern = "L2A_", x = prodname)) == 1) { + s2level <- "L2A" + } + return(s2level) +} + +#" This function gets tile from S2 image +#" +#" @param prodname character. original name for the S2 image +#" +#" @return tilename character +#" @importFrom tools file_path_sans_ext +#" @export +get_tile <- function(prodname) { + prodname <- basename(prodname) + tilename <- tools::file_path_sans_ext(gsub("_.*", "", gsub(".*_T", "", prodname))) + return(tilename) +} + +#" This function gets acquisition date from S2 image +#" +#" @param prodname character. original name for the S2 image +#" +#" @return dateacq character +#" @export +get_date <- function(prodname) { + prodname <- basename(prodname) + dateacq <- as.Date(gsub("T.*", "", gsub(".*_20", "20", prodname)), format = "%Y%m%d") + return(dateacq) +} + +#" download S2 L1C data from Copernicus hub or Google cloud +#" +#" @param list_safe safe object. produced with sen2r::s2_list +#" @param l1c_path character. path for storage of L1C image +#" @param path_vector path for a vector file +#" @param time_interval dates. time interval for S2 query +#" @param googlecloud boolean. set to TRUE if google cloud SDK is installed and +#" @param forcegoogle boolean. set to TRUE if only google requested +#" sen2r configured as an alternative hub for S2 download +#" +#" @return prodname character. S2 Product name +#" @importFrom sen2r safe_is_online s2_list s2_download s2_order check_gcloud +#" @export +get_s2_l1c_image <- function(list_safe, l1c_path, path_vector, time_interval, + googlecloud = FALSE, forcegoogle = FALSE) { + # Check if available from Copernicus hub first + copernicus_avail <- sen2r::safe_is_online(list_safe) + # if available: download + prodname <- attr(list_safe, which = "name") + if (file.exists(file.path(l1c_path, prodname))) { + message("L1C file already downloaded") + message(file.path(l1c_path, prodname)) + } else { + if (copernicus_avail == TRUE && forcegoogle == FALSE) { + sen2r::s2_download(list_safe, outdir = l1c_path) + } else if (copernicus_avail == FALSE || forcegoogle == TRUE) { + # if not available and googlecloud==TRUE + if (googlecloud == TRUE) { + # check if google cloud SDK available from this computer + ggc <- sen2r::check_gcloud() + if (ggc == TRUE) { + message("downloading from Google cloud") + list_safe_ggc <- sen2r::s2_list(spatial_extent = sf::st_read(dsn = path_vector), + time_interval = time_interval, + server = "gcloud") + prodname <- attr(list_safe_ggc, which = "name") + if (file.exists(file.path(l1c_path, prodname))) { + message("L1C file already downloaded") + message(file.path(l1c_path, prodname)) + } else { + sen2r::s2_download(list_safe_ggc, outdir = l1c_path) + # check if QI_DATA exists in DATASTRIP, and create it if not the case + datastrip_path <- file.path(l1c_path, prodname, "DATASTRIP") + dsdir <- list.dirs(datastrip_path, recursive = FALSE) + if (length(match(list.dirs(dsdir, recursive = FALSE, full.names = FALSE), "QI_DATA")) == 0) { + dir.create(file.path(dsdir, "QI_DATA")) + } + } + } else if (ggc == FALSE) { + message("googlecloud set to TRUE but missing") + message("Please install Google cloud SDK") + message("https://cloud.google.com/sdk/docs/install") + message("and/or set configuration of sen2r following instructions") + message("https://www.r-bloggers.com/2021/06/downloading-sentinel-2-archives-from-google-cloud-with-sen2r/") + } + } + } + if (copernicus_avail == FALSE && googlecloud == FALSE) { + message("S2 image in Long Term Archive (LTA)") + message("Ordering image from LTA") + message("This may take 1 day, please run your script later") + orders2 <- sen2r::s2_order(list_safe) + message("An alternative is possible with Google cloud SDK") + message("https://cloud.google.com/sdk/docs/install") + message("and/or set configuration of sen2r following instructions") + message("https://www.r-bloggers.com/2021/06/downloading-sentinel-2-archives-from-google-cloud-with-sen2r/") + } + } + return(prodname) +} + +#" download S2 L2A data from Copernicus hub or convert L1C to L2A +#" +#" @param l2a_path character. path for storage of L2A image +#" @param spatial_extent path for a vector file +#" @param dateacq character. date of acquisition +#" @param deletel1c Boolean. set TRUE to delete L1C images +#" @param Sen2Cor Boolean. set TRUE to automatically perform atmospheric corrections using sen2Cor +#" @param googlecloud boolean. set to TRUE if google cloud SDK is installed and +#" sen2r configured as an alternative hub for S2 download +#" +#" @return pathl2a character. Path for L2A image +#" @importFrom sen2r s2_list s2_download +#" @importFrom R.utils getAbsolutePath + +#" @export +get_s2_l2a_image <- function(l2a_path, spatial_extent, dateacq, + deletel1c = FALSE, sen2cor = TRUE, + googlecloud = FALSE) { + + # Needs to be updated: define path for L1c data + l1c_path <- l2a_path + # define time interval + time_interval <- as.Date(c(dateacq, dateacq)) + # get list S2 products corresponding to study area and date of interest using sen2r package + if (googlecloud == TRUE) { + server <- c("scihub", "gcloud") + } else if (googlecloud == FALSE) { + server <- "scihub" + } + list_safe <- sen2r::s2_list(spatial_extent = sf::st_read(dsn = spatial_extent), + time_interval = time_interval, + server = server, availability = "check") + # download products + sen2r::s2_download(list_safe, outdir = l2a_path) + # name all products + prodname <- attr(list_safe, which = "name") + prodfullpath <- file.path(l2a_path, prodname) + if (sen2cor == TRUE) { + for (imgname in prodname) { + s2level <- get_s2_level(imgname) + if (s2level == "L1C") { + datepattern <- gsub(pattern = "-", replacement = "", x = dateacq) + pathl2a <- s2_from_l1c_to_l2a(prodname = imgname, l1c_path = l2a_path, l2a_path = l2a_path, + datepattern = datepattern, tmp_path = NULL) + if (deletel1c == TRUE) { + unlink(x = R.utils::getAbsolutePath(file.path(l1c_path, prodname)), + recursive = TRUE, force = TRUE) + # delete from full path and add atmospherically corrected + whichimg <- grep(x = prodfullpath, pattern = imgname) + dateacq <- get_date(imgname) + tilename <- get_tile(imgname) + pathl2a <- list.files(path = l2a_path, pattern = tilename, full.names = TRUE) + pathl2a <- pathl2a[grep(x = pathl2a, pattern = dateacq)] + pathl2a <- pathl2a[grep(x = basename(pathl2a), pattern = "L2A")] + prodfullpath[whichimg] <- pathl2a + } + } + } + } + + return(prodfullpath) +} + +#" convert image coordinates from index to X-Y +#" +#" @param Raster image raster object +#" @param image_index coordinates corresponding to the raster +ind2sub <- function(data_raster, image_index) { + c <- ((image_index - 1) %% data_raster@ncols) + 1 + r <- floor((image_index - 1) / data_raster@ncols) + 1 + my_list <- list("col" = c, "row" = r) + return(my_list) +} + +#" mosaicing a set of rasters +#" +#" @param list_rasters character. list of paths corresponding to rasters to mosaic +#" @param dst_mosaic character. path and name of mosaic produced +#" @param stretch boolean. Set TRUE to get 10% stretching at display for reflectance, mentioned in hdr only +#" +#" @return None +#" @importFrom gdalUtils mosaic_rasters +#" @importFrom raster hdr raster +#" @export +mosaic_rasters <- function(list_rasters, dst_mosaic, stretch = FALSE) { + + # produce mosaic + gdalUtils::mosaic_rasters(gdalfile = list_rasters, dst_dataset = dst_mosaic, + separate = FALSE, of = "Ehdr", verbose = TRUE) + + # convert hdr to ENVI format + raster::hdr(raster(dst_mosaic), format = "ENVI") + # add info to hdr based on initial rasters + hdr_init <- read_envi_header(get_hdr_name(list_rasters[1])) + hdr <- read_envi_header(get_hdr_name(dst_mosaic)) + hdr$`band names` <- hdr_init$`band names` + hdr$wavelength <- hdr_init$wavelength + if (stretch == TRUE) { + hdr$`default stretch` <- "0.000000 1000.000000 linear" + } + hdr$`z plot range` <- NULL + hdr$`data ignore value` <- "-Inf" + hdr$`sensor type` <- hdr_init$`sensor type` + hdr$`coordinate system string` <- read.table(paste(file_path_sans_ext(dst_mosaic), ".prj", sep = "")) + write_envi_header(hdr = hdr, hdrpath = get_hdr_name(dst_mosaic)) + return(invisible()) +} + +#" Reads ENVI hdr file +#" +#" @param hdrpath Path of the hdr file +#" +#" @return list of the content of the hdr file +#" @export +read_envi_header <- function(hdrpath) { + if (!grepl(".hdr$", hdrpath)) { + stop("File extension should be .hdr") + } + hdr <- readLines(hdrpath) + ## check ENVI at beginning of file + if (!grepl("ENVI", hdr[1])) { + stop("Not an ENVI header (ENVI keyword missing)") + } else { + hdr <- hdr [-1] + } + ## remove curly braces and put multi-line key-value-pairs into one line + hdr <- gsub("\\{([^}]*)\\}", "\\1", hdr) + l <- grep("\\{", hdr) + r <- grep("\\}", hdr) + + if (length(l) != length(r)) { + stop("Error matching curly braces in header (differing numbers).") + } + + if (any(r <= l)) { + stop("Mismatch of curly braces in header.") + } + + hdr[l] <- sub("\\{", "", hdr[l]) + hdr[r] <- sub("\\}", "", hdr[r]) + + for (i in rev(seq_along(l))) { + hdr <- c( + hdr [seq_len(l [i] - 1)], + paste(hdr [l [i]:r [i]], collapse = "\n"), + hdr [-seq_len(r [i])] + ) + } + + ## split key = value constructs into list with keys as names + hdr <- sapply(hdr, split_line, "=", USE.NAMES = FALSE) + names(hdr) <- tolower(names(hdr)) + + ## process numeric values + tmp <- names(hdr) %in% c( + "samples", "lines", "bands", "header offset", "data type", + "byte order", "default bands", "data ignore value", + "wavelength", "fwhm", "data gain values" + ) + hdr [tmp] <- lapply(hdr [tmp], function(x) { + as.numeric(unlist(strsplit(x, ", "))) + }) + + return(hdr) +} + +#" This function reads a list of files corresponding to S2 bands +#" S2 bands are expected to have uniform spatial resolution and footprint +#" @param s2_bands list. list of S2 bands obtained from get_s2_bands +#" @param path_vector path for a vector file +#" @param resampling numeric. resampling factor (default = 1, set to resampling = 2 to convert 20m into 10m resolution) +#" @param interpolation character. method for resampling. default = "bilinear" +#" +#" @return stack_s2 list. contains stack of S2 bands +#" +#" @importFrom stars read_stars +#" @importFrom sf st_bbox st_read st_crop +#" @export + +read_s2bands <- function(s2_bands, path_vector = NULL, + resampling = 1, interpolation = "bilinear") { + # get bounding box corresponding to footprint of image or image subset + bb_xycoords <- get_bb(path_raster = s2_bands[[1]], + path_vector = path_vector, buffer = 50) + + # prepare reading data for extent defined by bounding box + nxoff <- bb_xycoords$UL$col + nyoff <- bb_xycoords$UL$row + nxsize <- bb_xycoords$UR$col - bb_xycoords$UL$col + 1 + nysize <- bb_xycoords$LR$row - bb_xycoords$UR$row + 1 + nbufxsize <- resampling * nxsize + nbufysize <- resampling * nysize + if (resampling == 1) { + interpolation <- "nearest_neighbour" + } + # write interpolated individual bands in temp directory + tmpdir <- tempdir() + tmpfile <- list() + for (band in names(s2_bands)) { + stack_s2_tmp <- stars::read_stars(s2_bands[[band]], along = "band", + RasterIO = list(nXOff = nxoff, nYOff = nyoff, + nXSize = nxsize, nYSize = nysize, + nBufXSize = nbufxsize, nBufYSize = nbufysize, + resample = interpolation), proxy = FALSE) + if (!is.null(path_vector)) { + stack_s2_tmp <- sf::st_crop(x = stack_s2_tmp, y = st_bbox(st_read(dsn = path_vector, quiet = TRUE))) + } + tmpfile[[band]] <- file.path(tmpdir, tools::file_path_sans_ext(basename(s2_bands[[band]]))) + if (band == "Cloud") { + stars::write_stars(obj = stack_s2_tmp, dsn = tmpfile[[band]], + driver = "ENVI", type = "Byte", overwrite = TRUE) + } else { + stars::write_stars(obj = stack_s2_tmp, dsn = tmpfile[[band]], + driver = "ENVI", type = "Int16", overwrite = TRUE) + } + gc() + } + + stack_s2 <- stars::read_stars(tmpfile, along = "band", proxy = TRUE) + return(stack_s2) +} + +#" This function reads a raster stack, and gets footprint as pixel coordinates or vector file as input +#" @param path_raster character. path for raster file +#" @param path_vector character. path for vector file +#" @param bbpix list. coordinates of pixels corresponding to a bounding box +#" +#" @return starsobj stars object corresponding to raster or raster subset +#" +#" @importFrom stars read_stars +#" @importFrom sf st_bbox st_read st_crop +#" @export +read_raster <- function(path_raster, path_vector = NULL, bbpix = NULL) { + # get bounding box corresponding to footprint of image or image subset + if (is.null(bbpix)) { + bb_xycoords <- get_bb(path_raster = path_raster, + path_vector = path_vector, buffer = 0) + } else { + bb_xycoords <- bbpix + } + # prepare reading data for extent defined by bounding box + nxoff <- bb_xycoords$UL$col + nyoff <- bb_xycoords$UL$row + nxsize <- bb_xycoords$UR$col - bb_xycoords$UL$col + 1 + nysize <- bb_xycoords$LR$row - bb_xycoords$UR$row + 1 + nbufxsize <- nxsize + nbufysize <- nysize + starsobj <- stars::read_stars(path_raster, along = "band", + RasterIO = list(nXOff = nxoff, nYOff = nyoff, + nXSize = nxsize, nYSize = nysize, + nBufXSize = nbufxsize, nBufYSize = nbufysize), + proxy = FALSE) + return(starsobj) +} + +#" This function reprojects a shapefile and saves reprojected shapefile +#" +#" @param path_vector_init character. path for a shapefile to be reprojected +#" @param newprojection character. projection to be applied to path_vector_init +#" @param path_vector_reproj character. path for the reprojected shapefile +#" +#" @return path_vector character. path of the shapefile +#" - path_vector_init if the vector did not need reprojection +#" - path_vector_reproj if the vector needed reprojection +#" +#" @importFrom rgdal readOGR writeOGR +#" @importFrom sp spTransform +#" @importFrom raster projection +#" @export +reproject_shp <- function(path_vector_init, newprojection, path_vector_reproj) { + + dir_vector_init <- dirname(path_vector_init) + # shapefile extension + fileext <- file_ext(basename(path_vector_init)) + if (fileext == "shp") { + name_vector_init <- file_path_sans_ext(basename(path_vector_init)) + vector_init_ogr <- rgdal::readOGR(dir_vector_init, name_vector_init, verbose = FALSE) + } else if (fileext == "kml") { + vector_init_ogr <- rgdal::readOGR(path_vector_init, verbose = FALSE) + } + vector_init_proj <- raster::projection(vector_init_ogr) + + if (!vector_init_proj == newprojection) { + dir_vector_reproj <- dirname(path_vector_reproj) + name_vector_reproj <- file_path_sans_ext(basename(path_vector_reproj)) + vector_reproj <- sp::spTransform(vector_init_ogr, newprojection) + rgdal::writeOGR(obj = vector_reproj, dsn = dir_vector_reproj, layer = name_vector_reproj, + driver = "ESRI Shapefile", overwrite_layer = TRUE) + path_vector <- path_vector_reproj + } else { + path_vector <- path_vector_init + } + return(path_vector) +} + + +#" perform atmospheric corrections to convert L1C to L2A data with Sen2cor +#" +#" @param prodname character. produced with sen2r::s2_list +#" @param l1c_path character. path of directory where L1C image is stored +#" @param l2a_path character. path of directory where L2A image is stored +#" @param datepattern character. pattern corresponding to date of acquisition to identify L2A directory +#" @param tmp_path character. path of temporary directory where L2A image is stored +#" sen2r configured as an alternative hub for S2 download +#" +#" @return pathl2a character. S2 Product name +#" @importFrom sen2r safe_is_online s2_list s2_download s2_order +#" @importFrom R.utils getAbsolutePath +#" +#" @export +s2_from_l1c_to_l2a <- function(prodname, l1c_path, l2a_path, datepattern, tmp_path = NULL) { + + # define path for tmp directory + if (is.null(tmp_path)) { + tmp_path <- tempdir(check = TRUE) + } + tmp_prodlist <- prodname + # perform Sen2Cor atmospheric corrections + binpath <- sen2r::load_binpaths() + # 2- open a command prompt and directly run sen2cor with following command line + cmd <- paste(binpath$sen2cor, + "--output_dir", R.utils::getAbsolutePath(l2a_path), + R.utils::getAbsolutePath(file.path(l1c_path, prodname)), sep = " ") + system(cmd) + pathl2a <- list.files(path = l2a_path, pattern = datepattern, full.names = TRUE) + + return(pathl2a) +} + +#" This function saves cloud masks. +#" "cloudMask_Binary" is default binary mask with 0 where clouds are detected and 1 for clean pixels +#" "cloudMask_RAW" is the original cloud layer produced by atmospheric correction algorithm +#" --> may be useful to refine cloud mask +#" +#" @param s2_stars list. stars object containing raster data. Can be produced with function extract_from_s2_l2a +#" @param cloud_path character. +#" @param s2source character. +#" @param footprint character. path for vector file defining footprint of interest in the image +#" @param saveraw boolean. should the original cloud mask layer be saved? +#" @param maxchunk numeric. Size of individual chunks to be written (in Mb) +#" +#" @return list of cloudmasks (binary mask, and raw mask if required) +#" @importFrom sf st_read +#" @importFrom stars write_stars +#" @importFrom raster raster +#" @export +save_cloud_s2 <- function(s2_stars, cloud_path, s2source = "SAFE", + footprint = NULL, saveraw = FALSE, maxchunk = 256) { + + whichcloud <- which(names(s2_stars$attr) == "Cloud") + # Save cloud mask + if (saveraw == TRUE) { + cloudraw <- file.path(cloud_path, "CloudMask_RAW") + obj <- stars::read_stars(s2_stars$attr[whichcloud], proxy = TRUE) + sizeobj <- dim(obj)[1] * dim(obj)[2] / (1024**2) + nbchunks <- ceiling(sizeobj / maxchunk) + stars::write_stars(obj, + dsn = cloudraw, + driver = "ENVI", + type = "Byte", + chunk_size = c(dim(obj)[1], dim(obj)[2] / nbchunks), + progress = TRUE) + } else { + cloudraw <- NULL + } + # Save cloud mask as in biodivMapR (0 = clouds, 1 = pixel ok) + cloudmask <- stars::read_stars(s2_stars$attr[whichcloud], proxy = FALSE) + if (s2source == "SAFE" || s2source == "THEIA") { + cloudy <- which(cloudmask[[1]] > 0) + sunny <- which(cloudmask[[1]] == 0) + } else if (s2source == "LaSRC") { + cloudy <- which(is.na(cloudmask[[1]])) + sunny <- which(cloudmask[[1]] == 1) + } + + cloudmask[[1]][cloudy] <- 0 + cloudmask[[1]][sunny] <- 1 + cloudbin <- file.path(cloud_path, "CloudMask_Binary") + stars::write_stars(cloudmask, dsn = cloudbin, driver = "ENVI", type = "Byte", overwrite = TRUE) + cloudmasks <- list("BinaryMask" = cloudbin, "RawMask" = cloudraw) + # delete temporary file + file.remove(s2_stars$attr[whichcloud]) + if (file.exists(paste(s2_stars$attr[whichcloud], ".hdr", sep = ""))) file.remove(paste(s2_stars$attr[whichcloud], ".hdr", sep = "")) + gc() + return(cloudmasks) +} + +#" This function saves reflectance files +#" +#" @param s2_stars list. stars object containing raster data. Can be produced with function extract_from_s2_l2a +#" @param refl_path character. path for reflectance file to be stored +#" @param format character. file format for reflectance data +#" @param datatype character. data type (integer, float, 16bits, 32bits...) +#" @param s2sat character. Sentinel-2 mission ("2A" or "2B") +#" @param tile_s2 character. S2 tile name (2 numbers + 3 letters) +#" @param dateacq_s2 double. date of acquisition +#" @param MTD character. path for metadata file +#" @param MTD_MSI character. path for metadata MSI file +#" @param mtd_lasrc character. path for metadata LaSRC file +#" @param maxchunk numeric. Size of individual chunks to be written (in Mb) +#" +#" @return None +#" @importFrom stars write_stars st_apply +#" @importFrom XML xml +#" @export +save_reflectance_s2 <- function(s2_stars, refl_path, format = "ENVI", datatype = "Int16", + s2sat = NULL, tile_s2 = NULL, dateacq_s2 = NULL, + mtd = NULL, mtd_msi = NULL, mtd_lasrc = NULL, + maxchunk = 256) { + # identify if S2A or S2B, if possible + s2mission <- check_s2mission(s2sat = s2sat, tile_s2 = tile_s2, dateacq_s2 = dateacq_s2) + + # define central wavelength corresponding to each spectral band + if (s2mission == "2A") { + wl_s2 <- list("B02" = 496.6, "B03" = 560.0, "B04" = 664.5, + "B05" = 703.9, "B06" = 740.2, "B07" = 782.5, "B08" = 835.1, + "B8A" = 864.8, "B11" = 1613.7, "B12" = 2202.4) + } else if (s2mission == "2B") { + wl_s2 <- list("B02" = 492.1, "B03" = 559.0, "B04" = 665.0, + "B05" = 703.8, "B06" = 739.1, "B07" = 779.7, "B08" = 833.0, + "B8A" = 864.0, "B11" = 1610.4, "B12" = 2185.7) + } + if (s2mission == "2A") { + sensor <- "Sentinel_2A" + } else if (s2mission == "2B") { + sensor <- "Sentinel_2B" + } + + # apply offset when necessary + listbands_bis <- c("B2", "B3", "B4", "B5", "B6", "B7", "B8", "B8A", "B11", "B12") + if (!is.null(mtd_msi) && is.null(mtd_lasrc)) { + # read XML file containing info about geometry of acquisition + s2xml <- XML::xmlToList(mtd_msi) + xml_offset <- s2xml$General_Info$Product_Image_Characteristics$BOA_ADD_offset_VALUES_LIST + bands <- lapply(s2xml$General_Info$Product_Image_Characteristics$Spectral_Information_List, "[[", 4) + if (!is.null(xml_offset) && !is.null(bands)) { + bandid <- lapply(bands, "[[", 1) + bandname <- lapply(bands, "[[", 2) + offset <- data.frame("bandname" = unlist(bandname), + "bandid" = unlist(bandid), + "offset" = unlist(lapply(xml_offset, "[[", 1))) + selbands <- match(listbands_bis, offset$bandname) + offset <- offset[selbands, ] + boa_quantval <- as.numeric(s2xml$General_Info$Product_Image_Characteristics$QUANTIFICATION_VALUES_LIST$BOA_QUANTIFICATION_VALUE[1]) + } else { + offset <- data.frame("bandname" = listbands_bis, + "bandid" = c(1, 2, 3, 4, 5, 6, 7, 8, 11, 12), + "offset" = 0) + boa_quantval <- 10000 + } + } else if (!is.null(mtd_lasrc)) { + # read XML file containing info about geometry of acquisition + s2xml <- XML::xmlToList(mtd_lasrc) + attributes_lasrc <- s2xml$bands[[14]]$.attrs + attributes_lasrc_df <- data.frame(attributes_lasrc) + if (match("add_offset", rownames(attributes_lasrc_df)) > 0 && match("scale_factor", rownames(attributes_lasrc_df)) > 0) { + xml_offset <- as.numeric(attributes_lasrc[["add_offset"]]) + boa_quantval <- 1 / as.numeric(attributes_lasrc[["scale_factor"]]) + offset <- data.frame("bandname" = listbands_bis, + "bandid" = c(1, 2, 3, 4, 5, 6, 7, 8, 11, 12), + "offset" = xml_offset) + } else { + offset <- data.frame("bandname" = listbands_bis, + "bandid" = c(1, 2, 3, 4, 5, 6, 7, 8, 11, 12), + "offset" = 0) + boa_quantval <- 10000 + } + } else { + offset <- data.frame("bandname" = listbands_bis, + "bandid" = c(1, 2, 3, 4, 5, 6, 7, 8, 11, 12), + "offset" = 0) + boa_quantval <- 10000 + } + + # identify where spectral bands are in the stars object + stars_spectral <- list() + starsnames <- names(s2_stars$attr) + stars_spectral$bandname <- starsnames[which(!starsnames == "Cloud")] + stars_spectral$wavelength <- wl_s2[stars_spectral$bandname] + + sortedwl <- names(wl_s2) + reorder <- match(sortedwl, stars_spectral$bandname) + elim <- which(is.na(reorder)) + if (length(elim) > 0) { + reorder <- reorder[-elim] + } + pathr <- s2_stars$attr[reorder] + + names(pathr) <- NULL + s2_stars2 <- stars::read_stars(pathr, along = "band", proxy = TRUE) + stars_spectral$bandname <- stars_spectral$bandname[reorder] + stars_spectral$wavelength <- stars_spectral$wavelength[reorder] + + uniqueoffset <- as.numeric(unique(offset$offset)) + if (length(uniqueoffset) > 1) { + message("Warning: BOA offset differs between bands.") + message("offset will not be applied to the final S2 reflectance raster") + message("check metadata file to identify the offset applied on each band") + print(mtd_msi) + } else { + message("applying offset to reflectance data") + if (is.null(mtd_lasrc) || uniqueoffset == 0) { + offsets2 <- function(x) (round(x + uniqueoffset) * (10000 / boa_quantval)) + s2_stars2 <- stars::st_apply(X = s2_stars2, MARGIN = "band", FUN = offsets2) + } else { + offsets2 <- function(x) (round(10000 * ((x + uniqueoffset * boa_quantval) / boa_quantval))) + s2_stars2 <- stars::st_apply(X = s2_stars2, MARGIN = "band", FUN = offsets2) + } + } + write_stack_s2(stars_s2 = s2_stars2, stars_spectral = stars_spectral, refl_path = refl_path, + format = format, datatype = datatype, sensor = sensor, maxchunk = maxchunk) + # save metadata file as well if available + if (!is.null(mtd)) { + if (file.exists(mtd)) { + file.copy(from = mtd, to = file.path(dirname(refl_path), basename(mtd)), overwrite = TRUE) + } + } + # save metadata file as well if available + if (!is.null(mtd_msi)) { + if (file.exists(mtd_msi)) { + file.copy(from = mtd_msi, to = file.path(dirname(refl_path), basename(mtd_msi)), overwrite = TRUE) + } + } + # save LaSRC metadata file as well if available + if (!is.null(mtd_lasrc)) { + if (file.exists(mtd_lasrc)) { + file.copy(from = mtd_lasrc, to = file.path(dirname(refl_path), basename(mtd_lasrc)), overwrite = TRUE) + } + } + # delete temporary file + for (pathtemp in pathr) { + file.remove(pathtemp) + if (file.exists(paste(pathtemp, ".hdr", sep = ""))) file.remove(paste(pathtemp, ".hdr", sep = "")) + } + gc() + return(invisible()) +} + +#" ENVI functions +#" +#" based on https://github.com/cran/hyperSpec/blob/master/R/read.ENVI.R +#" added wavelength, fwhm, ... to header reading +#" Title +#" +#" @param x character. +#" @param separator character +#" @param trim_blank boolean. +#" +#" @return list. +#" @export +split_line <- function(x, separator, trim_blank = TRUE) { + tmp <- regexpr(separator, x) + key <- substr(x, 1, tmp - 1) + value <- substr(x, tmp + 1, nchar(x)) + if (trim_blank) { + blank_pattern <- "^[[:blank:]]*([^[:blank:]]+.*[^[:blank:]]+)[[:blank:]]*$" + key <- sub(blank_pattern, "\\1", key) + value <- sub(blank_pattern, "\\1", value) + } + value <- as.list(value) + names(value) <- key + return(value) +} + +#" save raster footprint as vector file +#" +#" @param path_raster character. path for a raster file +#" @param path_vector character. path for a vector file +#" @param driver character. driver for vector +#" +#" @return None +#" @importFrom raster raster extent projection +#" @importFrom sf st_as_sf st_write +#" @export +vectorize_raster_extent <- function(path_raster, path_vector, driver = "ESRI Shapefile") { + rast <- raster(path_raster) + e <- extent(rast) + # coerce to a SpatialPolygons object + p <- as(e, "SpatialPolygons") + projection(p) <- projection(rast) + p <- sf::st_as_sf(p) + sf::st_write(obj = p, path_vector, driver = driver) # create to a shapefile + return(invisible()) +} + +#" writes ENVI hdr file +#" +#" @param hdr content to be written +#" @param hdrpath Path of the hdr file +#" +#" @return None +#" @importFrom stringr str_count +#" @export +write_envi_header <- function(hdr, hdrpath) { + h <- lapply(hdr, function(x) { + if (length(x) > 1 || (is.character(x) && stringr::str_count(x, "\\w+") > 1)) { + x <- paste0("{", paste(x, collapse = ", "), "}") + } + # convert last numerics + x <- as.character(x) + }) + writeLines(c("ENVI", paste(names(hdr), h, sep = " = ")), con = hdrpath) + return(invisible()) +} + +#" This function writes a raster Stack object into a ENVI raster file +#" +#" @param stackobj list. raster stack +#" @param stackpath character. path where to store the stack +#" @param bands list. should include "bandname", and if possible "wavelength" +#" @param datatype character. should be INT2S or FLT4S for example +#" @param sensor character. Name of the sensor used to acquire the image +#" @param stretch boolean. Set TRUE to get 10% stretching at display for reflectance, mentioned in hdr only +#" +#" @return None +#" @importFrom utils read.table +#" @export +write_rasterstack_envi <- function(stackobj, stackpath, bands, datatype = "INT2S", + sensor = "Unknown", stretch = FALSE) { + + r <- raster::writeRaster(x = stackobj, filename = stackpath, format = "Ehdr", overwrite = TRUE, datatype = datatype) + raster::hdr(r, format = "ENVI") + # Edit hdr file to add metadata + hdr <- read_envi_header(get_hdr_name(stackpath)) + hdr$`band names` <- bands$bandname + if (length(bands$wavelength) == length(bands$bandname)) { + hdr$wavelength <- bands$wavelength + } else { + hdr$wavelength <- NULL + } + if (stretch == TRUE) { + hdr$`default stretch` <- "0.000000 1000.000000 linear" + } + hdr$`z plot range` <- NULL + hdr$`data ignore value` <- "-Inf" + hdr$`coordinate system string` <- read.table(paste(stackpath, ".prj", sep = "")) + proj <- strsplit(x = strsplit(x = projection(stackobj), split = " ")[[1]][1], split = "=")[[1]][2] + zone <- strsplit(x = strsplit(x = projection(stackobj), split = " ")[[1]][2], split = "=")[[1]][2] + datum <- strsplit(x = strsplit(x = projection(stackobj), split = " ")[[1]][3], split = "=")[[1]][2] + oldproj <- hdr$`map info` + newproj <- gsub(pattern = "projection", replacement = proj, x = oldproj) + newproj <- paste(newproj, zone, datum, sep = ", ") + hdr$`map info` <- newproj + hdr$`sensor type` <- sensor + write_envi_header(hdr = hdr, hdrpath = get_hdr_name(stackpath)) + + # remove unnecessary files + file2remove <- paste(stackpath, ".aux.xml", sep = "") + file.remove(file2remove) + file2remove <- paste(stackpath, ".prj", sep = "") + file.remove(file2remove) + file2remove <- paste(stackpath, ".stx", sep = "") + file.remove(file2remove) + return(invisible()) +} + + +#" This function writes a stars object into a raster file +#" +#" @param stars_s2 list. stars object containing raster data. Can be produced with function Crop_n_resample_S2 +#" @param stars_spectral list. band name to be saved in the stack and spectral bands corresponding to the image +#" @param refl_path character. path where to store the image +#" @param format character. default = ENVI BSQ. otherwise use gdal drivers +#" @param datatype character. should be Int16 or Float64 for example +#" @param sensor character. Name of the sensor used to acquire the image +#" @param maxchunk numeric. Size of individual chunks to be written (in Mb) +#" +#" @return None +#" @export +write_stack_s2 <- function(stars_s2, stars_spectral, refl_path, format = "ENVI", + datatype = "Int16", sensor = "Unknown", maxchunk = 256) { + + # write raster file from proxy using chunks + sizeobj <- 2 * dim(stars_s2)[1] * dim(stars_s2)[2] * dim(stars_s2)[3] / (1024**2) + nbchunks <- ceiling(sizeobj / maxchunk) + stars::write_stars(obj = stars_s2, + dsn = refl_path, + driver = format, + type = datatype, + chunk_size = c(dim(stars_s2)[1], ceiling(dim(stars_s2)[2] / nbchunks)), + progress = TRUE) + + if (format == "ENVI") { + adjust_envi_hdr(dsn = refl_path, bands = stars_spectral, + sensor = sensor, stretch = TRUE) + } + return(invisible()) +}