Mercurial > repos > ecology > cb_ivr

comparison cb_ivr.r @ 0:89f977a3a256 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/Marie59/champ_blocs commit 0d86db7d42b608c386a54500064f5f9c9d7019a4
author ecology
date Wed, 04 Jan 2023 13:20:53 +0000
parents
children 83a4d3b47184
comparison
equal deleted inserted replaced
-1:000000000000 0:89f977a3a256
1 # author: "Jonathan Richir"
2 # date: "19 April 2021" )
3
4 #Rscript
5
6 ###############################
7 ## ##
8 ###############################
9
10 #####Packages : dplyr
11 # tidyr
12 # readr
13 # writexl
14 # stringr
15 # readxl
16 # tibble
17 # lubridate
18 # cowplot
19 # magrittr
20 # rmarkdown
21 library(magrittr)
22 #####Load arguments
23
24 args <- commandArgs(trailingOnly = TRUE)
25
26 ### Import data
27
28 if (length(args) < 1) {
29 stop("This tool needs at least 1 argument")
30 }else {
31 fiche_val <- args[1]
32 input_data <- args[2]
33
34 }
35
36 #############################################################
37 # #
38 # Load and clean the data #
39 # #
40 #############################################################
41 ### load ivr data
42
43 ivr <- read.csv2(input_data, header = FALSE, fileEncoding = "Latin1")
44 names_ <- as.vector(unlist(ivr[1, ]))
45 names_ <- gsub(" ", ".", names_)
46 colnames(ivr) <- names_
47 ivr <- ivr[-1, ]
48 ivr <- ivr[, -17]
49
50 # NB inversion between id and ID.Fiche variable names
51 ivr <- dplyr::rename(ivr, XX = ID.Fiche)
52 ivr <- dplyr::rename(ivr, ID.Fiche = id)
53 ivr <- dplyr::rename(ivr, id = XX)
54
55
56 ### load excel files "Fiche terrain" the metadata
57
58 fiche <- read.csv2(fiche_val, fileEncoding = "Latin1") # fileEncoding = "Latin1" cfr é in variable names
59
60 date_fiche <- as.Date(stringr::str_sub(fiche$date.sortie, end = 10), origin = "1970-01-01")
61 fiche <- tibble::add_column(fiche, date_fiche, .after = "date.sortie")
62 rm(date_fiche)
63
64 ## ivr vs fiche terrain
65 ivr$id <- as.numeric(ivr[, c("id")])
66
67 fiche_red <- dplyr::filter(fiche, fiche$ID.Fiche %in% unique(ivr[, c("id")]))
68
69 id_count <- ivr %>% dplyr::group_by(id) %>% dplyr::count()
70 id_count <- dplyr::rename(id_count, "ID.Fiche" = "id")
71 id_count <- dplyr::ungroup(id_count)
72 id_count <- as.data.frame(id_count)
73
74 fiche_red <- dplyr::left_join(fiche_red, id_count)
75
76 # rep fiche terrain information
77 fiche_expanded <- fiche_red[rep(row.names(fiche_red), fiche_red$n), 1:ncol(fiche_red)]
78 fiche_expanded <- dplyr::rename(fiche_expanded, "id" = "ID.Fiche")
79
80 ## merge ivr data and ficheterrain information
81 ivr <- dplyr::bind_cols(ivr, fiche_expanded)
82 ivr <- dplyr::rename(ivr, "id.ivr" = "id...1")
83 ivr <- dplyr::rename(ivr, "id.fiche" = "id...17")
84
85 rm(fiche_expanded, fiche_red, id_count)
86
87 ivr <- ivr %>% tidyr::separate(date_fiche, c("Year", "Month", "Day"), sep = "-", remove = FALSE)
88
89 ## I create two new variables for Site names, one for data analysis and one for data reporting. Only works for actual ivr df with 22 sites !
90
91 # Name for data analysis
92 ivr <- tibble::add_column(ivr, Site = ivr$zone.habitat, .after = "ID.Fiche")
93 ivr$Site <- gsub(pattern = " \\(champ de blocs\\)", replacement = "", ivr$Site)
94 ivr$Site <- gsub(pattern = " \\(champ blocs\\)", replacement = "", ivr$Site)
95
96 for (x in seq_along(ivr$Site)) {
97 if (grepl(pattern = "Locmariaquer", ivr$Site[x]) == TRUE) {
98 ivr$Site[x] <- "GDMO_Locmariaquer"
99 } else if (grepl(pattern = "Beg Lann", ivr$Site[x]) == TRUE) {
100 ivr$Site[x] <- "GDMO_BegLann"
101 } else if (grepl(pattern = "Plateau du Four", ivr$Site[x]) == TRUE) {
102 ivr$Site[x] <- "FOUR_PlateauFour"
103 } else if (grepl(pattern = "Grouin", ivr$Site[x]) == TRUE) {
104 ivr$Site[x] <- "EGMP_GroinCou"
105 } else if (grepl(pattern = "Ensembert", ivr$Site[x]) == TRUE) {
106 ivr$Site[x] <- "EGMP_PasEmsembert"
107 } else if (grepl(pattern = "Brée-les-Bains", ivr$Site[x]) == TRUE) {
108 ivr$Site[x] <- "EGMP_BreeBains"
109 } else if (grepl(pattern = "Antiochat", ivr$Site[x]) == TRUE) {
110 ivr$Site[x] <- "EGMP_PerreAntiochat"
111 } else if (grepl(pattern = "Chassiron", ivr$Site[x]) == TRUE) {
112 ivr$Site[x] <- "EGMP_Chassiron"
113 } else if (grepl(pattern = "zone p", ivr$Site[x]) == TRUE) {
114 ivr$Site[x] <- "BASQ_FlotsBleusZP"
115 } else if (grepl(pattern = "zone f", ivr$Site[x]) == TRUE) {
116 ivr$Site[x] <- "BASQ_FlotsBleusZF"
117 } else if (grepl(pattern = "Saint-Michel", ivr$Site[x]) == TRUE) {
118 ivr$Site[x] <- "GONB_IlotStMichel"
119 } else if (grepl(pattern = "Quéménès", ivr$Site[x]) == TRUE) {
120 ivr$Site[x] <- "FINS_Quemenes"
121 } else if (grepl(pattern = "Goulenez", ivr$Site[x]) == TRUE) {
122 ivr$Site[x] <- "FINS_SeinGoulenez"
123 } else if (grepl(pattern = "Kilaourou", ivr$Site[x]) == TRUE) {
124 ivr$Site[x] <- "FINS_SeinKilaourou"
125 } else if (grepl(pattern = "Verdelet", ivr$Site[x]) == TRUE) {
126 ivr$Site[x] <- "ARMO_Verdelet"
127 } else if (grepl(pattern = "Piégu", ivr$Site[x]) == TRUE) {
128 ivr$Site[x] <- "ARMO_Piegu"
129 } else if (grepl(pattern = "Bilfot", ivr$Site[x]) == TRUE) {
130 ivr$Site[x] <- "ARMO_Bilfot"
131 } else if (grepl(pattern = "Plate", ivr$Site[x]) == TRUE) {
132 ivr$Site[x] <- "ARMO_IlePlate"
133 } else if (grepl(pattern = "Perharidy", ivr$Site[x]) == TRUE) {
134 ivr$Site[x] <- "PDMO_Perharidy"
135 } else if (grepl(pattern = "Keraliou", ivr$Site[x]) == TRUE) {
136 ivr$Site[x] <- "BRES_Keraliou"
137 } else if (grepl(pattern = "Mousterlin", ivr$Site[x]) == TRUE) {
138 ivr$Site[x] <- "FINS_Mousterlin"
139 } else if (grepl(pattern = "Nicolas", ivr$Site[x]) == TRUE) {
140 ivr$Site[x] <- "FINS_StNicolasGlenan"
141 }
142 if (grepl(pattern = "Roz", ivr$site[x]) == TRUE) {
143 ivr$Site[x] <- "FINS_AnseRoz"
144 }
145 }
146
147 # Name for report/plot
148
149 ivr <- tibble::add_column(ivr, Site_bis = ivr$Site, .after = "Site")
150
151 ivr$Site_bis <- ifelse(ivr$Site == "GDMO_Locmariaquer", "Locmariaquer", ivr$Site_bis)
152 ivr$Site_bis <- ifelse(ivr$Site == "GDMO_BegLann", "Beg Lann", ivr$Site_bis)
153 ivr$Site_bis <- ifelse(ivr$Site == "FOUR_PlateauFour", "Plateau du Four", ivr$Site_bis)
154 ivr$Site_bis <- ifelse(ivr$Site == "EGMP_GroinCou", "Grouin du Cou", ivr$Site_bis)
155 ivr$Site_bis <- ifelse(ivr$Site == "EGMP_PasEmsembert", "Le Pas d'Emsembert", ivr$Site_bis)
156 ivr$Site_bis <- ifelse(ivr$Site == "EGMP_BreeBains", "La Brée-les-Bains", ivr$Site_bis)
157 ivr$Site_bis <- ifelse(ivr$Site == "EGMP_PerreAntiochat", "Le Perré d'Antiochat", ivr$Site_bis)
158 ivr$Site_bis <- ifelse(ivr$Site == "EGMP_Chassiron", "Chassiron", ivr$Site_bis)
159 ivr$Site_bis <- ifelse(ivr$Site == "BASQ_FlotsBleusZP", "Les Flots Bleus / zone pêcheurs", ivr$Site_bis)
160 ivr$Site_bis <- ifelse(ivr$Site == "BASQ_FlotsBleusZF", "Les Flots Bleus / zone familles", ivr$Site_bis)
161 ivr$Site_bis <- ifelse(ivr$Site == "GONB_IlotStMichel", "Îlot Saint-Michel", ivr$Site_bis)
162 ivr$Site_bis <- ifelse(ivr$Site == "FINS_Quemenes", "Quéménès", ivr$Site_bis)
163 ivr$Site_bis <- ifelse(ivr$Site == "FINS_SeinGoulenez", "Île de Sein - Goulenez", ivr$Site_bis)
164 ivr$Site_bis <- ifelse(ivr$Site == "FINS_SeinKilaourou", "Île de Sein - Kilaourou", ivr$Site_bis)
165 ivr$Site_bis <- ifelse(ivr$Site == "ARMO_Verdelet", "Îlot du Verdelet", ivr$Site_bis)
166 ivr$Site_bis <- ifelse(ivr$Site == "ARMO_Piegu", "Piégu", ivr$Site_bis)
167 ivr$Site_bis <- ifelse(ivr$Site == "ARMO_Bilfot", "Pointe de Bilfot", ivr$Site_bis)
168 ivr$Site_bis <- ifelse(ivr$Site == "ARMO_IlePlate", "Île Plate", ivr$Site_bis)
169 ivr$Site_bis <- ifelse(ivr$Site == "PDMO_Perharidy", "Perharidy", ivr$Site_bis)
170 ivr$Site_bis <- ifelse(ivr$Site == "BRES_Keraliou", "Keraliou", ivr$Site_bis)
171 ivr$Site_bis <- ifelse(ivr$Site == "FINS_Mousterlin", "Pointe de Mousterlin", ivr$Site_bis)
172 ivr$Site_bis <- ifelse(ivr$Site == "FINS_StNicolasGlenan", "Saint-Nicolas des Glénan", ivr$Site_bis)
173 ivr$Site_bis <- ifelse(ivr$Site == "FINS_AnseRoz", "Pointe de l'Anse du Roz", ivr$Site_bis)
174
175 ## change some variable format to integer
176 ivr$Nb.Blocs.Non.Retournes <- as.integer(ivr$Nb.Blocs.Non.Retournes)
177 ivr$Nb.Blocs.Retournes <- as.integer(ivr$Nb.Blocs.Retournes)
178
179 ivr$Year <- as.integer(ivr$Year)
180 ivr$Month <- as.integer(ivr$Month)
181 ivr$Day <- as.integer(ivr$Day)
182 ivr$Numero.Quadrat <- as.integer(ivr$Numero.Quadrat)
183
184
185 ## save the final, commplete ivr df.
186
187 ivr <- ivr[, c(19:54, 1:18)]
188
189
190 ## percentage of unturned vs overturned boulders and IVR previous 0-5 discrete scale values calculation
191
192 # create two new variables first
193 site_year_month_day <- paste0(ivr$Site, ".", gsub("-", ".", as.character(ivr$date_fiche)))
194 ivr <- tibble::add_column(ivr, site_year_month_day, .after = "Site_bis")
195 rm(site_year_month_day)
196
197 site_year_month_day_qdnb <- paste0(ivr$Site, ".", gsub("-", ".", as.character(ivr$Date)), ".", ivr$Numero.Quadrat)
198 ivr <- tibble::add_column(ivr, site_year_month_day_qdnb, .after = "site_year_month_day")
199 rm(site_year_month_day_qdnb)
200
201 ivr <- dplyr::arrange(ivr, Site, Year, Month, Numero.Quadrat)
202
203 # remove data with NA value for Nb.Blocs.Retournes & Nb.Blocs.Non.Retournes
204 ivr_naomit <- ivr %>% dplyr::filter(!is.na(ivr$Nb.Blocs.Retournes))
205 ivr_naomit <- as.data.frame(ivr_naomit)
206 colnames(ivr_naomit) <- colnames(ivr)
207 ivr_naomit <- ivr_naomit %>% dplyr::filter(!is.na(ivr_naomit$Nb.Blocs.Non.Retournes))
208 ivr_naomit <- as.data.frame(ivr_naomit)
209
210 # also remove data with Nb.Blocs.Retournes = 0 & Nb.Blocs.Non.Retournes = 0, cfr equivalent to quadrat with no boulders ... makes no sense to consider quadrat without boulder for ivr determination.
211 ivr_rm <- dplyr::filter(ivr_naomit, ivr_naomit$Nb.Blocs.Retournes == 0 && ivr_naomit$Nb.Blocs.Non.Retournes == 0)
212 ivr_naomit <- ivr_naomit %>% dplyr::anti_join(ivr_rm)
213 rm(ivr_rm)
214
215 ivr_val_qu_ <- ivr_naomit
216
217
218
219 #############################################################
220 # #
221 # Calcul of the IVR #
222 # #
223 #############################################################
224
225 ### Percentage of turned boulder
226 for (i in 1:nrow(ivr_naomit)) {
227 (bm <- sum(ivr_naomit$Nb.Blocs.Non.Retournes[i], ivr_naomit$Nb.Blocs.Retournes[i]))
228 (ivr_val_qu_$blocs.retournes.fr.[i] <- (ivr_naomit$Nb.Blocs.Retournes[i] / bm) * 100)
229 (ivr_val_qu_$blocs.non.retournes.fr.[i] <- (ivr_naomit$Nb.Blocs.Non.Retournes[i] / bm) * 100)
230 }
231
232 rm(bm, i)
233
234
235 ivr_val_qu_$blocs.non.retournes.fr. <- ifelse(is.nan(ivr_val_qu_$blocs.non.retournes.fr.), NA, ivr_val_qu_$blocs.non.retournes.fr.)
236 ivr_val_qu_$blocs.retournes.fr. <- ifelse(is.nan(ivr_val_qu_$blocs.retournes.fr.), NA, ivr_val_qu_$blocs.retournes.fr.)
237
238
239 # ivr for loop by quadrat.
240 for (i in 1:seq_len(nrow(ivr_val_qu_))) {
241 if (ivr_val_qu_$Nb.Blocs.Non.Retournes[i] == 0 && ivr_val_qu_$Nb.Blocs.Retournes[i] == 0) {
242 ivr_ <- NA
243 }else {
244 if (ivr_val_qu_$blocs.retournes.fr.[i] < 5) {
245 ivr_ <- 0
246 } else if (ivr_val_qu_$blocs.retournes.fr.[i] >= 5 && ivr_val_qu_$blocs.retournes.fr.[i] < 25) {
247 ivr_ <- 1
248 } else if (ivr_val_qu_$blocs.retournes.fr.[i] >= 25 && ivr_val_qu_$blocs.retournes.fr.[i] < 45) {
249 ivr_ <- 2
250 } else if (ivr_val_qu_$blocs.retournes.fr.[i] >= 45 && ivr_val_qu_$blocs.retournes.fr.[i] < 65) {
251 ivr_ <- 3
252 } else if (ivr_val_qu_$blocs.retournes.fr.[i] >= 65 && ivr_val_qu_$blocs.retournes.fr.[i] < 85) {
253 ivr_ <- 4
254 } else {
255 ivr_ <- 5
256 }
257
258 ivr_val_qu_$valeur.ivr_quadrat[i] <- ivr_
259 }
260 }
261
262 rm(i, ivr_)
263
264 # reorder variables for logical purpose
265 ivr_val_qu_ <- ivr_val_qu_[, c(1:56, 58, 57, 59)]
266 indic_full <- ivr_val_qu_
267 saveRDS(ivr_val_qu_, "ivr_val_qu.RDS")
268 rm(ivr_naomit)
269
270
271 ## Calculate ivr statistics now
272 ivr_val_qu_stat_ <- ivr_val_qu_ %>% dplyr::group_by(id.ivr, Site, Site_bis, Year, Month, Day) %>% dplyr::summarize(ivr_moy = mean(valeur.ivr_quadrat), ivr_et = sd(valeur.ivr_quadrat), ivr_med = median(valeur.ivr_quadrat), ivr_min = min(valeur.ivr_quadrat), ivr_max = max(valeur.ivr_quadrat), fr.r.moy = mean(blocs.retournes.fr.), fr.r.et = sd(blocs.retournes.fr.), fr.r.med = median(blocs.retournes.fr.), fr.r.min = min(blocs.retournes.fr.), fr.r.max = max(blocs.retournes.fr.), fr.nr.moy = mean(blocs.non.retournes.fr.), fr.nr.et = sd(blocs.non.retournes.fr.), fr.nr.med = median(blocs.non.retournes.fr.), fr.nr.min = min(blocs.non.retournes.fr.), fr.nr.max = max(blocs.non.retournes.fr.), nb. = dplyr::n())
273
274 Date <- as.Date(paste0(ivr_val_qu_stat_$Year, "-", ivr_val_qu_stat_$Month, "-", ivr_val_qu_stat_$Day), origin = "1970-01-01")
275 ivr_val_qu_stat_ <- tibble::add_column(ivr_val_qu_stat_, Date, .after = "Site_bis")
276 rm(Date)
277
278 ivr_val_qu_stat_ <- as.data.frame(ivr_val_qu_stat_)
279 indic <- ivr_val_qu_stat_
280
281
282 #############################################################
283 # #
284 # Plot the IVR per site #
285 # #
286 #############################################################
287
288 ## plot ivr (NB: Year, Month, Day variable names are replace by Annee, Mois, Jour, cfr previous label use in the script)
289 ivr_val_qu_stat_ <- dplyr::rename(ivr_val_qu_stat_, Annee = Year)
290 ivr_val_qu_stat_ <- dplyr::rename(ivr_val_qu_stat_, Mois = Month)
291 ivr_val_qu_stat_ <- dplyr::rename(ivr_val_qu_stat_, Jour = Day)
292
293
294 # new IVR scale with continuous 0 to 5 environmental status levels based on % of overturned boulders /20, plus other site data
295
296 for (i in c(1:length(unique(ivr_val_qu_stat_$Site)))) {
297
298 ivr_val_eg <- dplyr::filter(ivr_val_qu_stat_, ivr_val_qu_stat_$Site == unique(ivr_val_qu_stat_$Site)[i])
299
300 ivr_plot <- ggplot2::ggplot() +
301 ggplot2::geom_point(ggplot2::aes(x = ivr_val_qu_stat_$Date, y = ivr_val_qu_stat_$fr.r.moy / 20), col = "grey") +
302 ggplot2::geom_rect(ggplot2::aes(xmin = min(ivr_val_qu_stat_$Date), xmax = max(ivr_val_qu_stat_$Date), ymin = - 0.5, ymax = 5 / 20, fill = "#FF0000"), alpha = 0.3) +
303 ggplot2::geom_rect(ggplot2::aes(xmin = min(ivr_val_qu_stat_$Date), xmax = max(ivr_val_qu_stat_$Date), ymin = 5 / 20, ymax = 25 / 20, fill = "#F59404"), alpha = 0.3) +
304 ggplot2::geom_rect(ggplot2::aes(xmin = min(ivr_val_qu_stat_$Date), xmax = max(ivr_val_qu_stat_$Date), ymin = 25 / 20, ymax = 45 / 20, fill = "#FAFA15"), alpha = 0.3) +
305 ggplot2::geom_rect(ggplot2::aes(xmin = min(ivr_val_qu_stat_$Date), xmax = max(ivr_val_qu_stat_$Date), ymin = 45 / 20, ymax = 65 / 20, fill = "#18E125"), alpha = 0.3) +
306 ggplot2::geom_rect(ggplot2::aes(xmin = min(ivr_val_qu_stat_$Date), xmax = max(ivr_val_qu_stat_$Date), ymin = 65 / 20, ymax = 85 / 20, fill = "#04F5F5"), alpha = 0.3) +
307 ggplot2::geom_rect(ggplot2::aes(xmin = min(ivr_val_qu_stat_$Date), xmax = max(ivr_val_qu_stat_$Date), ymin = 85 / 20, ymax = 5.5, fill = "#1A1AE8"), alpha = 0.3) +
308 ggplot2::scale_fill_manual(values = c("#F59404", "#FAFA15", "#FF0000", "#04F5F5", "#18E125", "#1A1AE8")) +
309 ggplot2::geom_pointrange(ggplot2::aes(x = ivr_val_eg$Date, y = ivr_val_eg$fr.r.moy / 20, ymin = ivr_val_eg$fr.r.moy / 20 - ivr_val_eg$fr.r.et / 20, ymax = ivr_val_eg$fr.r.moy / 20 + ivr_val_eg$fr.r.et / 20), col = "black") +
310 ggplot2::xlab("Date") +
311 ggplot2::ylab("IVR") +
312 ggplot2::ggtitle(unique(ivr_val_eg$Site_bis)) +
313 ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1), legend.position = "none")
314
315 ggplot2::ggsave(paste0("ivr_", unique(ivr_val_eg$Site), ".png"), ivr_plot, height = 3, width = 3.5)
316
317
318 }
319
320 report <- args[3]
321 loop_file <- source(args[4])