Mercurial > repos > ecology > srs_preprocess_s2

comparison comparison_div.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 : stars
8 # utils
9 # biodivmapr
10 # raster
11 # sf
12 # mapview
13 # leafpop
14 # RColorBrewer
15 # labdsv
16 # rgdal
17 # ggplot2
18 # gridExtra
19 ##remotes::install_github("jbferet/biodivMapR")
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 plots_zip <- args[4]
33 choice <- as.character(args[5])
34 source(args[6])
35 # type of PCA:
36 # PCA: no rescaling of the data
37 # SPCA: rescaling of the data
38 typepca <- as.character(args[7])
39 }
40
41 ################################################################################
42 ## DEFINE PARAMETERS FOR DATASET TO BE PROCESSED ##
43 ################################################################################
44 if (data_raster == "") {
45 #Create a directory where to unzip your folder of data
46 dir.create("data_dir")
47 unzip(data, exdir = "data_dir")
48 # Path to raster
49 data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl")
50 input_image_file <- file.path("data_dir/results/Reflectance", data_raster[1])
51 input_header_file <- file.path("data_dir/results/Reflectance", data_raster[2])
52
53 } else {
54 input_image_file <- file.path(getwd(), data_raster, fsep = "/")
55 input_header_file <- file.path(getwd(), rasterheader, fsep = "/")
56 }
57
58 ################################################################################
59 ## PROCESS IMAGE ##
60 ################################################################################
61 # 1- Filter data in order to discard non vegetated / shaded / cloudy pixels
62
63 print("PERFORM PCA ON RASTER")
64 pca_output <- biodivMapR::perform_PCA(Input_Image_File = input_image_file, Input_Mask_File = input_mask_file,
65 Output_Dir = output_dir, TypePCA = typepca, FilterPCA = filterpca, nbCPU = nbcpu, MaxRAM = maxram)
66
67 pca_files <- pca_output$PCA_Files
68 pix_per_partition <- pca_output$Pix_Per_Partition
69 nb_partitions <- pca_output$nb_partitions
70 # path for the updated mask
71 input_mask_file <- pca_output$MaskPath
72
73 # 3- Select principal components from the PCA raster
74 # Select components from the PCA/SPCA/MNF raster
75 sel_compo <- c("1\n", "2\n", "3\n", "4\n", "5\n", "6\n", "7\n", "8")
76 image_name <- tools::file_path_sans_ext(basename(input_image_file))
77 output_dir_full <- file.path(output_dir, image_name, typepca, "PCA")
78
79 write.table(sel_compo, paste0(output_dir_full, "/Selected_Components.txt"))
80 sel_pc <- file.path(output_dir_full, "Selected_Components.txt")
81
82
83 ################################################################################
84 ## MAP ALPHA AND BETA DIVERSITY ##
85 ################################################################################
86 print("MAP SPECTRAL SPECIES")
87
88 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)
89
90 ################################################################################
91 ## COMPUTE ALPHA AND BETA DIVERSITY FROM FIELD PLOTS ##
92 ################################################################################
93 ## read selected features from dimensionality reduction
94
95 ## path for selected components
96
97 # location of the directory where shapefiles used for validation are saved
98 dir.create("VectorDir")
99 unzip(plots_zip, exdir = "VectorDir")
100
101 # list vector data
102 path_vector <- biodivMapR::list_shp("VectorDir")
103 name_vector <- tools::file_path_sans_ext(basename(path_vector))
104
105 # location of the spectral species raster needed for validation
106 path_spectralspecies <- kmeans_info$SpectralSpecies
107 # get diversity indicators corresponding to shapefiles (no partitioning of spectral dibversity based on field plots so far...)
108 biodiv_indicators <- biodivMapR::diversity_from_plots(Raster_SpectralSpecies = path_spectralspecies, Plots = path_vector, nbclusters = nbclusters, Raster_Functional = pca_files, Selected_Features = FALSE)
109
110 shannon_rs <- c(biodiv_indicators$Shannon)[[1]]
111 fric <- c(biodiv_indicators$FunctionalDiversity$FRic)
112 feve <- c(biodiv_indicators$FunctionalDiversity$FEve)
113 fdiv <- c(biodiv_indicators$FunctionalDiversity$FDiv)
114 # if no name for plots
115 biodiv_indicators$Name_Plot <- seq(1, length(biodiv_indicators$Shannon[[1]]), by = 1)
116
117
118 ####################################################
119 # write RS indicators #
120 ####################################################
121 # write a table for Shannon index
122
123 # write a table for all spectral diversity indices corresponding to alpha diversity
124 results <- data.frame(name_vector, biodiv_indicators$Richness, biodiv_indicators$Fisher,
125 biodiv_indicators$Shannon, biodiv_indicators$Simpson,
126 biodiv_indicators$FunctionalDiversity$FRic,
127 biodiv_indicators$FunctionalDiversity$FEve,
128 biodiv_indicators$FunctionalDiversity$FDiv)
129
130 names(results) <- c("ID_Plot", "Species_Richness", "Fisher", "Shannon", "Simpson", "fric", "feve", "fdiv")
131 write.table(results, file = "Diversity.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
132
133 if (choice == "Y") {
134 # write a table for Bray Curtis dissimilarity
135 bc_mean <- biodiv_indicators$BCdiss
136 bray_curtis <- data.frame(name_vector, bc_mean)
137 colnames(bray_curtis) <- c("ID_Plot", bray_curtis[, 1])
138 write.table(bray_curtis, file = "BrayCurtis.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
139
140 ####################################################
141 # illustrate results
142 ####################################################
143 # apply ordination using PCoA (same as done for map_beta_div)
144
145 mat_bc_dist <- as.dist(bc_mean, diag = FALSE, upper = FALSE)
146 betapco <- labdsv::pco(mat_bc_dist, k = 3)
147
148 # assign a type of vegetation to each plot, assuming that the type of vegetation
149 # is defined by the name of the shapefile
150
151 nbsamples <- shpname <- c()
152 for (i in 1:length(path_vector)) {
153 shp <- path_vector[i]
154 nbsamples[i] <- length(rgdal::readOGR(shp, verbose = FALSE))
155 shpname[i] <- tools::file_path_sans_ext(basename(shp))
156 }
157
158 type_vegetation <- c()
159 for (i in 1: length(nbsamples)) {
160 for (j in 1:nbsamples[i]) {
161 type_vegetation <- c(type_vegetation, shpname[i])
162 }
163 }
164
165 #data frame including a selection of alpha diversity metrics and beta diversity expressed as coordinates in the PCoA space
166 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)
167
168 #plot field data in the PCoA space, with size corresponding to shannon index
169 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"))
170
171 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"))
172
173 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"))
174
175 #extract legend
176 get_legend <- function(a_gplot) {
177 tmp <- ggplot2::ggplot_gtable(ggplot2::ggplot_build(a_gplot))
178 leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
179 legend <- tmp$grobs[[leg]]
180 return(legend)
181 }
182
183 legend <- get_legend(g3)
184 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))
185
186
187 filename <- ggplot2::ggsave("BetaDiversity_PcoA1_vs_PcoA2_vs_PcoA3.png", gall, scale = 0.65, width = 12, height = 9, units = "in", dpi = 200, limitsize = TRUE)
188
189 filename
190 }