Mercurial > repos > ecology > srs_preprocess_s2

comparison biodiv_indices_global.r @ 0:33a1e15f7252 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:37:05 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:33a1e15f7252
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)