Mercurial > repos > ecology > srs_process_data

comparison alpha_beta.r @ 0:82ee94fbefd7 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:53:22 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:82ee94fbefd7
1 #Rscript
2
3 ###########################################
4 ## Mapping alpha and beta diversity ##
5 ###########################################
6
7 #####Packages : stars
8 # utils
9 # biodivmapr
10 # raster
11 # sf
12 # mapview
13 # leafpop
14 # RColorBrewer
15 # labdsv
16 # rgdal
17 # ggplot2
18 # gridExtra
19
20 #####Load arguments
21
22 args <- commandArgs(trailingOnly = TRUE)
23
24 #####Import the S2 data
25
26 if (length(args) < 1) {
27 stop("This tool needs at least 1 argument")
28 }else {
29 data_raster <- args[1]
30 rasterheader <- args[2]
31 data <- args[3]
32 # type of PCA:
33 # PCA: no rescaling of the data
34 # SPCA: rescaling of the data
35 typepca <- as.character(args[4])
36 alpha <- as.logical(args[5])
37 beta <- as.logical(args[6])
38 funct <- as.logical(args[7])
39 all <- as.logical(args[8])
40 nbcpu <- as.integer(args[9])
41 source(args[10])
42 }
43
44 ################################################################################
45 ## DEFINE PARAMETERS FOR DATASET TO BE PROCESSED ##
46 ################################################################################
47 if (data_raster == "") {
48 #Create a directory where to unzip your folder of data
49 dir.create("data_dir")
50 unzip(data, exdir = "data_dir")
51 # Path to raster
52 data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl")
53 input_image_file <- file.path("data_dir/results/Reflectance", data_raster[1])
54 input_header_file <- file.path("data_dir/results/Reflectance", data_raster[2])
55
56 } else {
57 input_image_file <- file.path(getwd(), data_raster, fsep = "/")
58 input_header_file <- file.path(getwd(), rasterheader, fsep = "/")
59 }
60
61 ################################################################################
62 ## PROCESS IMAGE ##
63 ################################################################################
64 # 1- Filter data in order to discard non vegetated / shaded / cloudy pixels
65
66 print("PERFORM PCA ON RASTER")
67 pca_output <- biodivMapR::perform_PCA(Input_Image_File = input_image_file, Input_Mask_File = input_mask_file,
68 Output_Dir = output_dir, TypePCA = typepca, FilterPCA = filterpca, nbCPU = nbcpu, MaxRAM = maxram)
69
70 pca_files <- pca_output$PCA_Files
71 pix_per_partition <- pca_output$Pix_Per_Partition
72 nb_partitions <- pca_output$nb_partitions
73 # path for the updated mask
74 input_mask_file <- pca_output$MaskPath
75
76
77 selected_pcs <- seq(1, dim(raster::stack(input_image_file))[3])
78
79 selected_pcs <- all(selected_pcs)
80 ################################################################################
81 ## MAP ALPHA AND BETA DIVERSITY ##
82 ################################################################################
83 print("MAP SPECTRAL SPECIES")
84
85 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)
86
87 image_name <- tools::file_path_sans_ext(basename(input_image_file))
88 if (alpha == TRUE || beta == TRUE || all == TRUE) {
89 ## alpha
90 print("MAP ALPHA DIVERSITY")
91 index_alpha <- c("Shannon")
92 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)
93
94 alpha_zip <- file.path(output_dir, image_name, typepca, "ALPHA", "Shannon_10_Fullres.zip")
95 alpha_path <- file.path(output_dir, image_name, typepca, "ALPHA")
96 unzip(alpha_zip, exdir = alpha_path)
97 alpha_path <- file.path(output_dir, image_name, typepca, "ALPHA", "Shannon_10_Fullres")
98 alpha_raster <- raster::raster(alpha_path)
99 get_alpha <- convert_raster(alpha_raster)
100
101 if (alpha == TRUE || all == TRUE) {
102 colnames(get_alpha) <- c("Alpha", "longitude", "latitude")
103 plot_indices(get_alpha, titre = "Alpha")
104
105 write.table(get_alpha, file = "alpha.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
106 }
107 if (beta == TRUE || all == TRUE) {
108 ## beta
109 print("MAP BETA DIVERSITY")
110 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)
111
112 beta_path <- file.path(output_dir, image_name, typepca, "BETA", "BetaDiversity_BCdiss_PCO_10")
113 beta_raster <- raster::raster(beta_path)
114 get_beta <- convert_raster(beta_raster)
115
116 colnames(get_beta) <- c("Beta", "longitude", "latitude")
117 plot_indices(get_beta, titre = "Beta")
118
119 write.table(get_beta, file = "beta.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
120 }
121 }
122
123
124 ################################################################################
125 ## COMPUTE ALPHA AND BETA DIVERSITY FROM FIELD PLOTS ##
126 ################################################################################
127
128 if (funct == TRUE || all == TRUE) {
129 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)
130
131 funct_zip <- file.path(output_dir, image_name, typepca, "FUNCTIONAL", "FunctionalDiversity_Map_MeanFilter_Fullres.zip")
132 funct_path <- file.path(output_dir, image_name, typepca, "FUNCTIONAL")
133 unzip(funct_zip, exdir = funct_path)
134 funct_path <- file.path(output_dir, image_name, typepca, "FUNCTIONAL", "FunctionalDiversity_Map_MeanFilter_Fullres")
135 funct_raster <- raster::raster(funct_path)
136 get_funct <- convert_raster(funct_raster)
137
138 colnames(get_funct) <- c("Functionnal", "longitude", "latitude")
139 plot_indices(get_funct, titre = "Functionnal")
140
141 write.table(get_funct, file = "Functionnal.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
142 }