|
0
|
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
|
