Mercurial > repos > ecology > srs_process_data
diff indices_spectral.r @ 0:cf69ad260611 draft default tip
planemo upload for repository https://github.com/Marie59/Sentinel_2A/srs_tools commit b32737c1642aa02cc672534e42c5cb4abe0cd3e7
author | ecology |
---|---|
date | Mon, 09 Jan 2023 13:36:02 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/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