Mercurial > repos > ecology > srs_process_data
comparison biodiv_indices_global.r @ 0:cf69ad260611 draft default tip
planemo upload for repository https://github.com/Marie59/Sentinel_2A/srs_tools commit b32737c1642aa02cc672534e42c5cb4abe0cd3e7
author | ecology |
---|---|
date | Mon, 09 Jan 2023 13:36:02 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:cf69ad260611 |
---|---|
1 #Rscript | |
2 | |
3 ########################################### | |
4 ## Mapping biodiversity indicators ## | |
5 ########################################### | |
6 | |
7 #####Packages : raster, | |
8 # rgdal, | |
9 # sp, | |
10 # rasterdiv, | |
11 # ggplot2, | |
12 | |
13 #####Load arguments | |
14 | |
15 args <- commandArgs(trailingOnly = TRUE) | |
16 | |
17 #####Import the S2 data | |
18 | |
19 if (length(args) < 1) { | |
20 stop("This tool needs at least 1 argument") | |
21 }else { | |
22 data_raster <- args[1] | |
23 data_header <- args[2] | |
24 data <- args[3] | |
25 alpha <- as.numeric(args[4]) | |
26 source(args[5]) | |
27 | |
28 } | |
29 | |
30 ######################################################################## | |
31 ## COMPUTE BIODIVERSITY INDICES ## | |
32 ######################################################################## | |
33 | |
34 if (data_raster == "") { | |
35 #Create a directory where to unzip your folder of data | |
36 dir.create("data_dir") | |
37 unzip(data, exdir = "data_dir") | |
38 # Path to raster | |
39 data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl") | |
40 data_raster <- file.path("data_dir/results/Reflectance", data_raster[1]) | |
41 } | |
42 | |
43 # Read raster | |
44 copndvi <- raster::raster(data_raster) | |
45 | |
46 copndvilr <- raster::reclassify(copndvi, cbind(252, 255, NA), right = TRUE) | |
47 | |
48 #Resample using raster::aggregate and a linear factor of 10 | |
49 copndvilr <- raster::aggregate(copndvilr, fact = 20) | |
50 #Set float numbers as integers to further speed up the calculation | |
51 storage.mode(copndvilr[]) <- "integer" | |
52 | |
53 # Convert raster to SpatialPointsDataFrame | |
54 r_pts <- convert_raster(copndvilr) | |
55 | |
56 #Shannon's Diversity | |
57 sha <- rasterdiv::Shannon(copndvilr, window = 9, np = 1) | |
58 sha_df <- data.frame(raster::rasterToPoints(sha, spatial = TRUE)) | |
59 sha_name <- "Shannon" | |
60 r_pts[, sha_name] <- sha_df[, 1] | |
61 | |
62 #Renyi's Index | |
63 ren <- rasterdiv::Renyi(copndvilr, window = 9, alpha = alpha, np = 1) | |
64 ren_df <- data.frame(raster::rasterToPoints(ren[[1]])) | |
65 ren_name <- "Renyi" | |
66 r_pts[, ren_name] <- ren_df[, 3] | |
67 | |
68 #Berger-Parker's Index | |
69 ber <- rasterdiv::BergerParker(copndvilr, window = 9, np = 1) | |
70 ber_df <- data.frame(raster::rasterToPoints(ber, spatial = TRUE)) | |
71 ber_name <- "Berger-Parker" | |
72 r_pts[, ber_name] <- ber_df[, 1] | |
73 | |
74 #Pielou's Evenness | |
75 pie <- rasterdiv::Pielou(copndvilr, window = 9, np = 1) | |
76 pie_df <- data.frame(raster::rasterToPoints(pie)) | |
77 if (length(pie_df[, 1]) == length(r_pts[, 1])) { | |
78 pie_name <- "Pielou" | |
79 r_pts[, pie_name] <- pie_df[, 1] | |
80 } | |
81 | |
82 #Hill's numbers | |
83 hil <- rasterdiv::Hill(copndvilr, window = 9, alpha = alpha, np = 1) | |
84 hil_df <- data.frame(raster::rasterToPoints(hil[[1]])) | |
85 hil_name <- "Hill" | |
86 r_pts[, hil_name] <- hil_df[, 3] | |
87 | |
88 #Parametric Rao's quadratic entropy with alpha ranging from 1 to 5 | |
89 prao <- rasterdiv::paRao(copndvilr, window = 9, alpha = alpha, dist_m = "euclidean", np = 1) | |
90 prao_df <- data.frame(raster::rasterToPoints(prao$window.9[[1]])) | |
91 prao_name <- "Prao" | |
92 r_pts[, prao_name] <- prao_df[, 3] | |
93 | |
94 #Cumulative Residual Entropy | |
95 cre <- rasterdiv::CRE(copndvilr, window = 9, np = 1) | |
96 cre_df <- data.frame(raster::rasterToPoints(cre)) | |
97 if (length(cre_df[, 1]) == length(r_pts[, 1])) { | |
98 cre_name <- "CRE" | |
99 r_pts[, cre_name] <- cre_df[, 1] | |
100 } | |
101 | |
102 if (length(cre_df[, 1]) == length(r_pts[, 1]) || length(pie_df[, 1]) == length(r_pts[, 1])) { | |
103 list_indice <- list("Shannon", "Renyi", "Berger-Parker", "Pielou", "Hill", "Prao", "CRE") | |
104 } else { | |
105 list_indice <- list("Shannon", "Renyi", "Berger-Parker", "Hill", "Prao") | |
106 } | |
107 ## Plotting all the graph and writing a tabular | |
108 for (indice in list_indice) { | |
109 plot_indices(data = r_pts, titre = indice) | |
110 } | |
111 | |
112 write.table(r_pts, file = "BiodivIndex.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE) |