Mercurial > repos > ecology > srs_preprocess_s2

comparison indices_spectral.r @ 0:a990ea26442f draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/Marie59/Sentinel_2A/srs_tools commit b32737c1642aa02cc672534e42c5cb4abe0cd3e7
author ecology
date Mon, 09 Jan 2023 13:54:29 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:a990ea26442f
1 #Rscript
2
3 ###########################################
4 ## Mapping alpha and beta diversity ##
5 ###########################################
6
7 #####Packages : expint,
8 # pracma,
9 # R.utils,
10 # raster,
11 # sp,
12 # matrixStats,
13 # ggplot2,
14 # expandFunctions,
15 # stringr,
16 # XML,
17 # rgdal,
18 # stars,
19 #####Load arguments
20
21 args <- commandArgs(trailingOnly = TRUE)
22
23 #####Import the S2 data
24
25 if (length(args) < 1) {
26 stop("This tool needs at least 1 argument")
27 }else {
28 data_raster <- args[1]
29 data_header <- args[2]
30 data <- args[3]
31 source(args[4])
32 source(args[5])
33 source(args[6])
34 indice_choice <- strsplit(args[7], ",")[[1]]
35 source(args[8])
36 output_raster <- as.character(args[9])
37
38 }
39
40 ########################################################################
41 ## COMPUTE SPECTRAL INDEX : NDVI ##
42 ########################################################################
43
44 if (data != "") {
45 #Create a directory where to unzip your folder of data
46 dir.create("data_dir")
47 unzip(data, exdir = "data_dir")
48
49 # Read raster
50 data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl")
51 data_raster <- file.path("data_dir/results/Reflectance", data_raster[1])
52 refl <- raster::brick(data_raster)
53 refl2 <- raster::raster(data_raster)
54 } else {
55 # Read raster
56 refl <- raster::brick(data_raster)
57 refl2 <- raster::raster(data_raster)
58 }
59 # get raster band name and clean format. Expecting band name and wav
60 # reflFactor = 10000 when reflectance is coded as INT16
61 refl <- raster::aggregate(refl, fact = 10)
62
63 # Convert raster to SpatialPointsDataFrame
64 refl2 <- raster::aggregate(refl2, fact = 10)
65 r_pts <- convert_raster(refl2)
66 table_ind <- r_pts
67 # create directory for Spectralelength to be documented in image
68 hdr_refl <- read_envi_header(get_hdr_name(data_raster))
69 sensorbands <- hdr_refl$wavelength
70 # compute a set of spectral indices defined by indexlist from S2 data indices
71 si_path <- file.path("SpectralIndices")
72 dir.create(path = si_path, showWarnings = FALSE, recursive = TRUE)
73 # Save spectral indices
74
75 indices <- lapply(indice_choice, function(x) {
76 indices_list <- computespectralindices_raster(refl = refl, sensorbands = sensorbands,
77 sel_indices = x,
78 reflfactor = 10000, stackout = FALSE)
79
80 index_path <- file.path(si_path, paste(basename(data_raster), "_", x, sep = ""))
81 spec_indices <- stars::write_stars(stars::st_as_stars(indices_list$spectralindices[[1]]), dsn = index_path, driver = "ENVI", type = "Float32")
82
83 # write band name in HDR
84 hdr <- read_envi_header(get_hdr_name(index_path))
85 hdr$`band names` <- x
86 write_envi_header(hdr = hdr, hdrpath = get_hdr_name(index_path))
87 # Writting the tabular and the plot
88 r_pts[, x] <- as.data.frame(sapply(spec_indices, c))
89 plot_indices(data = r_pts, titre = x)
90 return(r_pts)
91 })
92
93 new_table <- as.data.frame(indices)
94 new_table <- new_table[, !grepl("longitude", names(new_table))]
95 new_table <- new_table[, !grepl("latitude", names(new_table))]
96 new_table <- new_table[, !grepl(basename(data_raster), names(new_table))]
97
98 table_ind <- cbind(table_ind, new_table)
99 if (length(indice_choice) == 1) {
100 colnames(table_ind) <- c(basename(data_raster), "longitude", "latitude", indice_choice)
101 }
102
103 write.table(table_ind, file = "Spec_Index.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
104
105 # Get the raster layer of the indice as an output