|
Commit message:
planemo upload for repository https://github.com/Marie59/Sentinel_2A/srs_tools commit b32737c1642aa02cc672534e42c5cb4abe0cd3e7 |
|
added:
Lib_preprocess_S2.r alpha_beta.r biodiv_indices_global.r comparison_div.r comparison_div.xml functions.r indices_spectral.r macro.xml pca_raster.r preprocess_S2.r prosail-master/R/Lib_PROSAIL.R prosail-master/R/Lib_PROSAIL_HybridInversion.R prosail-master/R/Lib_SpectralIndices.R test-data/12a0b625-9ad5-4251-a57a-305e22edef2e.xml test-data/Metadata_validation.txt test-data/Mission.csv test-data/NDVI.tabular test-data/S2A_MSIL2A_20200306T015621_N0214_R117_T51JXN_20200306T034744.zip test-data/S2A_Subset test-data/S2A_Subset.hdr test-data/S2A_T33NUD_Plots.zip val_metadata.r |
| b |
| diff -r 000000000000 -r cf69ad260611 Lib_preprocess_S2.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Lib_preprocess_S2.r Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| b'@@ -0,0 +1,1335 @@\n+# == == == == == == == == == == == == == == == == == == == == == == == == == == ==\n+# preprocS2\n+# Lib_preprocess_S2.R\n+# == == == == == == == == == == == == == == == == == == == == == == == == == == ==\n+# PROGRAMMERS:\n+# Jean-Baptiste FERET <jb.feret@teledetection.fr>\n+# Copyright 2021/08 Jean-Baptiste FERET\n+# == == == == == == == == == == == == == == == == == == == == == == == == == == ==\n+# This Library contains functions to preprocess Sentinel-2 images downloaded from\n+# different data hubs, such as THEIA, PEPS or SCIHUB\n+# == == == == == == == == == == == == == == == == == == == == == == == == == == ==\n+\n+#" This function adjusts information from ENVI header\n+#"\n+#" @param dsn character. path where to store the stack\n+#" @param bands list. should include "bandname", and if possible "wavelength"\n+#" @param sensor character. Name of the sensor used to acquire the image\n+#" @param stretch boolean. Set TRUE to get 10% stretching at display for reflectance, mentioned in hdr only\n+#"\n+#" @return None\n+#" @importFrom utils read.table\n+#" @importFrom raster hdr raster\n+#" @export\n+adjust_envi_hdr <- function(dsn, bands, sensor = "Unknown", stretch = FALSE) {\n+\n+ # Edit hdr file to add metadata\n+ hdr <- read_envi_header(get_hdr_name(dsn))\n+ hdr$`band names` <- bands$bandname\n+ if (length(bands$wavelength) == length(bands$bandname)) {\n+ hdr$wavelength <- bands$wavelength\n+ }else {\n+ hdr$wavelength <- NULL\n+ }\n+ if (stretch == TRUE) {\n+ hdr$`default stretch` <- "0.000000 1000.000000 linear"\n+ }\n+ hdr$`z plot range` <- NULL\n+ hdr$`data ignore value` <- "-Inf"\n+ hdr$`sensor type` <- sensor\n+ write_envi_header(hdr = hdr, hdrpath = get_hdr_name(dsn))\n+\n+ # remove unnecessary files\n+ file2remove <- paste(dsn, ".aux.xml", sep = "")\n+ if (file.exists(file2remove)) file.remove(file2remove)\n+ file2remove <- paste(dsn, ".prj", sep = "")\n+ if (file.exists(file2remove)) file.remove(file2remove)\n+ file2remove <- paste(dsn, ".stx", sep = "")\n+ if (file.exists(file2remove)) file.remove(file2remove)\n+ return(invisible())\n+}\n+\n+#" This function saves reflectance files\n+#"\n+#" @param s2sat character. Sentinel-2 mission ("2A" or "2B")\n+#" @param tile_s2 character. S2 tile name (2 numbers + 3 letters)\n+#" @param dateacq_s2 double. date of acquisition\n+#"\n+#" @return s2mission character. name of the S2 mission (2A or 2B)\n+#" @importFrom sen2r safe_getMetadata check_scihub_connection s2_list\n+#" @export\n+check_s2mission <- function(s2sat, tile_s2, dateacq_s2) {\n+\n+ # is mission already defined by user?\n+ if (!is.null(s2sat)) {\n+ if (s2sat == "2A") {\n+ s2mission <- "2A"\n+ }else if (s2sat == "2B") {\n+ s2mission <- "2B"\n+ }else {\n+ message("Could not identify if image from Sentinel-2A or -2B")\n+ message("Defining central wavelength of spectral bands based on S2A")\n+ s2mission <- "2A"\n+ }\n+ }else {\n+ message("Could not identify if image from Sentinel-2A or -2B")\n+ message("Defining central wavelength of spectral bands based on S2A")\n+ s2mission <- "2A"\n+ }\n+ return(s2mission)\n+}\n+\n+#" this function aims at computing directory size\n+#" @param path character. path for directory\n+#" @param recursive boolean . set T if recursive\n+#"\n+#" @return size_files numeric. size in bytes\n+#" - image stack\n+#" - path for individual band files corresponding to the stack\n+#" - path for vector (reprojected if needed)\n+#"\n+#" @importFrom raster raster\n+#" @importFrom tools file_path_sans_ext file_ext\n+#" @export\n+dir_size <- function(path, recursive = TRUE) {\n+ stopifnot(is.character(path))\n+ files <- list.files(path, full.names = TRUE, recursive = recursive)\n+ vect_size <- sapply(files, function(x) file.size(x))\n+ size_files <- sum(vect_size)\n+ return(size_files)\n+}\n+\n+#" This function reads S2 data from L2A directories downloaded from\n+#" various data hubs including THEIA, PEPS & SCIHUB (SAFE format & LaSRC)\n+#" @param path_dir_s2 character. path for S2 directory\n+#" @p'..b'", paste(x, collapse = ", "), "}")\n+ }\n+ # convert last numerics\n+ x <- as.character(x)\n+ })\n+ writeLines(c("ENVI", paste(names(hdr), h, sep = " = ")), con = hdrpath)\n+ return(invisible())\n+}\n+\n+#" This function writes a raster Stack object into a ENVI raster file\n+#"\n+#" @param stackobj list. raster stack\n+#" @param stackpath character. path where to store the stack\n+#" @param bands list. should include "bandname", and if possible "wavelength"\n+#" @param datatype character. should be INT2S or FLT4S for example\n+#" @param sensor character. Name of the sensor used to acquire the image\n+#" @param stretch boolean. Set TRUE to get 10% stretching at display for reflectance, mentioned in hdr only\n+#"\n+#" @return None\n+#" @importFrom utils read.table\n+#" @export\n+write_rasterstack_envi <- function(stackobj, stackpath, bands, datatype = "INT2S",\n+ sensor = "Unknown", stretch = FALSE) {\n+\n+ r <- raster::writeRaster(x = stackobj, filename = stackpath, format = "Ehdr", overwrite = TRUE, datatype = datatype)\n+ raster::hdr(r, format = "ENVI")\n+ # Edit hdr file to add metadata\n+ hdr <- read_envi_header(get_hdr_name(stackpath))\n+ hdr$`band names` <- bands$bandname\n+ if (length(bands$wavelength) == length(bands$bandname)) {\n+ hdr$wavelength <- bands$wavelength\n+ } else {\n+ hdr$wavelength <- NULL\n+ }\n+ if (stretch == TRUE) {\n+ hdr$`default stretch` <- "0.000000 1000.000000 linear"\n+ }\n+ hdr$`z plot range` <- NULL\n+ hdr$`data ignore value` <- "-Inf"\n+ hdr$`coordinate system string` <- read.table(paste(stackpath, ".prj", sep = ""))\n+ proj <- strsplit(x = strsplit(x = projection(stackobj), split = " ")[[1]][1], split = "=")[[1]][2]\n+ zone <- strsplit(x = strsplit(x = projection(stackobj), split = " ")[[1]][2], split = "=")[[1]][2]\n+ datum <- strsplit(x = strsplit(x = projection(stackobj), split = " ")[[1]][3], split = "=")[[1]][2]\n+ oldproj <- hdr$`map info`\n+ newproj <- gsub(pattern = "projection", replacement = proj, x = oldproj)\n+ newproj <- paste(newproj, zone, datum, sep = ", ")\n+ hdr$`map info` <- newproj\n+ hdr$`sensor type` <- sensor\n+ write_envi_header(hdr = hdr, hdrpath = get_hdr_name(stackpath))\n+\n+ # remove unnecessary files\n+ file2remove <- paste(stackpath, ".aux.xml", sep = "")\n+ file.remove(file2remove)\n+ file2remove <- paste(stackpath, ".prj", sep = "")\n+ file.remove(file2remove)\n+ file2remove <- paste(stackpath, ".stx", sep = "")\n+ file.remove(file2remove)\n+ return(invisible())\n+}\n+\n+\n+#" This function writes a stars object into a raster file\n+#"\n+#" @param stars_s2 list. stars object containing raster data. Can be produced with function Crop_n_resample_S2\n+#" @param stars_spectral list. band name to be saved in the stack and spectral bands corresponding to the image\n+#" @param refl_path character. path where to store the image\n+#" @param format character. default = ENVI BSQ. otherwise use gdal drivers\n+#" @param datatype character. should be Int16 or Float64 for example\n+#" @param sensor character. Name of the sensor used to acquire the image\n+#" @param maxchunk numeric. Size of individual chunks to be written (in Mb)\n+#"\n+#" @return None\n+#" @export\n+write_stack_s2 <- function(stars_s2, stars_spectral, refl_path, format = "ENVI",\n+ datatype = "Int16", sensor = "Unknown", maxchunk = 256) {\n+\n+ # write raster file from proxy using chunks\n+ sizeobj <- 2 * dim(stars_s2)[1] * dim(stars_s2)[2] * dim(stars_s2)[3] / (1024**2)\n+ nbchunks <- ceiling(sizeobj / maxchunk)\n+ stars::write_stars(obj = stars_s2,\n+ dsn = refl_path,\n+ driver = format,\n+ type = datatype,\n+ chunk_size = c(dim(stars_s2)[1], ceiling(dim(stars_s2)[2] / nbchunks)),\n+ progress = TRUE)\n+\n+ if (format == "ENVI") {\n+ adjust_envi_hdr(dsn = refl_path, bands = stars_spectral,\n+ sensor = sensor, stretch = TRUE)\n+ }\n+ return(invisible())\n+}\n' |
| b |
| diff -r 000000000000 -r cf69ad260611 alpha_beta.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/alpha_beta.r Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| @@ -0,0 +1,142 @@ +#Rscript + +########################################### +## Mapping alpha and beta diversity ## +########################################### + +#####Packages : stars +# utils +# biodivmapr +# raster +# sf +# mapview +# leafpop +# RColorBrewer +# labdsv +# rgdal +# ggplot2 +# gridExtra + +#####Load arguments + +args <- commandArgs(trailingOnly = TRUE) + +#####Import the S2 data + +if (length(args) < 1) { + stop("This tool needs at least 1 argument") +}else { + data_raster <- args[1] + rasterheader <- args[2] + data <- args[3] + # type of PCA: + # PCA: no rescaling of the data + # SPCA: rescaling of the data + typepca <- as.character(args[4]) + alpha <- as.logical(args[5]) + beta <- as.logical(args[6]) + funct <- as.logical(args[7]) + all <- as.logical(args[8]) + nbcpu <- as.integer(args[9]) + source(args[10]) +} + +################################################################################ +## DEFINE PARAMETERS FOR DATASET TO BE PROCESSED ## +################################################################################ +if (data_raster == "") { + #Create a directory where to unzip your folder of data + dir.create("data_dir") + unzip(data, exdir = "data_dir") + # Path to raster + data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl") + input_image_file <- file.path("data_dir/results/Reflectance", data_raster[1]) + input_header_file <- file.path("data_dir/results/Reflectance", data_raster[2]) + +} else { + input_image_file <- file.path(getwd(), data_raster, fsep = "/") + input_header_file <- file.path(getwd(), rasterheader, fsep = "/") +} + +################################################################################ +## PROCESS IMAGE ## +################################################################################ +# 1- Filter data in order to discard non vegetated / shaded / cloudy pixels + +print("PERFORM PCA ON RASTER") +pca_output <- biodivMapR::perform_PCA(Input_Image_File = input_image_file, Input_Mask_File = input_mask_file, + Output_Dir = output_dir, TypePCA = typepca, FilterPCA = filterpca, nbCPU = nbcpu, MaxRAM = maxram) + +pca_files <- pca_output$PCA_Files +pix_per_partition <- pca_output$Pix_Per_Partition +nb_partitions <- pca_output$nb_partitions +# path for the updated mask +input_mask_file <- pca_output$MaskPath + + +selected_pcs <- seq(1, dim(raster::stack(input_image_file))[3]) + +selected_pcs <- all(selected_pcs) +################################################################################ +## MAP ALPHA AND BETA DIVERSITY ## +################################################################################ +print("MAP SPECTRAL SPECIES") + +kmeans_info <- biodivMapR::map_spectral_species(Input_Image_File = input_image_file, Output_Dir = output_dir, PCA_Files = pca_files, Input_Mask_File = input_mask_file, SelectedPCs = selected_pcs, Pix_Per_Partition = pix_per_partition, nb_partitions = nb_partitions, TypePCA = typepca, nbCPU = nbcpu, MaxRAM = maxram, nbclusters = nbclusters) + +image_name <- tools::file_path_sans_ext(basename(input_image_file)) +if (alpha == TRUE || beta == TRUE || all == TRUE) { +## alpha + print("MAP ALPHA DIVERSITY") + index_alpha <- c("Shannon") + alpha_div <- biodivMapR::map_alpha_div(Input_Image_File = input_image_file, Output_Dir = output_dir, TypePCA = typepca, window_size = window_size, nbCPU = nbcpu, MaxRAM = maxram, Index_Alpha = index_alpha, nbclusters = nbclusters, FullRes = TRUE, LowRes = FALSE, MapSTD = FALSE) + + alpha_zip <- file.path(output_dir, image_name, typepca, "ALPHA", "Shannon_10_Fullres.zip") + alpha_path <- file.path(output_dir, image_name, typepca, "ALPHA") + unzip(alpha_zip, exdir = alpha_path) + alpha_path <- file.path(output_dir, image_name, typepca, "ALPHA", "Shannon_10_Fullres") + alpha_raster <- raster::raster(alpha_path) + get_alpha <- convert_raster(alpha_raster) + + if (alpha == TRUE || all == TRUE) { + colnames(get_alpha) <- c("Alpha", "longitude", "latitude") + plot_indices(get_alpha, titre = "Alpha") + + write.table(get_alpha, file = "alpha.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE) +} + if (beta == TRUE || all == TRUE) { +## beta + print("MAP BETA DIVERSITY") + beta_div <- biodivMapR::map_beta_div(Input_Image_File = input_image_file, Output_Dir = output_dir, TypePCA = typepca, window_size = window_size, nb_partitions = nb_partitions, nbCPU = nbcpu, MaxRAM = maxram, nbclusters = nbclusters) + + beta_path <- file.path(output_dir, image_name, typepca, "BETA", "BetaDiversity_BCdiss_PCO_10") + beta_raster <- raster::raster(beta_path) + get_beta <- convert_raster(beta_raster) + + colnames(get_beta) <- c("Beta", "longitude", "latitude") + plot_indices(get_beta, titre = "Beta") + + write.table(get_beta, file = "beta.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE) + } +} + + +################################################################################ +## COMPUTE ALPHA AND BETA DIVERSITY FROM FIELD PLOTS ## +################################################################################ + +if (funct == TRUE || all == TRUE) { + mapper <- biodivMapR::map_functional_div(Original_Image_File = input_image_file, Functional_File = pca_files, Selected_Features = selected_pcs, Output_Dir = output_dir, window_size = window_size, nbCPU = nbcpu, MaxRAM = maxram, TypePCA = typepca, FullRes = TRUE, LowRes = FALSE, MapSTD = FALSE) + + funct_zip <- file.path(output_dir, image_name, typepca, "FUNCTIONAL", "FunctionalDiversity_Map_MeanFilter_Fullres.zip") + funct_path <- file.path(output_dir, image_name, typepca, "FUNCTIONAL") + unzip(funct_zip, exdir = funct_path) + funct_path <- file.path(output_dir, image_name, typepca, "FUNCTIONAL", "FunctionalDiversity_Map_MeanFilter_Fullres") + funct_raster <- raster::raster(funct_path) + get_funct <- convert_raster(funct_raster) + + colnames(get_funct) <- c("Functionnal", "longitude", "latitude") + plot_indices(get_funct, titre = "Functionnal") + + write.table(get_funct, file = "Functionnal.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE) +} |
| b |
| diff -r 000000000000 -r cf69ad260611 biodiv_indices_global.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/biodiv_indices_global.r Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| @@ -0,0 +1,112 @@ +#Rscript + +########################################### +## Mapping biodiversity indicators ## +########################################### + +#####Packages : raster, +# rgdal, +# sp, +# rasterdiv, +# ggplot2, + +#####Load arguments + +args <- commandArgs(trailingOnly = TRUE) + +#####Import the S2 data + +if (length(args) < 1) { + stop("This tool needs at least 1 argument") +}else { + data_raster <- args[1] + data_header <- args[2] + data <- args[3] + alpha <- as.numeric(args[4]) + source(args[5]) + +} + +######################################################################## +## COMPUTE BIODIVERSITY INDICES ## +######################################################################## + +if (data_raster == "") { + #Create a directory where to unzip your folder of data + dir.create("data_dir") + unzip(data, exdir = "data_dir") + # Path to raster + data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl") + data_raster <- file.path("data_dir/results/Reflectance", data_raster[1]) +} + +# Read raster +copndvi <- raster::raster(data_raster) + +copndvilr <- raster::reclassify(copndvi, cbind(252, 255, NA), right = TRUE) + +#Resample using raster::aggregate and a linear factor of 10 +copndvilr <- raster::aggregate(copndvilr, fact = 20) +#Set float numbers as integers to further speed up the calculation +storage.mode(copndvilr[]) <- "integer" + +# Convert raster to SpatialPointsDataFrame +r_pts <- convert_raster(copndvilr) + +#Shannon's Diversity +sha <- rasterdiv::Shannon(copndvilr, window = 9, np = 1) +sha_df <- data.frame(raster::rasterToPoints(sha, spatial = TRUE)) +sha_name <- "Shannon" +r_pts[, sha_name] <- sha_df[, 1] + +#Renyi's Index +ren <- rasterdiv::Renyi(copndvilr, window = 9, alpha = alpha, np = 1) +ren_df <- data.frame(raster::rasterToPoints(ren[[1]])) +ren_name <- "Renyi" +r_pts[, ren_name] <- ren_df[, 3] + +#Berger-Parker's Index +ber <- rasterdiv::BergerParker(copndvilr, window = 9, np = 1) +ber_df <- data.frame(raster::rasterToPoints(ber, spatial = TRUE)) +ber_name <- "Berger-Parker" +r_pts[, ber_name] <- ber_df[, 1] + +#Pielou's Evenness +pie <- rasterdiv::Pielou(copndvilr, window = 9, np = 1) +pie_df <- data.frame(raster::rasterToPoints(pie)) +if (length(pie_df[, 1]) == length(r_pts[, 1])) { + pie_name <- "Pielou" + r_pts[, pie_name] <- pie_df[, 1] +} + +#Hill's numbers +hil <- rasterdiv::Hill(copndvilr, window = 9, alpha = alpha, np = 1) +hil_df <- data.frame(raster::rasterToPoints(hil[[1]])) +hil_name <- "Hill" +r_pts[, hil_name] <- hil_df[, 3] + +#Parametric Rao's quadratic entropy with alpha ranging from 1 to 5 +prao <- rasterdiv::paRao(copndvilr, window = 9, alpha = alpha, dist_m = "euclidean", np = 1) +prao_df <- data.frame(raster::rasterToPoints(prao$window.9[[1]])) +prao_name <- "Prao" +r_pts[, prao_name] <- prao_df[, 3] + +#Cumulative Residual Entropy +cre <- rasterdiv::CRE(copndvilr, window = 9, np = 1) +cre_df <- data.frame(raster::rasterToPoints(cre)) +if (length(cre_df[, 1]) == length(r_pts[, 1])) { + cre_name <- "CRE" + r_pts[, cre_name] <- cre_df[, 1] +} + +if (length(cre_df[, 1]) == length(r_pts[, 1]) || length(pie_df[, 1]) == length(r_pts[, 1])) { +list_indice <- list("Shannon", "Renyi", "Berger-Parker", "Pielou", "Hill", "Prao", "CRE") +} else { +list_indice <- list("Shannon", "Renyi", "Berger-Parker", "Hill", "Prao") +} +## Plotting all the graph and writing a tabular +for (indice in list_indice) { + plot_indices(data = r_pts, titre = indice) +} + +write.table(r_pts, file = "BiodivIndex.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE) |
| b |
| diff -r 000000000000 -r cf69ad260611 comparison_div.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/comparison_div.r Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| b'@@ -0,0 +1,190 @@\n+#Rscript\n+\n+###########################################\n+## Mapping alpha and beta diversity ##\n+###########################################\n+\n+#####Packages : stars\n+# utils\n+# biodivmapr\n+# raster\n+# sf\n+# mapview\n+# leafpop\n+# RColorBrewer\n+# labdsv\n+# rgdal\n+# ggplot2\n+# gridExtra\n+##remotes::install_github("jbferet/biodivMapR")\n+#####Load arguments\n+\n+args <- commandArgs(trailingOnly = TRUE)\n+\n+#####Import the S2 data\n+\n+if (length(args) < 1) {\n+ stop("This tool needs at least 1 argument")\n+}else {\n+ data_raster <- args[1]\n+ rasterheader <- args[2]\n+ data <- args[3]\n+ plots_zip <- args[4]\n+ choice <- as.character(args[5])\n+ source(args[6])\n+ # type of PCA:\n+ # PCA: no rescaling of the data\n+ # SPCA: rescaling of the data\n+ typepca <- as.character(args[7])\n+}\n+\n+################################################################################\n+## DEFINE PARAMETERS FOR DATASET TO BE PROCESSED ##\n+################################################################################\n+if (data_raster == "") {\n+ #Create a directory where to unzip your folder of data\n+ dir.create("data_dir")\n+ unzip(data, exdir = "data_dir")\n+ # Path to raster\n+ data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl")\n+ input_image_file <- file.path("data_dir/results/Reflectance", data_raster[1])\n+ input_header_file <- file.path("data_dir/results/Reflectance", data_raster[2])\n+\n+} else {\n+ input_image_file <- file.path(getwd(), data_raster, fsep = "/")\n+ input_header_file <- file.path(getwd(), rasterheader, fsep = "/")\n+}\n+\n+################################################################################\n+## PROCESS IMAGE ##\n+################################################################################\n+# 1- Filter data in order to discard non vegetated / shaded / cloudy pixels\n+\n+print("PERFORM PCA ON RASTER")\n+pca_output <- biodivMapR::perform_PCA(Input_Image_File = input_image_file, Input_Mask_File = input_mask_file,\n+ Output_Dir = output_dir, TypePCA = typepca, FilterPCA = filterpca, nbCPU = nbcpu, MaxRAM = maxram)\n+\n+pca_files <- pca_output$PCA_Files\n+pix_per_partition <- pca_output$Pix_Per_Partition\n+nb_partitions <- pca_output$nb_partitions\n+# path for the updated mask\n+input_mask_file <- pca_output$MaskPath\n+\n+# 3- Select principal components from the PCA raster\n+# Select components from the PCA/SPCA/MNF raster\n+sel_compo <- c("1\\n", "2\\n", "3\\n", "4\\n", "5\\n", "6\\n", "7\\n", "8")\n+image_name <- tools::file_path_sans_ext(basename(input_image_file))\n+output_dir_full <- file.path(output_dir, image_name, typepca, "PCA")\n+\n+write.table(sel_compo, paste0(output_dir_full, "/Selected_Components.txt"))\n+sel_pc <- file.path(output_dir_full, "Selected_Components.txt")\n+\n+\n+################################################################################\n+## MAP ALPHA AND BETA DIVERSITY ##\n+################################################################################\n+print("MAP SPECTRAL SPECIES")\n+\n+kmeans_info <- biodivMapR::map_spectral_species(Input_Image_File = input_image_file, Output_Dir = output_dir, PCA_Files = pca_files, Input_Mask_File = input_mask_file, Pix_Per_Partition = pix_per_partition, nb_partitions = nb_partitions, nbCPU = nbcpu, MaxRAM = maxram, nbclusters = nbclusters, TypePCA = typepca)\n+\n+################################################################################\n+## COMPUTE ALPHA AND BETA DIVERSITY FROM FIELD PLOTS ##\n+################################################################################\n+## read selected features from dimensionality reduction\n+\n+## path for selected components\n+\n+# location'..b'div <- c(biodiv_indicators$FunctionalDiversity$FDiv)\n+# if no name for plots\n+biodiv_indicators$Name_Plot <- seq(1, length(biodiv_indicators$Shannon[[1]]), by = 1)\n+\n+\n+####################################################\n+# write RS indicators #\n+####################################################\n+# write a table for Shannon index\n+\n+# write a table for all spectral diversity indices corresponding to alpha diversity\n+results <- data.frame(name_vector, biodiv_indicators$Richness, biodiv_indicators$Fisher,\n+ biodiv_indicators$Shannon, biodiv_indicators$Simpson,\n+ biodiv_indicators$FunctionalDiversity$FRic,\n+ biodiv_indicators$FunctionalDiversity$FEve,\n+ biodiv_indicators$FunctionalDiversity$FDiv)\n+\n+names(results) <- c("ID_Plot", "Species_Richness", "Fisher", "Shannon", "Simpson", "fric", "feve", "fdiv")\n+write.table(results, file = "Diversity.tabular", sep = "\\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)\n+\n+if (choice == "Y") {\n+# write a table for Bray Curtis dissimilarity\n+bc_mean <- biodiv_indicators$BCdiss\n+bray_curtis <- data.frame(name_vector, bc_mean)\n+colnames(bray_curtis) <- c("ID_Plot", bray_curtis[, 1])\n+write.table(bray_curtis, file = "BrayCurtis.tabular", sep = "\\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)\n+\n+####################################################\n+# illustrate results\n+####################################################\n+# apply ordination using PCoA (same as done for map_beta_div)\n+\n+mat_bc_dist <- as.dist(bc_mean, diag = FALSE, upper = FALSE)\n+betapco <- labdsv::pco(mat_bc_dist, k = 3)\n+\n+# assign a type of vegetation to each plot, assuming that the type of vegetation\n+# is defined by the name of the shapefile\n+\n+nbsamples <- shpname <- c()\n+for (i in 1:length(path_vector)) {\n+ shp <- path_vector[i]\n+ nbsamples[i] <- length(rgdal::readOGR(shp, verbose = FALSE))\n+ shpname[i] <- tools::file_path_sans_ext(basename(shp))\n+}\n+\n+type_vegetation <- c()\n+for (i in 1: length(nbsamples)) {\n+ for (j in 1:nbsamples[i]) {\n+ type_vegetation <- c(type_vegetation, shpname[i])\n+ }\n+}\n+\n+#data frame including a selection of alpha diversity metrics and beta diversity expressed as coordinates in the PCoA space\n+results <- data.frame("vgtype" = type_vegetation, "pco1" = betapco$points[, 1], "pco2" = betapco$points[, 2], "pco3" = betapco$points[, 3], "shannon" = shannon_rs, "fric" = fric, "feve" = feve, "fdiv" = fdiv)\n+\n+#plot field data in the PCoA space, with size corresponding to shannon index\n+g1 <- ggplot2::ggplot(results, ggplot2::aes(x = pco1, y = pco2, color = vgtype, size = shannon)) + ggplot2::geom_point(alpha = 0.6) + ggplot2::scale_color_manual(values = c("#e6140a", "#e6d214", "#e68214", "#145ae6"))\n+\n+g2 <- ggplot2::ggplot(results, ggplot2::aes(x = pco1, y = pco3, color = vgtype, size = shannon)) + ggplot2::geom_point(alpha = 0.6) + ggplot2::scale_color_manual(values = c("#e6140a", "#e6d214", "#e68214", "#145ae6"))\n+\n+g3 <- ggplot2::ggplot(results, ggplot2::aes(x = pco2, y = pco3, color = vgtype, size = shannon)) + ggplot2::geom_point(alpha = 0.6) + ggplot2::scale_color_manual(values = c("#e6140a", "#e6d214", "#e68214", "#145ae6"))\n+\n+#extract legend\n+get_legend <- function(a_gplot) {\n+ tmp <- ggplot2::ggplot_gtable(ggplot2::ggplot_build(a_gplot))\n+ leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")\n+ legend <- tmp$grobs[[leg]]\n+ return(legend)\n+}\n+\n+legend <- get_legend(g3)\n+gall <- gridExtra::grid.arrange(gridExtra::arrangeGrob(g1 + ggplot2::theme(legend.position = "none"), g2 + ggplot2::theme(legend.position = "none"), g3 + ggplot2::theme(legend.position = "none"), nrow = 1), legend, nrow = 2, heights = c(3, 2))\n+\n+\n+filename <- ggplot2::ggsave("BetaDiversity_PcoA1_vs_PcoA2_vs_PcoA3.png", gall, scale = 0.65, width = 12, height = 9, units = "in", dpi = 200, limitsize = TRUE)\n+\n+filename\n+}\n' |
| b |
| diff -r 000000000000 -r cf69ad260611 comparison_div.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/comparison_div.xml Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| @@ -0,0 +1,138 @@ +<tool id="srs_process_data" name="Compare diversity indicators" version="@VERSION@" profile="20.01"> + <description>with remote sensing data</description> + <macros> + <import>macro.xml</import> + </macros> + <expand macro="SRS_requirements"/> + <command detect_errors="exit_code"><![CDATA[ + #import re + #if $method.type == 'envi_bil': + #set input_raster = $method.input_raster + #set input_raster_identifier = re.sub('[^\s\w\-]', '_', str($input_raster.element_identifier)) + #set input_header = $method.input_header + #set input_header_identifier = re.sub('[^\s\w\-]+[^.hdr]', '_', str($input_header.element_identifier)) + cp '${input_raster}' '${input_raster_identifier}' && + cp '${input_header}' '${input_header_identifier}' && + #end if + Rscript + '$__tool_directory__/comparison_div.r' + #if $method.type == 'envi_bil': + '$input_raster_identifier' + '$input_header_identifier' + '' + #else: + '' + '' + '$method.input' + #end if + '$input_plot' + '$choice' + '$__tool_directory__/functions.r' + '$typepca' + '$output_div' + '$output_beta' + '$plots' + ]]> + </command> + <inputs> + <conditional name="method"> + <param name="type" type="select" label="In which format are your data ?"> + <option value="zipper">The data you are using are in a zip folder Reflectance</option> + <option value="envi_bil">Your already have the files in ENVI BIL format</option> + </param> + <when value="zipper"> + <param name="input" type="data" format="zip" multiple="true" label="Input data"/> + </when> + <when value="envi_bil"> + <param name="input_raster" type="data" format="bil" label="Input raster" help="It can be the raw data in bil or the PCA raster layer in bil"/> + <param name="input_header" type="data" format="hdr" label="Input header"/> + </when> + </conditional> + <param name="input_plot" type="data" format="data" label="Plots folder zip"/> + <param name="choice" type="select" label="Do you want to compute beta diversity (needs mutliple locations) ?" display="radio"> + <option value="Y">Yes</option> + <option value="N">No</option> + </param> + <param name="typepca" type="select" label="Do you want to do a PCA or a SPCA ?" display="radio" help="If you choose PCA there is no rescaling of the data as oppposed as if you choose SPCA"> + <option value="SPCA">SPCA</option> + <option value="PCA">PCA</option> + </param> + </inputs> + <outputs> + <data name="output_div" from_work_dir="Diversity.tabular" format="tabular" label="Global diversity"> + </data> + <data name="output_beta" from_work_dir="BrayCurtis.tabular" format="tabular" label="Bray Curtis"> + <filter>choice == 'Y'</filter> + </data> + <collection type="list" name="plots" label="Comparison plot"> + <discover_datasets pattern="(?P<designation>.+)\.png" visible="false" format="png"/> + <filter>choice =='Y'</filter> + </collection> + </outputs> + <tests> + <test> + <param name="type" value="envi_bil"/> + <param name="input_raster" value="S2A_Subset"/> + <param name="input_header" value="S2A_Subset.hdr"/> + <param name="input_plot" value="S2A_T33NUD_Plots.zip"/> + <param name="choice" value="Y"/> + <output name="output_div"> + <assert_contents> + <has_n_lines n="25"/> + </assert_contents> + </output> + <output name="output_beta"> + <assert_contents> + <has_n_lines n="25"/> + </assert_contents> + </output> + <output_collection name="plots" type="list" count="1"/> + </test> + </tests> + <help><![CDATA[ +======================================================================== +Process satellite remote sensing data to produce biodiversity indicators +======================================================================== + + +**What it does** + +Féret and Asner (2014) developed a method for **tropical forest** diversity mapping based on very high spatial resolution airborne imaging spectroscopy. + +The goal of this tool using the package biodivMapR is to compute diversity indices over each spatial polygon of a shapefile of plots, if available, in order to compare field inventories with diversity indices estimated from remotely-sensed images. + +**Input description** + +It expects an image file as input, with a specific data format. ENVI HDR image with BIL interleave required. +The image is an ENVI raster including : + +- A binary file (which has no extension here). + +- A header file (with .hdr extension). + +The header file is a text file including all necessary metadata which can be read with a text editor. It includes image dimensions, projection, and the name and central wavelength for each spectral band. + +In order to get such input we advise to use the tool preprocessing sentinel 2 data. + ++--------------+----------+--------------+ +| BIL | ENVI HDR | Shapefiles | ++==============+==========+==============+ +| raster stack | Metadata | plots.zip | ++--------------+----------+--------------+ +| ... | ... | ... | ++--------------+----------+--------------+ + +**Output** + +- Two tabulars : + - One matrix for Bray-Curtis indicator + - One table for the following indicators; Species richness, shannon, fisher, simpson, richness, eveness, divergence + +- One comparison png plot in the Pcoa space that summarizes α- and β-diversity in scatterplots and illustrates that the combination of the three components computed with PCoA allows proper differentiation among vegetation types: + - PCoA#1 allows differentiating medium and high diversity forests from low diversity forest and low vegetation, but does not discriminate medium and high diversity forests. + - PCoA#2 allows differentiating low diversity forest from medium/high diversity forests and low vegetation + - PCoA#3 allows differentiating medium diversity forests from high diversity forests and low vegetation. + + ]]> </help> + <expand macro="SRS_BDMRref"/> +</tool> |
| b |
| diff -r 000000000000 -r cf69ad260611 functions.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/functions.r Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| @@ -0,0 +1,79 @@ +#Rscript + +######################### +## Functions ## +######################### + +#####Packages : raster +# sp +# ggplot2 + +####Set paramaters for all tools using BiodivMapR + +# path for the Mask raster corresponding to image to process +# expected to be in ENVI HDR format, 1 band, integer 8bits +# expected values in the raster: 0 = masked, 1 = selected +# set to FALSE if no mask available +input_mask_file <- FALSE + +# relative or absolute path for the Directory where results will be stored +# For each image processed, a subdirectory will be created after its name +output_dir <- "RESULTS" + +# SPATIAL RESOLUTION +# resolution of spatial units for alpha and beta diversity maps (in pixels), relative to original image +# if Res.Map = 10 for images with 10 m spatial resolution, then spatial units will be 10 pixels x 10m = 100m x 100m surfaces +# rule of thumb: spatial units between 0.25 and 4 ha usually match with ground data +# too small window_size results in low number of pixels per spatial unit, hence limited range of variation of diversity in the image +window_size <- 10 + +# PCA FILTERING: Set to TRUE if you want second filtering based on PCA outliers to be processed. Slower +filterpca <- TRUE + +################################################################################ +## DEFINE PARAMETERS FOR METHOD ## +################################################################################ +nbcpu <- 4 +maxram <- 0.5 +nbclusters <- 50 + +################################################################################ +## PROCESS IMAGE ## +################################################################################ +# 1- Filter data in order to discard non vegetated / shaded / cloudy pixels +ndvi_thresh <- 0.5 +blue_thresh <- 500 +nir_thresh <- 1500 +continuum_removal <- TRUE + + + +#### Convert raster to dataframe + +# Convert raster to SpatialPointsDataFrame +convert_raster <- function(data_raster) { +r_pts <- raster::rasterToPoints(data_raster, spatial = TRUE) + +# reproject sp object +geo_prj <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" +r_pts <- sp::spTransform(r_pts, sp::CRS(geo_prj)) + + +# Assign coordinates to @data slot, display first 6 rows of data.frame +r_pts@data <- data.frame(r_pts@data, longitude = sp::coordinates(r_pts)[, 1], + latitude = sp::coordinates(r_pts)[, 2]) + +return(r_pts@data) +} + + +#### Potting indices + +plot_indices <- function(data, titre) { + graph_indices <- ggplot2::ggplot(data) + + ggplot2::geom_point(ggplot2::aes_string(x = data[, 2], y = data[, 3], color = data[, titre]), shape = "square", size = 2) + ggplot2::scale_colour_gradient(low = "blue", high = "orange", na.value = "grey50") + + ggplot2::xlab("Longitude") + ggplot2::ylab("Latitude") + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1), plot.title = ggplot2::element_text(color = "black", size = 12, face = "bold")) + ggplot2::ggtitle(titre) + +ggplot2::ggsave(paste0(titre, ".png"), graph_indices, width = 12, height = 10, units = "cm") +} |
| b |
| diff -r 000000000000 -r cf69ad260611 indices_spectral.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/indices_spectral.r Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| @@ -0,0 +1,105 @@ +#Rscript + +########################################### +## Mapping alpha and beta diversity ## +########################################### + +#####Packages : expint, +# pracma, +# R.utils, +# raster, +# sp, +# matrixStats, +# ggplot2, +# expandFunctions, +# stringr, +# XML, +# rgdal, +# stars, +#####Load arguments + +args <- commandArgs(trailingOnly = TRUE) + +#####Import the S2 data + +if (length(args) < 1) { + stop("This tool needs at least 1 argument") +}else { + data_raster <- args[1] + data_header <- args[2] + data <- args[3] + source(args[4]) + source(args[5]) + source(args[6]) + indice_choice <- strsplit(args[7], ",")[[1]] + source(args[8]) + output_raster <- as.character(args[9]) + +} + +######################################################################## +## COMPUTE SPECTRAL INDEX : NDVI ## +######################################################################## + +if (data != "") { + #Create a directory where to unzip your folder of data + dir.create("data_dir") + unzip(data, exdir = "data_dir") + + # Read raster + data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl") + data_raster <- file.path("data_dir/results/Reflectance", data_raster[1]) + refl <- raster::brick(data_raster) + refl2 <- raster::raster(data_raster) +} else { + # Read raster + refl <- raster::brick(data_raster) + refl2 <- raster::raster(data_raster) +} +# get raster band name and clean format. Expecting band name and wav +# reflFactor = 10000 when reflectance is coded as INT16 +refl <- raster::aggregate(refl, fact = 10) + +# Convert raster to SpatialPointsDataFrame +refl2 <- raster::aggregate(refl2, fact = 10) +r_pts <- convert_raster(refl2) +table_ind <- r_pts +# create directory for Spectralelength to be documented in image +hdr_refl <- read_envi_header(get_hdr_name(data_raster)) +sensorbands <- hdr_refl$wavelength +# compute a set of spectral indices defined by indexlist from S2 data indices +si_path <- file.path("SpectralIndices") +dir.create(path = si_path, showWarnings = FALSE, recursive = TRUE) +# Save spectral indices + +indices <- lapply(indice_choice, function(x) { + indices_list <- computespectralindices_raster(refl = refl, sensorbands = sensorbands, + sel_indices = x, + reflfactor = 10000, stackout = FALSE) + + index_path <- file.path(si_path, paste(basename(data_raster), "_", x, sep = "")) + spec_indices <- stars::write_stars(stars::st_as_stars(indices_list$spectralindices[[1]]), dsn = index_path, driver = "ENVI", type = "Float32") + + # write band name in HDR + hdr <- read_envi_header(get_hdr_name(index_path)) + hdr$`band names` <- x + write_envi_header(hdr = hdr, hdrpath = get_hdr_name(index_path)) + # Writting the tabular and the plot + r_pts[, x] <- as.data.frame(sapply(spec_indices, c)) + plot_indices(data = r_pts, titre = x) + return(r_pts) +}) + +new_table <- as.data.frame(indices) +new_table <- new_table[, !grepl("longitude", names(new_table))] +new_table <- new_table[, !grepl("latitude", names(new_table))] +new_table <- new_table[, !grepl(basename(data_raster), names(new_table))] + +table_ind <- cbind(table_ind, new_table) +if (length(indice_choice) == 1) { + colnames(table_ind) <- c(basename(data_raster), "longitude", "latitude", indice_choice) +} + +write.table(table_ind, file = "Spec_Index.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE) + +# Get the raster layer of the indice as an output |
| b |
| diff -r 000000000000 -r cf69ad260611 macro.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/macro.xml Mon Jan 09 13:36:02 2023 +0000 |
| b |
| @@ -0,0 +1,102 @@ +<macros> + <token name="@VERSION@">0.0.1</token> + <xml name="SRS_requirements"> + <requirements> + <requirement type="package" version="4.2.2">r-base</requirement> + <requirement type="package" version="2.12.2">r-r.utils</requirement> + <requirement type="package" version="1.0_7">r-sf</requirement> + <requirement type="package" version="1.5_32">r-rgdal</requirement> + <requirement type="package" version="0.6_1">r-rgeos</requirement> + <requirement type="package" version="0.5_5">r-stars</requirement> + <requirement type="package" version="1.5.0">r-stringr</requirement> + <requirement type="package" version="1.9.8">r-biodivmapr</requirement> + <requirement type="package" version="0.36">r-xfun</requirement> + <requirement type="package" version="0.51.4">r-rastervis</requirement> + <requirement type="package" version="3.4.0">r-ggplot2</requirement> + <requirement type="package" version="2.11.0">r-mapview</requirement> + <requirement type="package" version="2.3">r-gridextra</requirement> + <requirement type="package" version="3.0.0">r-emstreer</requirement> + <requirement type="package" version="3.2.3">r-filesstrings</requirement> + <!--requirement type="package" version="9.0.1">proj</requirement--> + <yield/> + </requirements> + </xml> + <xml name="SRS_input"> + <param name="input" type="data" format="tabular,csv,txt" label="Input table"/> + <param name="colnames" type="boolean" label="First line is a header line" checked="true"/> + </xml> + <xml name="SRS_bibref"> + <citations> + <citation type="bibtex"> + @Manual{, + title = {obisindicators: Develop marine biodiversity indicators from OBIS data}, + author = {Ben Ben and Pieter Provoost and Tylar Murray}, + year = {2022}, + note = {https://marinebon.org/obisindicators, + https://github.com/marinebon/obisindicators}, + } + </citation> + </citations> + </xml> + <xml name="SRS_S2ref"> + <citations> + <citation type="bibtex"> + @Manual{, + title = {preprocS2: preprocessing of Sentinel-2 data downloaded from various data hubs / produced from various atmospheric correction methods}, + author = {Jean-Baptiste Feret}, + year = {2022}, + note = {R package version 1.1.3}, + url = {https://gitlab.com/jbferet/preprocS2}, + } + </citation> + </citations> + </xml> + <xml name="SRS_BDMRref"> + <citations> + <citation type="bibtex"> + @Manual{, + title = {biodivMapR: biodivMapR: an R package for a- and ß-diversity mapping using remotely-sensed images}, + author = {Jean-Baptiste Feret and Florian {de Boissieu}}, + year = {2022}, + note = {R package version 1.9.4}, + url = {https://github.com/jbferet/biodivMapR}, + } + </citation> + </citations> + </xml> + <xml name="SRS_rasterdivref"> + <citations> + <citation type="doi">doi:10.1016/j.ecolind.2016.07.039</citation> + <citation type="doi">doi:10.1101/2021.01.23.427872</citation> + </citations> + </xml> + <xml name="SRS_prosailref"> + <citations> + <citation type="bibtex"> + @Manual{, + title = {prosail: PROSAIL leaf and canopy radiative transfer model and inversion routines}, + author = {Jean-Baptiste Feret and Florian {de Boissieu}}, + year = {2022}, + note = {R package version 1.1.1}, + url = {https://gitlab.com/jbferet/prosail}, + } + </citation> + </citations> + </xml> + <xml name="SRS_metaref"> + <citations> + <citation type="doi">doi.org/10.5281/zenodo.5907920</citation> + </citations> + </xml> + <xml name="SRS_vegetpheno"> + <citations> + <citation type="doi">doi.org/10.1016/j.rse.2013.09.022</citation> + </citations> + </xml> + <xml name="topic"> + <edam_topics> + <edam_topic>topic_0610</edam_topic> + <edam_topic>topic_3050</edam_topic> + </edam_topics> + </xml> +</macros> |
| b |
| diff -r 000000000000 -r cf69ad260611 pca_raster.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pca_raster.r Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| @@ -0,0 +1,75 @@ +#Rscript + +########################################### +## Getting PCA raster ## +########################################### + +#####Packages : stars +# utils +# biodivmapr +# raster +# sf +# mapview +# leafpop +# RColorBrewer +# labdsv +# rgdal +# ggplot2 +# gridExtra +## remotes::install_github("jbferet/biodivMapR") +#####Load arguments + +args <- commandArgs(trailingOnly = TRUE) + +#####Import the S2 data + +if (length(args) < 1) { + stop("This tool needs at least 1 argument") +}else { + data_raster <- args[1] + rasterheader <- args[2] + data <- args[3] + typepca <- as.character(args[4]) + source(args[5]) +} + +################################################################################ +## DEFINE PARAMETERS FOR DATASET TO BE PROCESSED ## +################################################################################ +# expected to be in ENVI HDR + +if (data_raster == "") { + #Create a directory where to unzip your folder of data + dir.create("data_dir") + unzip(data, exdir = "data_dir") + # Path to raster + data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl") + input_image_file <- file.path("data_dir/results/Reflectance", data_raster[1]) + input_header_file <- file.path("data_dir/results/Reflectance", data_raster[2]) + +} else { + input_image_file <- file.path(getwd(), data_raster, fsep = "/") + input_header_file <- file.path(getwd(), rasterheader, fsep = "/") +} + +################################################################################ +## PROCESS IMAGE ## +################################################################################ +# 1- Filter data in order to discard non vegetated / shaded / cloudy pixels +print("PERFORM PCA ON RASTER") +pca_output <- biodivMapR::perform_PCA(Input_Image_File = input_image_file, Input_Mask_File = input_mask_file, + Output_Dir = output_dir, TypePCA = typepca, FilterPCA = filterpca, nbCPU = nbcpu, MaxRAM = maxram) + + +pca_path <- file.path(output_dir, basename(data_raster), typepca, "PCA", "OutputPCA_8_PCs") +pca_raster <- raster::raster(pca_path) +get_pca <- convert_raster(pca_raster) + +colnames(get_pca) <- c("PCA", "longitude", "latitude") +plot_indices(get_pca, titre = "PCA") + +write.table(get_pca, file = "PCA.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE) +#### Get the raster layer files +pca_files <- file.path("RESULTS", basename(data_raster), typepca, "PCA") +to_dir_short <- output_dir +file.copy(pca_files, to_dir_short) #copy files from long to short paths |
| b |
| diff -r 000000000000 -r cf69ad260611 preprocess_S2.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/preprocess_S2.r Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| @@ -0,0 +1,121 @@ +#Rscript + +########################################### +## Preprocessing Sentinel 2 data ## +########################################### + +#####Packages : sen2r, +# jqr, +# protolite, +# raster, +# sf, +# rgeos, +# sp, +# raster, +# stars, +# stringr, +# progress, +# rgdal, +# R.utils, +# gdalUtils, +# fasterize, +# XML, +# XML2 + +#####Load arguments + +args <- commandArgs(trailingOnly = TRUE) + +if (length(args) < 1) { + stop("This tool needs at least 1 argument") +}else { + data <- args[1] + source(args[2]) + data_source <- as.character(args[3]) + sat_type <- as.character(args[4]) +} + +##____________________________________________________________________## +## Define where data is stored and where to write results ## +##--------------------------------------------------------------------## + +#Create a directory where to unzip your folder of data +dir.create("data_dir") +unzip(data, exdir = "data_dir") + +# Result directory +result_path <- "results" +dir.create(path = result_path, showWarnings = FALSE, recursive = TRUE) + +#Csv file for output useless but needed for linter +write.csv(data_source, "Mission.csv") + +# define raster path +if (data_source == "SAFE") { + path_s2 <- file.path("data_dir", list.files("data_dir", pattern = ".SAFE")) + #To define the level and know if a correction is needed (convert not ready yet) + level_info <- get_s2_level(path_s2) + if (level_info == "L1C") { + stop("! This tool works for data of L2A level and NOT for the L1C level which is currently a work in progress !") + } +}else { + path_s2 <- file.path("data_dir") +} + +##____________________________________________________________________## +## Extract, resample & stack data ## +##--------------------------------------------------------------------## +# define resolution +resolution <- 10 +# define source of data +s2source <- data_source + +s2obj <- extract_from_s2_l2a(path_dir_s2 = path_s2, + path_vector = NULL, + s2source = s2source, + resolution = resolution) + +##____________________________________________________________________## +## Write CLOUD MASK ## +##--------------------------------------------------------------------## + +# directory for cloud mask +cloud_path <- file.path(result_path, "CloudMask") +dir.create(path = cloud_path, showWarnings = FALSE, recursive = TRUE) +# Filename for cloud mask +cloudmasks <- save_cloud_s2(s2_stars = s2obj$s2_stack, + cloud_path = cloud_path, + s2source = s2source, saveraw = TRUE) + +zip_cloud <- file.path("Cloud.zip") +zip::zip(zip_cloud, cloud_path) +##____________________________________________________________________## +## Write REFLECTANCE ## +##--------------------------------------------------------------------## + +# directory for Reflectance +refl_dir <- file.path(result_path, "Reflectance") +dir.create(path = refl_dir, showWarnings = FALSE, recursive = TRUE) + +if (data_source == "SAFE") { + # filename for Reflectance + refl_path <- file.path(refl_dir, paste(basename(s2obj$s2_bands$GRANULE), "_Refl", sep = "")) + + # Save Reflectance file as ENVI image with BIL interleaves + tile_s2 <- substring(strsplit(basename(s2obj$s2_bands$GRANULE), "_")[[1]][2], 2) + dateacq_s2 <- as.Date(substring(strsplit(basename(s2obj$s2_bands$GRANULE), "_")[[1]][4], 1, 8), format = "%Y%m%d") +}else { + # filename for Reflectance + refl_path <- file.path(refl_dir, paste(basename(s2obj$s2_bands$path_tile_s2), "_Refl", sep = "")) + + # Save Reflectance file as ENVI image with BIL interleaves + tile_s2 <- substring(strsplit(basename(s2obj$s2_bands$path_tile_s2), "_")[[1]][2], 2) + dateacq_s2 <- as.Date(substring(strsplit(basename(s2obj$s2_bands$path_tile_s2), "_")[[1]][4], 1, 8), format = "%Y%m%d") +} + +save_data <- save_reflectance_s2(s2_stars = s2obj$s2_stack, refl_path = refl_path, + s2sat = sat_type, tile_s2 = tile_s2, dateacq_s2 = dateacq_s2, + format = "ENVI", datatype = "Int16", mtd = s2obj$s2_bands$metadata, mtd_msi = s2obj$s2_bands$metadata_MSI) + +zip_files <- file.path("Refl.zip") +zip::zip(zip_files, refl_dir) |
| b |
| diff -r 000000000000 -r cf69ad260611 prosail-master/R/Lib_PROSAIL.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prosail-master/R/Lib_PROSAIL.R Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| b'@@ -0,0 +1,1233 @@\n+# ============================================================================= =\n+# prosail\n+# Lib_PROSAIL.R\n+# ============================================================================= =\n+# PROGRAMMERS:\n+# Jean-Baptiste FERET <jb.feret@teledetection.fr >\n+# Florian de BOISSIEU <fdeboiss@gmail.com >\n+# copyright 2019 / 11 Jean-Baptiste FERET\n+# ============================================================================= =\n+# This Library includes functions dedicated to PROSAIL simulation\n+# SAIL versions available are 4SAIL and 4SAIL2\n+# ============================================================================= =\n+\n+#" computes bidirectional reflectance factor based on outputs from PROSAIL and sun position\n+#"\n+#" The direct and diffuse light are taken into account as proposed by:\n+#" Francois et al. (2002) conversion of 400-1100 nm vegetation albedo\n+#" measurements into total shortwave broadband albedo using a canopy\n+#" radiative transfer model, Agronomie\n+#" Es = direct\n+#" Ed = diffuse\n+#"\n+#" @param rdot numeric. Hemispherical-directional reflectance factor in viewing direction\n+#" @param rsot numeric. Bi-directional reflectance factor\n+#" @param tts numeric. Solar zenith angle\n+#" @param specatm_sensor list. direct and diffuse radiation for clear conditions\n+#" @return brf numeric. Bidirectional reflectance factor\n+#" @export\n+compute_brf <- function(rdot, rsot, tts, specatm_sensor) {\n+\n+ ############################## #\n+ ##\tdirect / diffuse light\t##\n+ ############################## #\n+ es <- specatm_sensor$Direct_Light\n+ ed <- specatm_sensor$Diffuse_Light\n+ rd <- pi / 180\n+ skyl <- 0.847 - 1.61 * sin((90 - tts) * rd) + 1.04 * sin((90 - tts) * rd) * sin((90 - tts) * rd) # diffuse radiation (Francois et al., 2002)\n+ pardiro <- (1 - skyl) * es\n+ pardifo <- skyl * ed\n+ brf <- (rdot * pardifo + rsot * pardiro) / (pardiro + pardifo)\n+ return(brf)\n+}\n+\n+#" Performs PROSAIL simulation based on a set of combinations of input parameters\n+#" @param spec_sensor list. Includes optical constants required for PROSPECT\n+#" refractive index, specific absorption coefficients and corresponding spectral bands\n+#" @param input_prospect list. PROSPECT input variables\n+#" @param n numeric. Leaf structure parameter\n+#" @param chl numeric. chlorophyll content (microg.cm-2)\n+#" @param car numeric. carotenoid content (microg.cm-2)\n+#" @param ant numeric. anthocyain content (microg.cm-2)\n+#" @param brown numeric. brown pigment content (Arbitrary units)\n+#" @param ewt numeric. Equivalent Water Thickness (g.cm-2)\n+#" @param lma numeric. Leaf Mass per Area (g.cm-2)\n+#" @param prot numeric. protein content (g.cm-2)\n+#" @param cbc numeric. nonprotcarbon-based constituent content (g.cm-2)\n+#" @param alpha numeric. Solid angle for incident light at surface of leaf (simulation of roughness)\n+#" @param typelidf numeric. Type of leaf inclination distribution function\n+#" @param lidfa numeric.\n+#" if typelidf == 1, controls the average leaf slope\n+#" if typelidf == 2, corresponds to average leaf angle\n+#" @param lidfb numeric.\n+#" if typelidf == 1, unused\n+#" if typelidf == 2, controls the distribution"s bimodality\n+#" @param lai numeric. Leaf Area Index\n+#" @param q numeric. Hot Spot parameter\n+#" @param tts numeric. Sun zeith angle\n+#" @param tto numeric. Observer zeith angle\n+#" @param psi numeric. Azimuth Sun / Observer\n+#" @param rsoil numeric. Soil reflectance\n+#" @param fraction_brown numeric. Fraction of brown leaf area\n+#" @param diss numeric. Layer dissociation factor\n+#" @param cv numeric. vertical crown cover percentage\n+#" = % ground area covered with crowns as seen from nadir direction\n+#" @param zeta numeric. Tree shape factor\n+#" = ratio of crown diameter to crown height\n+#" @param sailversion character. choose between 4SAIL and 4SAIL2\n+#" @param brownvegetation list. Defines optical properties for brown vegetation, if not nULL\n+#" - WVL, reflectance, Transmittance\n+#" - Set to '..b'r) flux\n+#" @param l numeric.\n+#" @param t numeric. Leaf Area Index\n+#" @return jout numeric.\n+#" @export\n+jfunc3 <- function(k, l, t) {\n+ out <- (1. - exp(-(k + l) * t)) / (k + l)\n+ return(out)\n+}\n+\n+\n+#" j4 function for treating (near) conservative scattering\n+#"\n+#" @param m numeric. Extinction coefficient for direct (solar or observer) flux\n+#" @param t numeric. Leaf Area Index\n+#" @return jout numeric.\n+#" @export\n+jfunc4 <- function(m, t) {\n+\n+ del <- m * t\n+ out <- 0 * del\n+ out[del > 1e-3] <- (1 - exp(-del)) / (m * (1 + exp(-del)))\n+ out[del <= 1e-3] <- 0.5 * t * (1. - del * del / 12.)\n+ return(out)\n+}\n+\n+\n+#" compute volume scattering functions and interception coefficients\n+#" for given solar zenith, viewing zenith, azimuth and leaf inclination angle.\n+#"\n+#" @param tts numeric. solar zenith\n+#" @param tto numeric. viewing zenith\n+#" @param psi numeric. azimuth\n+#" @param ttl numeric. leaf inclination angle\n+#" @return res list. includes chi_s, chi_o, frho, ftau\n+#" @export\n+volscatt <- function(tts, tto, psi, ttl) {\n+ #********************************************************************************\n+ #*\tchi_s\t= interception functions\n+ #*\tchi_o\t= interception functions\n+ #*\tfrho\t= function to be multiplied by leaf reflectance rho\n+ #*\tftau\t= functions to be multiplied by leaf transmittance tau\n+ #********************************************************************************\n+ #\tWout Verhoef, april 2001, for CROMA\n+\n+ rd <- pi / 180\n+ costs <- cos(rd * tts)\n+ costo <- cos(rd * tto)\n+ sints <- sin(rd * tts)\n+ sinto <- sin(rd * tto)\n+ cospsi <- cos(rd * psi)\n+ psir <- rd * psi\n+ costl <- cos(rd * ttl)\n+ sintl <- sin(rd * ttl)\n+ cs <- costl * costs\n+ co <- costl * costo\n+ ss <- sintl * sints\n+ so <- sintl * sinto\n+\n+ #c ..............................................................................\n+ #c betas -bts- and betao -bto- computation\n+ #c Transition angles (beta) for solar (betas) and view (betao) directions\n+ #c if thetav + thetal > pi / 2, bottom side of the leaves is observed for leaf azimut\n+ #c interval betao + phi<leaf azimut<2pi-betao + phi.\n+ #c if thetav + thetal<pi / 2, top side of the leaves is always observed, betao=pi\n+ #c same consideration for solar direction to compute betas\n+ #c ..............................................................................\n+\n+ cosbts <- 5\n+ if (abs(ss) > 1e-6) {\n+ cosbts <- -cs / ss\n+ }\n+ cosbto <- 5\n+ if (abs(so) > 1e-6) {\n+ cosbto <- -co / so\n+ }\n+\n+ if (abs(cosbts) < 1) {\n+ bts <- acos(cosbts)\n+ ds <- ss\n+ } else {\n+ bts <- pi\n+ ds <- cs\n+ }\n+ chi_s <- 2. / pi * ((bts - pi * .5) * cs + sin(bts) * ss)\n+ if (abs(cosbto) < 1) {\n+ bto <- acos(cosbto)\n+ doo <- so\n+ } else if (tto < 90) {\n+ bto <- pi\n+ doo <- co\n+ } else {\n+ bto <- 0\n+ doo <- -co\n+ }\n+ chi_o <- 2. / pi * ((bto - pi * .5) * co + sin(bto) * so)\n+\n+ #c ..............................................................................\n+ #c computation of auxiliary azimut angles bt1, bt2, bt3 used\n+ #c for the computation of the bidirectional scattering coefficient w\n+ #c .............................................................................\n+\n+ btran1 <- abs(bts - bto)\n+ btran2 <- pi - abs(bts + bto - pi)\n+\n+ if (psir <= btran1) {\n+ bt1 <- psir\n+ bt2 <- btran1\n+ bt3 <- btran2\n+ } else {\n+ bt1 <- btran1\n+ if (psir <= btran2) {\n+ bt2 <- psir\n+ bt3 <- btran2\n+ } else {\n+ bt2 <- btran2\n+ bt3 <- psir\n+ }\n+ }\n+ t1 <- 2. * cs * co + ss * so * cospsi\n+ t2 <- 0\n+ if (bt2 > 0) {\n+ t2 <- sin(bt2) * (2. * ds * doo + ss * so * cos(bt1) * cos(bt3))\n+ }\n+\n+ denom <- 2. * pi * pi\n+ frho <- ((pi - bt2) * t1 + t2) / denom\n+ ftau <- (-bt2 * t1 + t2) / denom\n+\n+ if (frho < 0) {\n+ frho <- 0\n+ }\n+ if (ftau < 0) {\n+ ftau <- 0\n+ }\n+ res <- list("chi_s" = chi_s, "chi_o" = chi_o, "frho" = frho, "ftau" = ftau)\n+ return(res)\n+}\n' |
| b |
| diff -r 000000000000 -r cf69ad260611 prosail-master/R/Lib_PROSAIL_HybridInversion.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prosail-master/R/Lib_PROSAIL_HybridInversion.R Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| b'@@ -0,0 +1,554 @@\n+# ============================================================================= =\n+# prosail\n+# Lib_PROSAIL_HybridInversion.R\n+# ============================================================================= =\n+# PROGRAMMERS:\n+# Jean-Baptiste FERET <jb.feret@teledetection.fr>\n+# Florian de BOISSIEU <fdeboiss@gmail.com>\n+# Copyright 2019 / 11 Jean-Baptiste FERET\n+# ============================================================================= =\n+# This Library includes functions dedicated to PROSAIL inversion using hybrid\n+# approach based on SVM regression\n+# ============================================================================= =\n+\n+\n+#" This function applies SVR model on raster data in order to estimate\n+#" vegetation biophysical properties\n+#"\n+#" @param raster_path character. path for a raster file\n+#" @param hybridmodel list. hybrid models produced from train_prosail_inversion\n+#" each element of the list corresponds to a set of hybrid models for a vegetation parameter\n+#" @param pathout character. path for directory where results are written\n+#" @param selectedbands list. list of spectral bands to be selected from raster (identified by name of vegetation parameter)\n+#" @param bandname character. spectral bands corresponding to the raster\n+#" @param maskraster character. path for binary mask defining ON (1) and OFF (0) pixels in the raster\n+#" @param multiplyingfactor numeric. multiplying factor used to write reflectance in the raster\n+#" --> PROSAIL simulates reflectance between 0 and 1, and raster data expected in the same range\n+#"\n+#" @return None\n+#" @importFrom progress progress_bar\n+#" @importFrom stars read_stars\n+#" @importFrom raster raster brick blockSize readStart readStop getValues writeStart writeStop writeValues\n+#" @import rgdal\n+#" @export\n+apply_prosail_inversion <- function(raster_path, hybridmodel, pathout,\n+ selectedbands, bandname,\n+ maskraster = FALSE, multiplyingfactor = 10000) {\n+\n+ # explain which biophysical variables will be computed\n+ bpvar <- names(hybridmodel)\n+ print("The following biophysical variables will be computed")\n+ print(bpvar)\n+\n+ # get image dimensions\n+ if (attr(rgdal::GDALinfo(raster_path, returnStats = FALSE), "driver") == "ENVI") {\n+ hdr <- read_envi_header(get_hdr_name(raster_path))\n+ dimsraster <- list("rows" = hdr$lines, "cols" = hdr$samples, "bands" = hdr$bands)\n+ } else {\n+ dimsraster <- dim(read_stars(raster_path))\n+ dimsraster <- list("rows" = as.numeric(dimsraster[2]), "cols" = as.numeric(dimsraster[1]), "bands" = as.numeric(dimsraster[3]))\n+ }\n+\n+ # Produce a map for each biophysical property\n+ for (parm in bpvar) {\n+ print(paste("Computing", parm, sep = " "))\n+ # read by chunk to avoid memory problem\n+ blk <- blockSize(brick(raster_path))\n+ # reflectance file\n+ r_in <- readStart(brick(raster_path))\n+ # mask file\n+ r_inmask <- FALSE\n+ if (maskraster == FALSE) {\n+ selectpixels <- "ALL"\n+ } else if (!maskraster == FALSE) {\n+ if (file.exists(maskraster)) {\n+ r_inmask <- readStart(raster(maskraster))\n+ } else if (!file.exists(maskraster)) {\n+ message("WARNING: Mask file does not exist:")\n+ print(maskraster)\n+ message("Processing all image")\n+ selectpixels <- "ALL"\n+ }\n+ }\n+ # initiate progress bar\n+ pgbarlength <- length(hybridmodel[[parm]]) * blk$n\n+ pb <- progress_bar$new(\n+ format = "Hybrid inversion on raster [:bar] :percent in :elapsedfull, estimated time remaining :eta",\n+ total = pgbarlength, clear = FALSE, width = 100)\n+\n+ # output files\n+ bpvarpath <- file.path(pathout, paste(basename(raster_path), parm, sep = "_"))\n+ bpvarsdpath <- file.path(pathout, paste(basename(raster_path), parm, "STD", sep = "_"))\n+ r_outmean <- writeStart(raster(raster_path), filename = bpvarpath, format = "ENVI", overwrite = TRUE)\n+ r_outsd <- wr'..b' force4lowlai = force4lowlai)\n+ }\n+ if (sailversion == "4SAIL2") {\n+ # Definition of Cv && update LAI\n+ maxlai <- min(c(maxval$lai), 4)\n+ inputprosail$Cv <- NA * inputprosail$lai\n+ inputprosail$Cv[which(inputprosail$lai > maxlai)] <- 1\n+ inputprosail$Cv[which(inputprosail$lai <= maxlai)] <- (1 / maxlai) + inputprosail$lai[which(inputprosail$lai <= maxlai)] / (maxlai + 1)\n+ inputprosail$Cv <- inputprosail$Cv * matrix(rnorm(length(inputprosail$Cv), mean = 1, sd = 0.1))\n+ inputprosail$Cv[which(inputprosail$Cv < 0)] <- 0\n+ inputprosail$Cv[which(inputprosail$Cv > 1)] <- 1\n+ inputprosail$Cv[which(inputprosail$lai > maxlai)] <- 1\n+ inputprosail$fraction_brown <- 0 + 0 * inputprosail$lai\n+ inputprosail$diss <- 0 + 0 * inputprosail$lai\n+ inputprosail$Zeta <- 0.2 + 0 * inputprosail$lai\n+ inputprosail$lai <- inputprosail$lai * inputprosail$Cv\n+ }\n+\n+ # generate LUT of BRF corresponding to inputprosail, for a sensor\n+ brf_lut <- Generate_LUT_BRF(sailversion = sailversion, inputprosail = inputprosail,\n+ specprospect = specprospect, specsoil = specsoil, specatm = specatm)\n+\n+ # write parameters LUT\n+ output <- matrix(unlist(inputprosail), ncol = length(inputprosail), byrow = FALSE)\n+ filename <- file.path(path_results, "PROSAIL_LUT_InputParms.txt")\n+ write.table(x = format(output, digits = 3), file = filename, append = FALSE, quote = FALSE,\n+ col.names = names(inputprosail), row.names = FALSE, sep = "\\t")\n+ # Write BRF LUT corresponding to parameters LUT\n+ filename <- file.path(path_results, "PROSAIL_LUT_reflectance.txt")\n+ write.table(x = format(t(brf_lut), digits = 5), file = filename, append = FALSE, quote = FALSE,\n+ col.names = specprospect$lambda, row.names = FALSE, sep = "\\t")\n+\n+ # Which bands will be used for inversion?\n+ if (is.null(bands2select)) {\n+ bands2select <- list()\n+ for (parm in parms2estimate) {\n+ bands2select[[parm]] <- seq(1, length(specprospect$lambda))\n+ }\n+ }\n+ # Add gaussian noise to reflectance LUT: one specific LUT per parameter\n+ if (is.null(noiselevel)) {\n+ noiselevel <- list()\n+ for (parm in parms2estimate) {\n+ noiselevel[[parm]] <- 0.01\n+ }\n+ }\n+\n+ # produce LIT with noise\n+ brf_lut_noise <- list()\n+ for (parm in parms2estimate) {\n+ brf_lut_noise[[parm]] <- brf_lut[bands2select[[parm]], ] + brf_lut[bands2select[[parm]], ] * matrix(rnorm(nrow(brf_lut[bands2select[[parm]], ]) * ncol(brf_lut[bands2select[[parm]], ]),\n+ 0, noiselevel[[parm]]), nrow = nrow(brf_lut[bands2select[[parm]], ]))\n+ }\n+\n+ ###===================================================================###\n+ ### PERFORM HYBRID INVERSION ###\n+ ###===================================================================###\n+ # train SVR for each variable and each run\n+ modelsvr <- list()\n+ for (parm in parms2estimate) {\n+ colparm <- which(parm == names(inputprosail))\n+ inputvar <- inputprosail[[colparm]]\n+ modelsvr[[parm]] <- prosail_hybrid_train(brf_lut = brf_lut_noise[[parm]], inputvar = inputvar,\n+ figplot = figplot, nbensemble = nbmodels, withreplacement = replacement)\n+ }\n+ return(modelsvr)\n+}\n+\n+#" writes ENVI hdr file\n+#"\n+#" @param hdr content to be written\n+#" @param hdrpath Path of the hdr file\n+#"\n+#" @return None\n+#" @importFrom stringr str_count\n+#" @export\n+write_envi_header <- function(hdr, hdrpath) {\n+ h <- lapply(hdr, function(x) {\n+ if (length(x) > 1 || (is.character(x) && stringr::str_count(x, "\\\\w + ") > 1)) {\n+ x <- paste0("{", paste(x, collapse = ","), "}")\n+ }\n+ # convert last numerics\n+ x <- as.character(x)\n+ })\n+ writeLines(c("ENVI", paste(names(hdr), h, sep = "=")), con = hdrpath)\n+ return(invisible())\n+}\n' |
| b |
| diff -r 000000000000 -r cf69ad260611 prosail-master/R/Lib_SpectralIndices.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prosail-master/R/Lib_SpectralIndices.R Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| b'@@ -0,0 +1,679 @@\n+# ============================================================================== =\n+# prosail\n+# Lib_spectralindices.R\n+# ============================================================================== =\n+# PROGRAMMERS:\n+# Jean-Baptiste FERET <jb.feret@teledetection.fr>\n+# Florian de BOISSIEU <fdeboiss@gmail.com>\n+# Copyright 2019/11 Jean-Baptiste FERET\n+# ============================================================================== =\n+# This Library includes aims at computing spectral indices from reflectance data\n+# ============================================================================== =\n+\n+#" This function computes Area under curve for continuum removed reflectances\n+#" See Malenovsk\xc3\xbd et al. (2013) for details\n+#" http://dx.doi.org/10.1016/j.rse.2012.12.015\n+#"\n+#" @param refl RasterBrick, RasterStack or list. Raster bands in the order of sensorbands.\n+#" @param sensorbands numeric. vector containing central wavelength for each spectral band in the image\n+#" @param aucminmax list. wavelengths of lower and upper boundaries ("CRmin" and "CRmax")\n+#" @param reflfactor numeric. multiplying factor used to write reflectance in image (==10000 for S2)\n+#"\n+#" @return aucval raster\n+#" @export\n+auc <- function(refl, sensorbands, aucminmax, reflfactor = 1) {\n+\n+ aucbands <- list()\n+ aucbands[["CRmin"]] <- sensorbands[get_closest_bands(sensorbands, aucminmax[["CRmin"]])]\n+ aucbands[["CRmax"]] <- sensorbands[get_closest_bands(sensorbands, aucminmax[["CRmax"]])]\n+ bands <- get_closest_bands(sensorbands, aucbands)\n+ for (i in bands[["CRmin"]]:bands[["CRmax"]]) {\n+ if (is.na(match(i, bands))) {\n+ aucbands[[paste("B", i, sep = "")]] <- sensorbands[i]\n+ }\n+ }\n+ # compute continuum removal for all spectral bands\n+ cr <- cr_wl(refl = refl, sensorbands = sensorbands,\n+ crbands = aucbands, reflfactor = reflfactor)\n+\n+ wl <- sort(unlist(aucbands), decreasing = FALSE)\n+ aucval <- 0.5 * (1 - cr[[1]]) * (wl[2] - wl[1])\n+ for (i in 2:length(cr)) {\n+ aucval <- aucval + 0.5 * (2 - cr[[i - 1]] - cr[[i]]) * (wl[i + 1] - wl[i])\n+ }\n+ aucval <- aucval + 0.5 * (1 - cr[[length(cr)]]) * (wl[i + 2] - wl[i + 1])\n+ return(aucval)\n+}\n+\n+#" This function extracts boundaries to be used to compute continuum from reflectance data\n+#"\n+#" @param refl RasterBrick, RasterStack or list. Raster bands in the order of sensorbands.\n+#" @param sensorbands numeric. vector containing central wavelength for each spectral band in the image\n+#" @param crbands list. list of spectral bands (central wavelength) including CRmin and CRmax\n+#" @param reflfactor numeric. multiplying factor used to write reflectance in image ( == 10000 for S2)\n+#"\n+#" @return crminmax list. list of rasters corresponding to minimum and maximum wavelengths\n+#" @export\n+crbound <- function(refl, sensorbands, crbands, reflfactor = 1) {\n+\n+ # get closest spectral bands from CR1 and CR2\n+ bands <- get_closest_bands(sensorbands, list(crbands[["CRmin"]], crbands[["CRmax"]]))\n+ wl <- sensorbands[bands]\n+ # get equation for line going from CR1 to CR2\n+ crminmax <- readrasterbands(refl = refl, bands = bands, reflfactor = reflfactor)\n+ names(crminmax) <- paste("wl_", as.character(wl), sep = "")\n+ return(crminmax)\n+}\n+\n+#" This function extracts boundaries to be used to compute continuum from reflectance data\n+#"\n+#" @param refl RasterBrick, RasterStack or list. Raster bands in the order of sensorbands.\n+#" @param sensorbands numeric. vector containing central wavelength for each spectral band in the image\n+#" @param crbands list. list of spectral bands (central wavelength) including CRmin and CRmax\n+#" @param reflfactor numeric. multiplying factor used to write reflectance in image ( == 10000 for S2)\n+#"\n+#" @return outlier_iqr numeric. band numbers of original sensor corresponding to S2\n+#" @importFrom progress progress_bar\n+#" @export\n+cr_wl <- function(refl, sensorbands, crbands, reflfactor = 1) {\n+\n+ # Make sure CRmin and CRma'..b'[["B8"]]] - refl[, sen2s2[["B6"]]]) / (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B6"]]])\n+ spectralindices$RE_NDVI <- re_ndvi\n+ }\n+ if ("RE_NDWI" %in% sel_indices) {\n+ re_ndwi <- (refl[, sen2s2[["B4"]]] - refl[, sen2s2[["B6"]]]) / (refl[, sen2s2[["B4"]]] + refl[, sen2s2[["B6"]]])\n+ spectralindices$RE_NDWI <- re_ndwi\n+ }\n+ if ("S2REP" %in% sel_indices) {\n+ s2rep <- 705 + 35 * (0.5 * (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B5"]]]) - refl[, sen2s2[["B6"]]]) / (refl[, sen2s2[["B7"]]] - refl[, sen2s2[["B6"]]])\n+ spectralindices$S2REP <- s2rep\n+ }\n+ if ("SAVI" %in% sel_indices) {\n+ savi <- 1.5 * (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B5"]]]) / (refl[, sen2s2[["B8"]]] + refl[, sen2s2[["B5"]]] + 0.5)\n+ spectralindices$SAVI <- savi\n+ }\n+ if ("SIPI" %in% sel_indices) {\n+ sipi <- (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B2"]]]) / (refl[, sen2s2[["B8"]]] - refl[, sen2s2[["B4"]]])\n+ spectralindices$SIPI <- sipi\n+ }\n+ if ("SR" %in% sel_indices) {\n+ sr <- refl[, sen2s2[["B8"]]] / refl[, sen2s2[["B4"]]]\n+ spectralindices$SR <- sr\n+ }\n+ if ("CR_SWIR" %in% sel_indices) {\n+ cr_swir <- refl[, sen2s2[["B11"]]] / (refl[, sen2s2[["B8A"]]] + (s2bands$B11 - s2bands$B8A) * (refl[, sen2s2[["B12"]]] - refl[, sen2s2[["B8A"]]]) / (s2bands$B12 - s2bands$B8A))\n+ spectralindices$CR_SWIR <- cr_swir\n+ }\n+ res <- list("spectralindices" = spectralindices, "listindices" = listindices)\n+ return(res)\n+}\n+\n+#" this function identifies the bands of a given sensor with closest match to its spectral characteristics\n+#"\n+#" @param sensorbands numeric. wavelength in nanometer of the sensor of interest\n+#" @param listbands numeric or list. Named vector or list of spectral bands corresponding to sensor\n+#"\n+#" @return numeric. band numbers of original sensor\n+#" @export\n+get_closest_bands <- function(sensorbands, listbands) {\n+ sapply(listbands, function(x) {\n+ b <- which.min(abs(sensorbands - x))\n+ names(b) <- ""\n+ b\n+ })\n+}\n+\n+#" This function computes interquartile range (IQR) criterion, which can be used\n+#" as a criterion for outlier detection\n+#"\n+#" @param distval numeric. vector of distribution of values\n+#" @param weightirq numeric. weighting factor appplied to IRQ to define lower and upper boudaries for outliers\n+#"\n+#" @return outlier_iqr numeric. band numbers of original sensor corresponding to S2\n+#" @importFrom stats IQR quantile\n+#" @export\n+iqr_outliers <- function(distval, weightirq = 1.5) {\n+ iqr <- IQR(distval, na.rm = TRUE)\n+ range_iqr <- c(quantile(distval, 0.25, na.rm = TRUE), quantile(distval, 0.75, na.rm = TRUE))\n+ outlier_iqr <- c(range_iqr[1] - weightirq * iqr, range_iqr[2] + weightirq * iqr)\n+ return(outlier_iqr)\n+}\n+\n+#" This function selects bands from a raster or stars object\n+#"\n+#" @param refl RasterBrick, RasterStack or list. Raster bands in the order of sensorbands.\n+#" @param bands numeric. rank of bands to be read in refl\n+#" @param reflfactor numeric. multiplying factor used to write reflectance in image ( == 10000 for S2)\n+#"\n+#" @return robj list. R object (default is raster, stars if refl is stars object)\n+#" @importFrom raster subset stack\n+#" @export\n+readrasterbands <- function(refl, bands, reflfactor = 1) {\n+\n+ # get equation for line going from CR1 to CR2\n+ classraster <- class(refl)\n+ if (classraster == "RasterBrick" || classraster == "RasterStack" || classraster == "stars") {\n+ # if !reflfactor == 1 then apply a reflectance factor\n+ if (classraster == "stars") {\n+ robj <- refl[bands]\n+ } else {\n+ robj <- raster::subset(refl, bands)\n+ }\n+ if (!reflfactor == 1) {\n+ robj <- robj / reflfactor\n+ }\n+ } else if (is.list(refl)) {\n+ robj <- raster::stack(refl[bands]) # checks that all rasters have same crs/extent\n+ if (!reflfactor == 1) {\n+ robj <- robj / reflfactor\n+ }\n+ } else {\n+ stop("refl is expected to be a RasterStack, RasterBrick, Stars object or a list of rasters")\n+ }\n+ return(robj)\n+}\n' |
| b |
| diff -r 000000000000 -r cf69ad260611 test-data/12a0b625-9ad5-4251-a57a-305e22edef2e.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/12a0b625-9ad5-4251-a57a-305e22edef2e.xml Mon Jan 09 13:36:02 2023 +0000 |
| b |
| b'@@ -0,0 +1,661 @@\n+<?xml version="1.0" encoding="UTF-8"?>\n+<gmd:MD_Metadata xmlns:gmd="http://www.isotc211.org/2005/gmd"\n+ xmlns:gml="http://www.opengis.net/gml/3.2"\n+ xmlns:gts="http://www.isotc211.org/2005/gts"\n+ xmlns:gco="http://www.isotc211.org/2005/gco"\n+ xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"\n+ xsi:schemaLocation="http://www.isotc211.org/2005/gmd http://www.isotc211.org/2005/gmd/gmd.xsd http://www.isotc211.org/2005/gmx http://www.isotc211.org/2005/gmx/gmx.xsd http://www.isotc211.org/2005/srv http://schemas.opengis.net/iso/19139/20060504/srv/srv.xsd">\n+ <gmd:fileIdentifier>\n+ <gco:CharacterString>12a0b625-9ad5-4251-a57a-305e22edef2e</gco:CharacterString>\n+ </gmd:fileIdentifier>\n+ <gmd:language>\n+ <gmd:LanguageCode codeList="http://www.loc.gov/standards/iso639-2/" codeListValue="fre"/>\n+ </gmd:language>\n+ <gmd:characterSet>\n+ <gmd:MD_CharacterSetCode codeListValue="utf8"\n+ codeList="http://standards.iso.org/iso/19139/resources/gmxCodelists.xml#MD_CharacterSetCode"/>\n+ </gmd:characterSet>\n+ <gmd:hierarchyLevel xmlns:gn="http://www.fao.org/geonetwork"\n+ xmlns:gmx="http://www.isotc211.org/2005/gmx"\n+ xmlns:srv="http://www.isotc211.org/2005/srv"\n+ xmlns:xlink="http://www.w3.org/1999/xlink">\n+ <gmd:MD_ScopeCode codeList="http://standards.iso.org/iso/19139/resources/gmxCodelists.xml#MD_ScopeCode"\n+ codeListValue="dataset"/>\n+ </gmd:hierarchyLevel>\n+ <gmd:contact>\n+ <gmd:CI_ResponsibleParty>\n+ <gmd:individualName>\n+ <gco:CharacterString>Santiago Poggio</gco:CharacterString>\n+ </gmd:individualName>\n+ <gmd:organisationName>\n+ <gco:CharacterString>IFEVA/ Cat\xc3\xa9dra de producci\xc3\xb3n vegetal, facultad de agronomia, Universidad de Buenos Aires/ CONICET</gco:CharacterString>\n+ </gmd:organisationName>\n+ <gmd:positionName gco:nilReason="missing">\n+ <gco:CharacterString/>\n+ </gmd:positionName>\n+ <gmd:contactInfo>\n+ <gmd:CI_Contact>\n+ <gmd:phone>\n+ <gmd:CI_Telephone>\n+ <gmd:voice gco:nilReason="missing">\n+ <gco:CharacterString/>\n+ </gmd:voice>\n+ <gmd:facsimile gco:nilReason="missing">\n+ <gco:CharacterString/>\n+ </gmd:facsimile>\n+ </gmd:CI_Telephone>\n+ </gmd:phone>\n+ <gmd:address>\n+ <gmd:CI_Address>\n+ <gmd:deliveryPoint gco:nilReason="missing">\n+ <gco:CharacterString/>\n+ </gmd:deliveryPoint>\n+ <gmd:city gco:nilReason="missing">\n+ <gco:CharacterString/>\n+ </gmd:city>\n+ <gmd:administrativeArea gco:nilReason="missing">\n+ <gco:CharacterString/>\n+ </gmd:administrativeArea>\n+ <gmd:postalCode gco:nilReason="missing">\n+ <gco:CharacterString/>\n+ </gmd:postalCode>\n+ <gmd:country gco:nilReason="missing">\n+ <gco:CharacterString/>\n+ </gmd:country>\n+ <gmd:electronicMailAddress>\n+ <gco:CharacterString>spoggio@agro.uba.ar</gco:CharacterString>\n+ </gmd:electronicMailAddress>\n+ </gmd:CI_Address>\n+ </gmd:address>\n+ </gmd:CI_Contact>\n+ </gmd:contactInfo>\n+ <gmd:role>\n+ <gmd:CI_RoleCode codeList="http://standards.iso.org/iso/19139/resources/gmxCodelists.xml#CI_RoleCode"\n+ codeListValue="originator"/>\n+ </gmd:role>\n+ </gmd:CI_Resp'..b' <gmd:northBoundLatitude>\n+ <gco:Decimal>48.696899414063</gco:Decimal>\n+ </gmd:northBoundLatitude>\n+ </gmd:EX_GeographicBoundingBox>\n+ </gmd:geographicElement>\n+ </gmd:EX_Extent>\n+ </gmd:extent>\n+ <gmd:supplementalInformation>\n+ <gco:CharacterString>France, Pleine-Foug\xc3\xa8res</gco:CharacterString>\n+ </gmd:supplementalInformation>\n+ </gmd:MD_DataIdentification>\n+ </gmd:identificationInfo>\n+ <gmd:distributionInfo>\n+ <gmd:MD_Distribution>\n+ <gmd:distributionFormat>\n+ <gmd:MD_Format>\n+ <gmd:name>\n+ <gco:CharacterString>ESRI Shapefile</gco:CharacterString>\n+ </gmd:name>\n+ <gmd:version>\n+ <gco:CharacterString>1.0</gco:CharacterString>\n+ </gmd:version>\n+ </gmd:MD_Format>\n+ </gmd:distributionFormat>\n+ <gmd:transferOptions>\n+ <gmd:MD_DigitalTransferOptions/>\n+ </gmd:transferOptions>\n+ </gmd:MD_Distribution>\n+ </gmd:distributionInfo>\n+ <gmd:dataQualityInfo>\n+ <gmd:DQ_DataQuality>\n+ <gmd:scope>\n+ <gmd:DQ_Scope>\n+ <gmd:level>\n+ <gmd:MD_ScopeCode codeListValue="series"\n+ codeList="http://standards.iso.org/iso/19139/resources/gmxCodelists.xml#MD_ScopeCode"/>\n+ </gmd:level>\n+ </gmd:DQ_Scope>\n+ </gmd:scope>\n+ <gmd:report xmlns:gn="http://www.fao.org/geonetwork"\n+ xmlns:gmx="http://www.isotc211.org/2005/gmx"\n+ xmlns:srv="http://www.isotc211.org/2005/srv"\n+ xmlns:xlink="http://www.w3.org/1999/xlink">\n+ <gmd:DQ_DomainConsistency>\n+ <gmd:result>\n+ <gmd:DQ_ConformanceResult>\n+ <gmd:specification>\n+ <gmd:CI_Citation>\n+ <gmd:title>\n+ <gco:CharacterString>L\xe2\x80\x99article 7, paragraphe 1, de la directive 2007/2/CE correspond aux modalit\xc3\xa9s techniques de l\xe2\x80\x99interop\xc3\xa9rabilit\xc3\xa9 : il s\xe2\x80\x99agit du r\xc3\xa8glement relatif \xc3\xa0 l\xe2\x80\x99interop\xc3\xa9rabilit\xc3\xa9 : r\xc3\xa8glement n\xc2\xb01253/2013 du 21 octobre 2013 modifiant et compl\xc3\xa9tant le r\xc3\xa8glement n\xc2\xb01089/2010 du 23 novembre 2010.</gco:CharacterString>\n+ </gmd:title>\n+ <gmd:date>\n+ <gmd:CI_Date>\n+ <gmd:date>\n+ <gco:Date>2013-10-21</gco:Date>\n+ </gmd:date>\n+ <gmd:dateType>\n+ <gmd:CI_DateTypeCode codeList="http://standards.iso.org/iso/19139/resources/gmxCodelists.xml#CI_DateTypeCode"\n+ codeListValue="publication"/>\n+ </gmd:dateType>\n+ </gmd:CI_Date>\n+ </gmd:date>\n+ </gmd:CI_Citation>\n+ </gmd:specification>\n+ <gmd:explanation gco:nilReason="missing">\n+ <gco:CharacterString/>\n+ </gmd:explanation>\n+ <gmd:pass>\n+ <gco:Boolean>true</gco:Boolean>\n+ </gmd:pass>\n+ </gmd:DQ_ConformanceResult>\n+ </gmd:result>\n+ </gmd:DQ_DomainConsistency>\n+ </gmd:report>\n+ <gmd:lineage>\n+ <gmd:LI_Lineage>\n+ <gmd:statement>\n+ <gco:CharacterString>Relev\xc3\xa9 sur le terrain des adventices du ma\xc3\xafs</gco:CharacterString>\n+ </gmd:statement>\n+ </gmd:LI_Lineage>\n+ </gmd:lineage>\n+ </gmd:DQ_DataQuality>\n+ </gmd:dataQualityInfo>\n+</gmd:MD_Metadata>\n\\ No newline at end of file\n' |
| b |
| diff -r 000000000000 -r cf69ad260611 test-data/Metadata_validation.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/Metadata_validation.txt Mon Jan 09 13:36:02 2023 +0000 |
| b |
| @@ -0,0 +1,5 @@ + +Validation of metadata according to ISO 19139 + +TRUE + according to ISO 19139 XML schemas! |
| b |
| diff -r 000000000000 -r cf69ad260611 test-data/Mission.csv --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/Mission.csv Mon Jan 09 13:36:02 2023 +0000 |
| b |
| @@ -0,0 +1,2 @@ +"","x" +"1","SAFE" |
| b |
| diff -r 000000000000 -r cf69ad260611 test-data/NDVI.tabular --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/NDVI.tabular Mon Jan 09 13:36:02 2023 +0000 |
| b |
| b'@@ -0,0 +1,10001 @@\n+S2A_Subset\tlongitude\tlatitude\tNDVI\n+234.29\t13.7115065213769\t3.17883791054367\t0.809475381273922\n+233.06\t13.7124063428304\t3.17883903826377\t0.814533326859802\n+231.57\t13.7133061646046\t3.17884016519658\t0.808063813546791\n+228.73\t13.7142059866993\t3.17884129134209\t0.814678963045019\n+230.6\t13.7151058091144\t3.17884241670032\t0.812649359747791\n+228.85\t13.7160056318495\t3.17884354127124\t0.806605795437084\n+229.59\t13.7169054549044\t3.17884466505487\t0.810899033966968\n+239.74\t13.7178052782789\t3.1788457880512\t0.808152735556578\n+228.84\t13.7187051019728\t3.17884691026023\t0.81364672485546\n+217.93\t13.7196049259859\t3.17884803168196\t0.821532180501305\n+221.14\t13.720504750318\t3.17884915231639\t0.813394469354217\n+226.34\t13.7214045749687\t3.17885027216352\t0.804053528518387\n+233.05\t13.722304399938\t3.17885139122333\t0.802659325445255\n+237.25\t13.7232042252255\t3.17885250949585\t0.809442899205139\n+234.61\t13.7241040508311\t3.17885362698105\t0.811580900151433\n+225.7\t13.7250038767545\t3.17885474367895\t0.820858536885687\n+243.54\t13.7259037029954\t3.17885585958953\t0.805156185849842\n+238.98\t13.7268035295538\t3.17885697471281\t0.814407496583588\n+224.3\t13.7277033564293\t3.17885808904877\t0.806504430597208\n+229.94\t13.7286031836218\t3.17885920259741\t0.814004957858205\n+230.32\t13.7295030111309\t3.17886031535874\t0.834195663335835\n+221.94\t13.7304028389565\t3.17886142733275\t0.830148226182174\n+221.75\t13.7313026670984\t3.17886253851944\t0.817970615615008\n+231.52\t13.7322024955563\t3.17886364891882\t0.826119282944004\n+230.78\t13.73310232433\t3.17886475853087\t0.837485155719081\n+237.15\t13.7340021534193\t3.1788658673556\t0.807358464103723\n+240.64\t13.7349019828239\t3.178866975393\t0.834367581945824\n+240.36\t13.7358018125437\t3.17886808264308\t0.814613959603203\n+234.42\t13.7367016425784\t3.17886918910583\t0.823334743713733\n+232.98\t13.7376014729278\t3.17887029478125\t0.812603727061501\n+223.55\t13.7385013035916\t3.17887139966934\t0.817053245524247\n+228.45\t13.7394011345697\t3.1788725037701\t0.820890651769377\n+234.76\t13.7403009658618\t3.17887360708353\t0.81472573939151\n+241.63\t13.7412007974676\t3.17887470960963\t0.800657983613983\n+242.54\t13.7421006293871\t3.17887581134839\t0.804551773444016\n+233.03\t13.7430004616199\t3.17887691229981\t0.818208709568542\n+234.33\t13.7439002941658\t3.1788780124639\t0.828198556142271\n+224.33\t13.7448001270246\t3.17887911184064\t0.815919459122387\n+219.08\t13.7456999601961\t3.17888021043005\t0.813250047953031\n+223.36\t13.74659979368\t3.17888130823211\t0.809531707397978\n+231.08\t13.7474996274761\t3.17888240524683\t0.806462838677604\n+228.14\t13.7483994615842\t3.17888350147421\t0.822623856627664\n+223.55\t13.7492992960042\t3.17888459691424\t0.824117178204669\n+235.86\t13.7501991307356\t3.17888569156692\t0.800303742692492\n+241.66\t13.7510989657784\t3.17888678543225\t0.803205534643691\n+242.49\t13.7519988011323\t3.17888787851024\t0.824298525547548\n+241.99\t13.752898636797\t3.17888897080087\t0.83021872061821\n+239.49\t13.7537984727725\t3.17889006230415\t0.824182760846283\n+231.58\t13.7546983090583\t3.17889115302007\t0.812953085714065\n+239.36\t13.7555981456544\t3.17889224294864\t0.7917462194466\n+247.4\t13.7564979825604\t3.17889333208985\t0.802430972507897\n+242.26\t13.7573978197762\t3.17889442044371\t0.814116487206741\n+247.77\t13.7582976573015\t3.1788955080102\t0.816578431872301\n+247.25\t13.7591974951362\t3.17889659478933\t0.819701083956995\n+244.15\t13.7600973332799\t3.17889768078111\t0.808361251426994\n+235.55\t13.7609971717325\t3.17889876598551\t0.814163958835801\n+236.54\t13.7618970104937\t3.17889985040256\t0.817860100896017\n+259.75\t13.7627968495633\t3.17890093403223\t0.795576202131322\n+244.36\t13.7636966889412\t3.17890201687454\t0.81506586995564\n+235.82\t13.764596528627\t3.17890309892948\t0.815240243897869\n+233.5\t13.7654963686205\t3.17890418019705\t0.828422528097906\n+236.56\t13.7663962089216\t3.17890526067725\t0.81758530183727\n+235.9\t13.7672960495299\t3.17890634037008\t0.821981600526329\n+228.58\t13.7681958904453\t3.17890741927553\t0.812693877551021\n+230.77\t13.7690957316676\t3.17890849739361\t0.814535088866951\n+238.03\t13.7699955731965\t3.1789095747243\t0.815205860065673\n+249.26\t13.7'..b'26547\n+252.93\t13.7413078109934\t3.08933075095534\t0.799096191091027\n+240.39\t13.7422075664652\t3.08933182160117\t0.795246724048254\n+233.2\t13.7431073222502\t3.08933289148189\t0.799465229488297\n+243.09\t13.7440070783481\t3.08933396059749\t0.799829666089557\n+245.05\t13.7449068347588\t3.08933502894797\t0.801686929850912\n+244.84\t13.7458065914821\t3.08933609653333\t0.803598670211712\n+238.05\t13.7467063485177\t3.08933716335356\t0.805478783937571\n+243.54\t13.7476061058653\t3.08933822940867\t0.808654736735149\n+246.47\t13.7485058635249\t3.08933929469866\t0.803937941383833\n+236.57\t13.749405621496\t3.08934035922353\t0.795263388114729\n+238.45\t13.7503053797786\t3.08934142298326\t0.799478301322523\n+240.44\t13.7512051383724\t3.08934248597787\t0.812258322247346\n+241.73\t13.7521048972772\t3.08934354820734\t0.802059113928082\n+240.44\t13.7530046564927\t3.08934460967169\t0.816870644261308\n+240.1\t13.7539044160187\t3.0893456703709\t0.819512410175021\n+243.78\t13.754804175855\t3.08934673030498\t0.800793870580633\n+235.42\t13.7557039360014\t3.08934778947393\t0.806453695521473\n+244.45\t13.7566036964577\t3.08934884787774\t0.802636719930074\n+247.95\t13.7575034572236\t3.08934990551641\t0.801756528819869\n+244.46\t13.7584032182988\t3.08935096238994\t0.797846271081076\n+248.5\t13.7593029796833\t3.08935201849833\t0.804278812974465\n+244.12\t13.7602027413767\t3.08935307384158\t0.808916346373018\n+238.45\t13.7611025033789\t3.08935412841969\t0.801174579218478\n+240.39\t13.7620022656895\t3.08935518223265\t0.81026928149062\n+241.73\t13.7629020283085\t3.08935623528047\t0.80523587309033\n+247.45\t13.7638017912355\t3.08935728756314\t0.8062458918902\n+240.04\t13.7647015544703\t3.08935833908066\t0.809991343746832\n+239.59\t13.7656013180128\t3.08935938983303\t0.808963465088155\n+243.88\t13.7665010818626\t3.08936043982025\t0.798495321622938\n+256.03\t13.7674008460196\t3.08936148904233\t0.79006788195551\n+243.66\t13.7683006104835\t3.08936253749924\t0.789433387636698\n+242.63\t13.7692003752542\t3.08936358519101\t0.809001174539713\n+240.7\t13.7701001403313\t3.08936463211761\t0.799121705889283\n+242.53\t13.7709999057148\t3.08936567827906\t0.798083990686864\n+248.27\t13.7718996714042\t3.08936672367535\t0.810274370230963\n+251.33\t13.7727994373995\t3.08936776830648\t0.802474215384129\n+232.11\t13.7736992037004\t3.08936881217245\t0.808259587020649\n+232.32\t13.7745989703066\t3.08936985527326\t0.813367121187447\n+235.12\t13.7754987372181\t3.08937089760891\t0.814957466918715\n+244.63\t13.7763985044344\t3.08937193917939\t0.808755561484318\n+236.07\t13.7772982719554\t3.0893729799847\t0.824709811248584\n+237.69\t13.7781980397809\t3.08937402002485\t0.811220648031232\n+235.16\t13.7790978079107\t3.08937505929982\t0.819760336359234\n+237.52\t13.7799975763445\t3.08937609780963\t0.815774586678555\n+234.72\t13.7808973450821\t3.08937713555427\t0.822390289793584\n+239.01\t13.7817971141232\t3.08937817253373\t0.81059286697826\n+237.01\t13.7826968834678\t3.08937920874802\t0.827556259697598\n+231\t13.7835966531154\t3.08938024419714\t0.814163760971504\n+238.92\t13.784496423066\t3.08938127888108\t0.809629840125502\n+243.67\t13.7853961933193\t3.08938231279984\t0.818684700316627\n+237.8\t13.786295963875\t3.08938334595342\t0.821450619607607\n+239.14\t13.787195734733\t3.08938437834182\t0.813330233279082\n+238.26\t13.7880955058929\t3.08938540996504\t0.812077927230885\n+229.87\t13.7889952773547\t3.08938644082308\t0.822445096923031\n+228.06\t13.789895049118\t3.08938747091593\t0.825128780403618\n+235.08\t13.7907948211827\t3.0893885002436\t0.819054446666415\n+232.1\t13.7916945935485\t3.08938952880608\t0.825441298034367\n+242.5\t13.7925943662152\t3.08939055660337\t0.807354815933473\n+251.48\t13.7934941391826\t3.08939158363548\t0.792630563636647\n+236.59\t13.7943939124504\t3.08939260990239\t0.805452978341491\n+235.81\t13.7952936860184\t3.08939363540411\t0.810522000654324\n+221.18\t13.7961934598864\t3.08939466014064\t0.81823117345099\n+234.21\t13.7970932340542\t3.08939568411198\t0.814564687668025\n+235.43\t13.7979930085216\t3.08939670731812\t0.807905486583901\n+234.22\t13.7988927832882\t3.08939772975906\t0.809455691014441\n+246.08\t13.799792558354\t3.08939875143481\t0.808684324455036\n+251.94\t13.8006923337187\t3.08939977234536\t0.827205483763792\n' |
| b |
| diff -r 000000000000 -r cf69ad260611 test-data/S2A_MSIL2A_20200306T015621_N0214_R117_T51JXN_20200306T034744.zip |
| b |
| Binary file test-data/S2A_MSIL2A_20200306T015621_N0214_R117_T51JXN_20200306T034744.zip has changed |
| b |
| diff -r 000000000000 -r cf69ad260611 test-data/S2A_Subset |
| b |
| Binary file test-data/S2A_Subset has changed |
| b |
| diff -r 000000000000 -r cf69ad260611 test-data/S2A_Subset.hdr --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/S2A_Subset.hdr Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| @@ -0,0 +1,35 @@ +ENVI +description = { + File Resize Result, x resize factor: 1.000000, y resize factor: 1.000000. + [Thu Jul 11 09:00:28 2019]} +samples = 1000 +lines = 1000 +bands = 10 +header offset = 0 +file type = ENVI Standard +data type = 2 +interleave = bil +sensor type = Unknown +byte order = 0 +x start = 998 +y start = 766 +map info = {UTM, 1.000, 1.000, 356780.000, 351500.000, 1.0000000000e+001, 1.0000000000e+001, 33, North, WGS-84, units=Meters} +coordinate system string = {PROJCS["WGS_1984_UTM_Zone_33N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["false_easting",500000.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",15.0],PARAMETER["scale_factor",0.9996],PARAMETER["latitude_of_origin",0.0],UNIT["Meter",1.0]]} +wavelength units = Nanometers +z plot range = {-32768.00, 32767.00} +data ignore value = 0 +default stretch = 0.000000 1000.000000 linear +band names = { + Resize (band 02:S2A_T33NUD_20180104_Subset), + Resize (band 03:S2A_T33NUD_20180104_Subset), + Resize (band 04:S2A_T33NUD_20180104_Subset), + Resize (band 05:S2A_T33NUD_20180104_Subset), + Resize (band 06:S2A_T33NUD_20180104_Subset), + Resize (band 07:S2A_T33NUD_20180104_Subset), + Resize (band 08:S2A_T33NUD_20180104_Subset), + Resize (band 08A:S2A_T33NUD_20180104_Subset), + Resize (band 11:S2A_T33NUD_20180104_Subset), + Resize (band 12:S2A_T33NUD_20180104_Subset)} +wavelength = { + 496.600006, 560.000000, 664.500000, 703.900024, 740.200012, 782.500000, + 835.099976, 864.799988, 1613.699951, 2202.399902} |
| b |
| diff -r 000000000000 -r cf69ad260611 test-data/S2A_T33NUD_Plots.zip |
| b |
| Binary file test-data/S2A_T33NUD_Plots.zip has changed |
| b |
| diff -r 000000000000 -r cf69ad260611 val_metadata.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/val_metadata.r Mon Jan 09 13:36:02 2023 +0000 |
| [ |
| @@ -0,0 +1,41 @@ +#Rscript + +############################################ +## Validate ISO 19139 metadata documen ## +############################################ + +#####Packages : ncdf4, +# geometa, +# httr +# xml +# xml2 +library(geometa) + +#####Load arguments + +args <- commandArgs(trailingOnly = TRUE) + +if (length(args) < 1) { + stop("This tool needs at least 1 argument") +}else { + input_data <- args[1] +} + +##------------------------------------------## +## Read ISO 19139 from a file or url ## +##------------------------------------------## + +# Test depuis catalogue Indores http://indores-tmp.in2p3.fr/geonetwork/srv/fre/catalog.search#/metadata/112ebeea-e79c-422c-8a43-a5a8323b446b +# <!--ISO 19139 XML compliance: NO--> +input_data <- xml2::read_xml(input_data) + +dir.create("results") +file.create("results/meta.xml") + +xml2::write_xml(input_data, file = "results/meta.xml") + +md <- geometa::readISO19139("results/meta.xml") + + +# validate iso +cat("\nValidation of metadata according to ISO 19139\n", md$validate(), file = "Metadata_validation.txt", fill = 1, append = FALSE) |