|
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 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 val_metadata.xml |
| b |
| diff -r 000000000000 -r 054b2522a933 Lib_preprocess_S2.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Lib_preprocess_S2.r Mon Jan 09 13:38:38 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 054b2522a933 alpha_beta.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/alpha_beta.r Mon Jan 09 13:38:38 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 054b2522a933 biodiv_indices_global.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/biodiv_indices_global.r Mon Jan 09 13:38:38 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 054b2522a933 comparison_div.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/comparison_div.r Mon Jan 09 13:38:38 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 054b2522a933 functions.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/functions.r Mon Jan 09 13:38:38 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 054b2522a933 indices_spectral.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/indices_spectral.r Mon Jan 09 13:38:38 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 054b2522a933 macro.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/macro.xml Mon Jan 09 13:38:38 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 054b2522a933 pca_raster.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pca_raster.r Mon Jan 09 13:38:38 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 054b2522a933 preprocess_S2.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/preprocess_S2.r Mon Jan 09 13:38:38 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 054b2522a933 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:38:38 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 054b2522a933 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:38:38 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 054b2522a933 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:38:38 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 054b2522a933 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:38:38 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 054b2522a933 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:38:38 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 054b2522a933 test-data/Mission.csv --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/Mission.csv Mon Jan 09 13:38:38 2023 +0000 |
| b |
| @@ -0,0 +1,2 @@ +"","x" +"1","SAFE" |
| b |
| diff -r 000000000000 -r 054b2522a933 test-data/NDVI.tabular --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/NDVI.tabular Mon Jan 09 13:38:38 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 054b2522a933 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 054b2522a933 test-data/S2A_Subset |
| b |
| Binary file test-data/S2A_Subset has changed |
| b |
| diff -r 000000000000 -r 054b2522a933 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:38:38 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 054b2522a933 test-data/S2A_T33NUD_Plots.zip |
| b |
| Binary file test-data/S2A_T33NUD_Plots.zip has changed |
| b |
| diff -r 000000000000 -r 054b2522a933 val_metadata.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/val_metadata.r Mon Jan 09 13:38:38 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) |
| b |
| diff -r 000000000000 -r 054b2522a933 val_metadata.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/val_metadata.xml Mon Jan 09 13:38:38 2023 +0000 |
| [ |
| @@ -0,0 +1,67 @@ +<tool id="srs_metadata" name="Validate ISO 19139" version="@VERSION@" profile="20.01"> + <description>metadata documents</description> + <macros> + <import>macro.xml</import> + </macros> + <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="3.5_21">r-raster</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.20">r-ncdf4</requirement> + <requirement type="package" version="0.7">r-geometa</requirement> + <requirement type="package" version="1.4.4">r-httr</requirement> + <requirement type="package" version="3.99_0.13">r-xml</requirement> + <requirement type="package" version="1.3.3">r-xml2</requirement> + </requirements> + <command detect_errors="exit_code"><![CDATA[ + Rscript + '$__tool_directory__/val_metadata.r' + '$input' + '$output' 1> Metadata_validation.txt + ]]> + </command> + <inputs> + <param name="input" type="data" format="xml" label="XML file to validate"/> + </inputs> + <outputs> + <data name="output" from_work_dir="Metadata_validation.txt" format="txt" label="Validate metadata"/> + </outputs> + <tests> + <test> + <param name="input" value="12a0b625-9ad5-4251-a57a-305e22edef2e.xml"/> + <output name="output" value="Metadata_validation.txt"/> + </test> + </tests> + <help><![CDATA[ +========================== +Preprocess Sentinel 2 data +========================== + + +**What it does** + +Read a ISO 19139 from a file or url into an object in the geometa model. + +**Input description** + +A metadata xml file or url. + ++-------------+ +| Metadata | ++=============+ +| file or url | ++-------------+ +| ... | ++-------------+ + +**Output** + +A text file saying TRUE or FALSE depending on if the xml metadata are conforme according to ISO 19139. If the xml metadata is not conform it will write down where are the issues. + ]]> </help> + <expand macro="SRS_metaref"/> +</tool> |