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