diff --git a/alpha_beta.r b/alpha_beta.r new file mode 100644 index 0000000..c6446ea --- /dev/null +++ b/alpha_beta.r @@ -0,0 +1,141 @@ +#Rscript + +########################################### +## Mapping alpha and beta diversity ## +########################################### + +#####Packages : stars +# utils +# biodivmapr +# raster +# sf +# mapview +# leafpop +# RColorBrewer +# labdsv +# rgdal +# ggplot2 +# gridExtra + +#####Load arguments + +args <- commandArgs(trailingOnly = TRUE) + +#####Import the S2 data + +if (length(args) < 1) { + stop("This tool needs at least 1 argument") +}else { + data_raster <- args[1] + rasterheader <- args[2] + data <- args[3] + # type of PCA: + # PCA: no rescaling of the data + # SPCA: rescaling of the data + typepca <- as.character(args[4]) + alpha <- as.logical(args[5]) + beta <- as.logical(args[6]) + funct <- as.logical(args[7]) + all <- as.logical(args[8]) + source(args[9]) +} + +################################################################################ +## DEFINE PARAMETERS FOR DATASET TO BE PROCESSED ## +################################################################################ +if (data_raster == "") { + #Create a directory where to unzip your folder of data + dir.create("data_dir") + unzip(data, exdir = "data_dir") + # Path to raster + data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl") + input_image_file <- file.path("data_dir/results/Reflectance", data_raster[1]) + input_header_file <- file.path("data_dir/results/Reflectance", data_raster[2]) + +} else { + input_image_file <- file.path(getwd(), data_raster, fsep = "/") + input_header_file <- file.path(getwd(), rasterheader, fsep = "/") +} + +################################################################################ +## PROCESS IMAGE ## +################################################################################ +# 1- Filter data in order to discard non vegetated / shaded / cloudy pixels + +print("PERFORM PCA ON RASTER") +pca_output <- biodivMapR::perform_PCA(Input_Image_File = input_image_file, Input_Mask_File = input_mask_file, + Output_Dir = output_dir, TypePCA = typepca, FilterPCA = filterpca, nbCPU = nbcpu, MaxRAM = maxram) + +pca_files <- pca_output$PCA_Files +pix_per_partition <- pca_output$Pix_Per_Partition +nb_partitions <- pca_output$nb_partitions +# path for the updated mask +input_mask_file <- pca_output$MaskPath + + +selected_pcs <- seq(1, dim(raster::stack(input_image_file))[3]) + +selected_pcs <- all(selected_pcs) +################################################################################ +## MAP ALPHA AND BETA DIVERSITY ## +################################################################################ +print("MAP SPECTRAL SPECIES") + +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, SelectedPCs = selected_pcs, Pix_Per_Partition = pix_per_partition, nb_partitions = nb_partitions, TypePCA = typepca, nbCPU = nbcpu, MaxRAM = maxram, nbclusters = nbclusters) + +image_name <- tools::file_path_sans_ext(basename(input_image_file)) +if (alpha == TRUE || beta == TRUE || all == TRUE) { +## alpha + print("MAP ALPHA DIVERSITY") + index_alpha <- c("Shannon") + alpha_div <- biodivMapR::map_alpha_div(Input_Image_File = input_image_file, Output_Dir = output_dir, TypePCA = typepca, window_size = window_size, nbCPU = nbcpu, MaxRAM = maxram, Index_Alpha = index_alpha, nbclusters = nbclusters, FullRes = TRUE, LowRes = FALSE, MapSTD = FALSE) + + alpha_zip <- file.path(output_dir, image_name, typepca, "ALPHA", "Shannon_10_Fullres.zip") + alpha_path <- file.path(output_dir, image_name, typepca, "ALPHA") + unzip(alpha_zip, exdir = alpha_path) + alpha_path <- file.path(output_dir, image_name, typepca, "ALPHA", "Shannon_10_Fullres") + alpha_raster <- raster::raster(alpha_path) + get_alpha <- convert_raster(alpha_raster) + + if (alpha == TRUE || all == TRUE) { + colnames(get_alpha) <- c("Alpha", "longitude", "latitude") + plot_indices(get_alpha, titre = "Alpha") + + write.table(get_alpha, file = "alpha.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE) +} + if (beta == TRUE || all == TRUE) { +## beta + print("MAP BETA DIVERSITY") + beta_div <- biodivMapR::map_beta_div(Input_Image_File = input_image_file, Output_Dir = output_dir, TypePCA = typepca, window_size = window_size, nb_partitions = nb_partitions, nbCPU = nbcpu, MaxRAM = maxram, nbclusters = nbclusters) + + beta_path <- file.path(output_dir, image_name, typepca, "BETA", "BetaDiversity_BCdiss_PCO_10") + beta_raster <- raster::raster(beta_path) + get_beta <- convert_raster(beta_raster) + + colnames(get_beta) <- c("Beta", "longitude", "latitude") + plot_indices(get_beta, titre = "Beta") + + write.table(get_beta, file = "beta.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE) + } +} + + +################################################################################ +## COMPUTE ALPHA AND BETA DIVERSITY FROM FIELD PLOTS ## +################################################################################ + +if (funct == TRUE || all == TRUE) { + mapper <- biodivMapR::map_functional_div(Original_Image_File = input_image_file, Functional_File = pca_files, Selected_Features = selected_pcs, Output_Dir = output_dir, window_size = window_size, nbCPU = nbcpu, MaxRAM = maxram, TypePCA = typepca, FullRes = TRUE, LowRes = FALSE, MapSTD = FALSE) + + funct_zip <- file.path(output_dir, image_name, typepca, "FUNCTIONAL", "FunctionalDiversity_Map_MeanFilter_Fullres.zip") + funct_path <- file.path(output_dir, image_name, typepca, "FUNCTIONAL") + unzip(funct_zip, exdir = funct_path) + funct_path <- file.path(output_dir, image_name, typepca, "FUNCTIONAL", "FunctionalDiversity_Map_MeanFilter_Fullres") + funct_raster <- raster::raster(funct_path) + get_funct <- convert_raster(funct_raster) + + colnames(get_funct) <- c("Functionnal", "longitude", "latitude") + plot_indices(get_funct, titre = "Functionnal") + + write.table(get_funct, file = "Functionnal.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE) +} diff --git a/alpha_beta.xml b/alpha_beta.xml new file mode 100644 index 0000000..2c02c19 --- /dev/null +++ b/alpha_beta.xml @@ -0,0 +1,161 @@ + + from remote sensing data + + macrobis.xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + type == 'alpha' or type == 'all' + + + type == 'beta' or type == 'all' + + + type == 'funct' or type == 'all' + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/cb_freq.r b/cb_freq.r index 8f7f6e0..e223080 100644 --- a/cb_freq.r +++ b/cb_freq.r @@ -26,351 +26,315 @@ if (length(args) < 1) { input_soussite <- args[3] fiche_val <- args[4] input_obs <- args[5] - input_qecb <- args[6] + fiche_2 <- args[6] + input_qecb <- args[7] } -freq.hab. <- read.csv2(input_hab, fileEncoding = "Latin1") -freq.hab.$status <- "valide" +freq_hab <- read.csv2(input_hab, fileEncoding = "Latin1") +freq_hab$status <- "valide" -freq.site <- read.csv2(input_site, fileEncoding = "Latin1") -freq.site$status <- "valide" +freq_site <- read.csv2(input_site, fileEncoding = "Latin1") +freq_site$status <- "valide" -freq.ss.site <- read.csv2(input_soussite, fileEncoding = "Latin1") -freq.ss.site$status <- "valide" +freq_ss_site <- read.csv2(input_soussite, fileEncoding = "Latin1") +freq_ss_site$status <- "valide" ficheterrain <- read.csv2(fiche_val, fileEncoding = "Latin1") ficheterrain$date.sortie <- as.Date(ficheterrain$date.sortie) # only keep frequentation data related to boulder field in freqh.hab. df -#unique(freq.hab.$Libellé.Habitat) -freq.hab. <- dplyr::filter(freq.hab., grepl(c("champs | Champs | blocs"), Libellé.Habitat)) -#unique(freq.hab.$Libellé.Habitat) +freq_hab <- dplyr::filter(freq_hab, grepl(c("champs | Champs | blocs"), Libellé.Habitat)) -# make coinciding site names between frequentation df. and qecb df. -qecbNew <- readRDS(input_qecb) +# make coinciding site names between frequentation df and qecb df -#intersect(sort(unique(qecbNew$code.site)), sort(unique(freq.site$Code.Site))) -#setdiff(sort(unique(qecbNew$code.site)), intersect(sort(unique(qecbNew$code.site)), sort(unique(freq.site$Code.Site)))) # two sites are missing in freq.site compared to qecb data: "PNMI_02" = Quéménès & "PNMI_07" = Île de Sein - Goulenez et Île de Sein - Kilaourou. According to Anna Capietto (email Mercredi 19 Mai 2021 15:43:53): "Nous (...) n’avons pas de données de fréquentation associées à ces suivis." (i.e. qecb & ivr). -#unique(dplyr::filter(qecbNew, code.site %in% setdiff(sort(unique(qecbNew$code.site)), intersect(sort(unique(qecbNew$code.site)), sort(unique(freq.site$Code.Site)))))[c("Site", "code.site")]) -freq.site <- dplyr::filter(freq.site, Code.Site %in% c(sort(unique(qecbNew$code.site)))) +qecbnew <- readRDS(input_qecb) -#intersect(sort(unique(qecbNew$code.site)), sort(unique(freq.ss.site$Code.Site))) -#intersect(sort(unique(qecbNew$code.site)), sort(unique(freq.ss.site$Code.Sous.Site))) -freq.ss.site <- dplyr::filter(freq.ss.site, Code.Site %in% c(sort(unique(qecbNew$code.site)))) +freq_site <- dplyr::filter(freq_site, Code.Site %in% c(sort(unique(qecbnew$code.site)))) + +freq_ss_site <- dplyr::filter(freq_ss_site, Code.Site %in% c(sort(unique(qecbnew$code.site)))) `%notin%` <- Negate(`%in%`) -freq.ss.site <- dplyr::filter(freq.ss.site, Code.Sous.Site %notin% c("ARMO_082", "ARMO_153", "EGMP_114_2", "EGMP_015")) +freq_ss_site <- dplyr::filter(freq_ss_site, Code.Sous.Site %notin% c("ARMO_082", "ARMO_153", "EGMP_114_2", "EGMP_015")) -# some other possibilities to merge df. based on a common variable? +# some other possibilities to merge df based on a common variable? # obviously Code.Site is the only variable with common observations/values -#sort(unique(qecbNew$code.site)) -#sort(unique(ficheterrain$code.site)) -#intersect(sort(unique(qecbNew$code.site)), sort(unique(ficheterrain$code.site))) -#sort(unique(freq.site$Code.Site)) # same sites -ficheterrain <- dplyr::filter(ficheterrain, code.site %in% c(sort(unique(qecbNew$code.site)))) +ficheterrain <- dplyr::filter(ficheterrain, code.site %in% c(sort(unique(qecbnew$code.site)))) # intersect between ficheterrain and data according to ID.Fiche variable, cfr Elodie Gamp discu. ficheterrain <- dplyr::rename(ficheterrain, Libellé.Sortie = libellé.sortie) -freq.site <- dplyr::left_join(freq.site, ficheterrain, by = c("ID.Fiche" +freq_site <- dplyr::left_join(freq_site, ficheterrain, by = c("ID.Fiche" , "Libellé.Sortie" )) -freq.ss.site <- dplyr::left_join(freq.ss.site, ficheterrain, by = c("ID.Fiche" +freq_ss_site <- dplyr::left_join(freq_ss_site, ficheterrain, by = c("ID.Fiche" , "Libellé.Sortie" )) -freq.hab. <- dplyr::left_join(freq.hab., ficheterrain, by = c("ID.Fiche" +freq_hab <- dplyr::left_join(freq_hab, ficheterrain, by = c("ID.Fiche" , "Libellé.Sortie" )) -# I noticed an issue with -999 data instead of NAs for freq.hab. df. I correct it here below. -min(freq.hab.$Nb.Total) -freq.hab.$Nb.Total <- ifelse(freq.hab.$Nb.Total == min(freq.hab.$Nb.Total), NA, freq.hab.$Nb.Total) +# I noticed an issue with -999 data instead of NAs for freq_hab df I correct it here below. +min(freq_hab$Nb.Total) +freq_hab$Nb.Total <- ifelse(freq_hab$Nb.Total == min(freq_hab$Nb.Total), NA, freq_hab$Nb.Total) -unique(freq.hab.[, c("code.site", "Code.Habitat", "Libellé.Habitat")]) +unique(freq_hab[, c("code.site", "Code.Habitat", "Libellé.Habitat")]) # only Flots Bleus is repeated: 1 site, 2 ss-sites -# insert site name based on qecb df. +# insert site name based on qecb df -(df. <- unique(qecbNew[, c("code.site", "Site", "Site_bis")])) +(df <- unique(qecbnew[, c("code.site", "Site", "Site_bis")])) ########################################################### -# Show df.join table throughout the process to check for site name correspondence. +# Show df_join table throughout the process to check for site name correspondence. #!!!!!!!!!!!!!!!!!!!!!!!!!!!! ########################################################### -(df.join <- unique(freq.hab.[, c("code.site", "site", "Code.Habitat", "Libellé.Habitat")])) -(df.join <- dplyr::left_join(df.join, df., by = "code.site")) +(df_join <- unique(freq_hab[, c("code.site", "site", "Code.Habitat", "Libellé.Habitat")])) +(df_join <- dplyr::left_join(df_join, df, by = "code.site")) ########################################################### #!!!!!!!!!!!!!!!!!!!!!!!!!!!! # check below kine after running code because lines to be removed might need to be updated !! or code another way with character strings -(df.join <- df.join[-c(9,10),]) # current +(df_join <- df_join[-c(9,10),]) # current #!!!!!!!!!!!!!!!!!!!!!!!!!!!! ########################################################### -freq.hab.$Site <- NA -freq.hab.$Site <- as.character(freq.hab.$Site) +freq_hab$Site <- NA +freq_hab$Site <- as.character(freq_hab$Site) library(data.table) # multiple ifelse statement doesn't work, so I used another function from package data.table -for (i in c(1:nrow(df.join))) { +for (i in c(1:nrow(df_join))) { #i <- 1 - setDT(freq.hab.)[code.site == df.join[i, "code.site"] & Code.Habitat == df.join[i, "Code.Habitat"], Site := df.join[i,"Site"]] + setDT(freq_hab)[code.site == df_join[i, "code.site"] & Code.Habitat == df_join[i, "Code.Habitat"], Site := df_join[i, "Site"]] } -freq.hab. <- data.frame(freq.hab.) +freq_hab <- data.frame(freq_hab) # NB: from table to dataframe space " " are replace by a ".", therefore colnames changes, e.g. below for "code.site" and "Code.Habitat". -#unique(freq.hab.[, c("code.site", "Code.Habitat", "Site")]) -freq.hab. <- tibble::add_column(freq.hab.[, c(1:(ncol(freq.hab.)-1))], Site = freq.hab.$Site, .after = "Libellé.Habitat") -freq.hab.$Site_bis <- NA -freq.hab.$Site_bis <- as.character(freq.hab.$Site_bis) + +freq_hab <- tibble::add_column(freq_hab[, c(1:(ncol(freq_hab)-1))], Site = freq_hab$Site, .after = "Libellé.Habitat") +freq_hab$Site_bis <- NA +freq_hab$Site_bis <- as.character(freq_hab$Site_bis) library(data.table) # multiple ifelse statement doesn't work, so I used another function from package data.table -for (i in c(1:nrow(df.join))) { +for (i in c(1:nrow(df_join))) { #i <- 1 - setDT(freq.hab.)[code.site == df.join[i, "code.site"] & Code.Habitat == df.join[i, "Code.Habitat"], Site_bis := df.join[i,"Site_bis"]] + setDT(freq_hab)[code.site == df_join[i, "code.site"] & Code.Habitat == df_join[i, "Code.Habitat"], Site_bis := df_join[i, "Site_bis"]] } -freq.hab. <- data.frame(freq.hab.) -#unique(freq.hab.[, c("code.site", "Code.Habitat", "Site", "Site_bis")]) -freq.hab. <- tibble::add_column(freq.hab.[, c(1:(ncol(freq.hab.)-1))], Site_bis = freq.hab.$Site_bis, .after = "Site") +freq_hab <- data.frame(freq_hab) -rm(df.join) +freq_hab <- tibble::add_column(freq_hab[, c(1:(ncol(freq_hab)-1))], Site_bis = freq_hab$Site_bis, .after = "Site") -#unique(freq.site[, c("sous.site", "zone.habitat")]) -(df.join <- unique(freq.site[, c("Code.Site", "Libellé.Site" +rm(df_join) + +(df_join <- unique(freq_site[, c("Code.Site", "Libellé.Site" #, "site" )])) -df.join <- dplyr::rename(df.join, code.site = Code.Site) -(df.join <- dplyr::left_join(df.join, df., by = "code.site")) -df.join$Site <- ifelse(df.join$code.site == "BASQ_01", "BASQ_FlotsBleus", df.join$Site) -df.join$Site <- ifelse(df.join$code.site == "ARMO_042-043", "ARMO_Piégu / Verdelet", df.join$Site) -df.join$Site_bis <- ifelse(df.join$code.site == "BASQ_01", "Les Flots Bleus", df.join$Site_bis) -df.join$Site_bis <- ifelse(df.join$code.site == "ARMO_042-043", "Îlot du Verdelet / Piégu", df.join$Site_bis) -df.join <- df.join[!duplicated(df.join), ] -freq.site$Site <- NA -freq.site$Site <- as.character(freq.site$Site) +df_join <- dplyr::rename(df_join, code.site = Code.Site) +(df_join <- dplyr::left_join(df_join, df, by = "code.site")) +df_join$Site <- ifelse(df_join$code.site == "BASQ_01", "BASQ_FlotsBleus", df_join$Site) +df_join$Site <- ifelse(df_join$code.site == "ARMO_042-043", "ARMO_Piégu / Verdelet", df_join$Site) +df_join$Site_bis <- ifelse(df_join$code.site == "BASQ_01", "Les Flots Bleus", df_join$Site_bis) +df_join$Site_bis <- ifelse(df_join$code.site == "ARMO_042-043", "Îlot du Verdelet / Piégu", df_join$Site_bis) +df_join <- df_join[!duplicated(df_join), ] +freq_site$Site <- NA +freq_site$Site <- as.character(freq_site$Site) library(data.table) -for (i in c(1:nrow(df.join))) { +for (i in c(1:nrow(df_join))) { #i <- 1 - setDT(freq.site)[Code.Site == df.join[i, "code.site"], Site := df.join[i,"Site"]] + setDT(freq_site)[Code.Site == df_join[i, "code.site"], Site := df_join[i,"Site"]] } -freq.site <- data.frame(freq.site) -#unique(freq.site[, c("Code.Site", "Libellé.Site", "Site")]) -freq.site <- tibble::add_column(freq.site[, c(1:(ncol(freq.site)-1))], Site = freq.site$Site, .after = "Libellé.Site") +freq_site <- data.frame(freq_site) + +freq_site <- tibble::add_column(freq_site[, c(1:(ncol(freq_site)-1))], Site = freq_site$Site, .after = "Libellé.Site") # check for Flots Bleus and Piégu / Verdelet ? -freq.site$Site_bis <- NA -freq.site$Site_bis <- as.character(freq.site$Site_bis) +freq_site$Site_bis <- NA +freq_site$Site_bis <- as.character(freq_site$Site_bis) library(data.table) -for (i in c(1:nrow(df.join))) { +for (i in c(1:nrow(df_join))) { #i <- 1 - setDT(freq.site)[Code.Site == df.join[i, "code.site"], Site_bis := df.join[i,"Site_bis"]] + setDT(freq_site)[Code.Site == df_join[i, "code.site"], Site_bis := df_join[i,"Site_bis"]] } -freq.site <- data.frame(freq.site) +freq_site <- data.frame(freq_site) -freq.site <- tibble::add_column(freq.site[, c(1:(ncol(freq.site)-1))], Site_bis = freq.site$Site_bis, .after = "Site") +freq_site <- tibble::add_column(freq_site[, c(1:(ncol(freq_site)-1))], Site_bis = freq_site$Site_bis, .after = "Site") # check for Flots Bleus and Piégu / Verdelet ? -#unique(dplyr::filter(freq.site, Code.Site == "BASQ_01")[, c("Code.Site", "Site_bis")]) -#unique(dplyr::filter(freq.site, Code.Site == "ARMO_042-043")[, c("Code.Site", "Site_bis")]) -rm(df.join) - -#unique(freq.ss.site[, c("sous.site", "zone.habitat")]) -(df.join <- unique(freq.ss.site[, c("Code.Site", "Libellé.Site", "Code.Sous.Site", "Libellé.Sous.Site")])) -df.join <- dplyr::rename(df.join, code.site = Code.Site) -(df.join <- dplyr::left_join(df.join, df., by = "code.site")) +rm(df_join) + + +(df_join <- unique(freq_ss_site[, c("Code.Site", "Libellé.Site", "Code.Sous.Site", "Libellé.Sous.Site")])) +df_join <- dplyr::rename(df_join, code.site = Code.Site) +(df_join <- dplyr::left_join(df_join, df, by = "code.site")) ########################################################### #!!!!!!!!!!!!!!!!!!!!!!!!!!!! # check below kine after running code because lines to be removed might need to be updated !! or code another way with character strings -(df.join <- df.join[-c(1,5),]) # current +(df_join <- df_join[-c(1,5),]) # current ########################################################### -freq.ss.site$Site <- NA -freq.ss.site$Site <- as.character(freq.ss.site$Site) +freq_ss_site$Site <- NA +freq_ss_site$Site <- as.character(freq_ss_site$Site) library(data.table) -for (i in c(1:nrow(df.join))) { +for (i in c(1:nrow(df_join))) { #i <- 1 - setDT(freq.ss.site)[Code.Site == df.join[i, "code.site"] & Code.Sous.Site == df.join[i, "Code.Sous.Site"], Site := df.join[i,"Site"]] + setDT(freq_ss_site)[Code.Site == df_join[i, "code.site"] & Code.Sous.Site == df_join[i, "Code.Sous.Site"], Site := df_join[i,"Site"]] } -freq.ss.site <- data.frame(freq.ss.site) -#unique(freq.ss.site[, c("Code.Site", "Libellé.Site", "Code.Sous.Site", "Libellé.Sous.Site", "Site")]) -freq.ss.site <- tibble::add_column(freq.ss.site[, c(1:(ncol(freq.ss.site)-1))], Site = freq.ss.site$Site, .after = "Libellé.Site") +freq_ss_site <- data.frame(freq_ss_site) + +freq_ss_site <- tibble::add_column(freq_ss_site[, c(1:(ncol(freq_ss_site)-1))], Site = freq_ss_site$Site, .after = "Libellé.Site") # check for Piégu / Verdelet ? -#unique(dplyr::filter(freq.ss.site, Code.Site == "ARMO_042-043")[, c("Code.Site", "Site")]) -freq.ss.site$Site_bis <- NA -freq.ss.site$Site_bis <- as.character(freq.ss.site$Site_bis) + +freq_ss_site$Site_bis <- NA +freq_ss_site$Site_bis <- as.character(freq_ss_site$Site_bis) library(data.table) -for (i in c(1:nrow(df.join))) { +for (i in c(1:nrow(df_join))) { #i <- 1 - setDT(freq.ss.site)[Code.Site == df.join[i, "code.site"] & Code.Sous.Site == df.join[i, "Code.Sous.Site"], Site_bis := df.join[i,"Site_bis"]] + setDT(freq_ss_site)[Code.Site == df_join[i, "code.site"] & Code.Sous.Site == df_join[i, "Code.Sous.Site"], Site_bis := df_join[i,"Site_bis"]] } -freq.ss.site <- data.frame(freq.ss.site) -#unique(freq.ss.site[, c("Code.Site", "Libellé.Site", "Code.Sous.Site", "Libellé.Sous.Site", "Site", "Site_bis")]) -freq.ss.site <- tibble::add_column(freq.ss.site[, c(1:(ncol(freq.ss.site)-1))], Site_bis = freq.ss.site$Site_bis, .after = "Site") +freq_ss_site <- data.frame(freq_ss_site) + +freq_ss_site <- tibble::add_column(freq_ss_site[, c(1:(ncol(freq_ss_site)-1))], Site_bis = freq_ss_site$Site_bis, .after = "Site") # check for Piégu / Verdelet ? -#unique(dplyr::filter(freq.ss.site, Code.Site == "ARMO_042-043")[, c("Code.Site", "Site", "Site_bis")]) -rm(df.join) -rm(df.) +rm(df_join) + +rm(df) # Conclusions from data mining (not shown) regarding frequentation data ; considerations required for better understanding of the script. not needed I guess in an application. -# for freq.hab., consider $Nb.Total var. -# for freq.site, consider $Nb.Total and $Nb.Pecheurs.Site var. (also eventually add $`Nb Pecheurs Arrivé` to $Nb.Pecheurs.Site to have a more exhaustif nb.); NB: $Nb.Pecheurs.Site is major part of $Nb.Total, but not always specified as such, so better to only consider $Nb.Total -# Ccl.: all count data in freq.ss.site are within freq.site, and freq.hab. data are little, and no correlation between both site and habitat count ! +# for freq_hab, consider $Nb.Total var. +# for freq_site, consider $Nb.Total and $Nb.Pecheurs.Site var. (also eventually add $`Nb Pecheurs Arrivé` to $Nb.Pecheurs.Site to have a more exhaustif nb.); NB: $Nb.Pecheurs.Site is major part of $Nb.Total, but not always specified as such, so better to only consider $Nb.Total +# Ccl.: all count data in freq_ss_site are within freq_site, and freq_hab data are little, and no correlation between both site and habitat count ! # see now if df are complementary with regards to frequentation # NB: Now we consider "Site" according to geographic level, i.e. Site vs Ss.Site vs Habitat (cfr e.g. FlotsBleus, with 2 Ss.Sites). -freq.hab.red. <- freq.hab.[, c("Code.Habitat", "Libellé.Habitat", "code.site", "Site", "Site_bis", "site", "date.sortie", "Nb.Total", "Nb.Adultes", "Nb.Enfants")] -#freq.hab.red. <- dplyr::rename(freq.hab.red., Nb.Total.freq.hab. = Nb.Total) -colnames(freq.hab.red.)[c(3,8:10)] <- c(paste0("FrHa_", names(freq.hab.red.[, c(3,8:10)]))) -freq.hab.red.$Site.date.sortie <- paste0(freq.hab.red.$Site, ".", freq.hab.red.$date.sortie) +freq_habred_ <- freq_hab[, c("Code.Habitat", "Libellé.Habitat", "code.site", "Site", "Site_bis", "site", "date.sortie", "Nb.Total", "Nb.Adultes", "Nb.Enfants")] + +colnames(freq_habred_)[c(3, 8:10)] <- c(paste0("frha_", names(freq_habred_[, c(3, 8:10)]))) +freq_habred_$Site.date.sortie <- paste0(freq_habred_$Site, ".", freq_habred_$date.sortie) # add replicate number when several counts the same date -freq.hab.red. <- freq.hab.red. %>% +freq_habred_ <- freq_habred_ %>% dplyr::group_by(Site.date.sortie) %>% dplyr::mutate(repl. = dplyr::row_number(Site.date.sortie)) -freq.hab.red. <- data.frame(freq.hab.red.) +freq_habred_ <- data.frame(freq_habred_) -freq.ss.site.red. <- freq.ss.site[, c("Code.Sous.Site", "Libellé.Sous.Site", "Code.Site", "Libellé.Site", "Site", "Site_bis", "site", "date.sortie", "Nb.Total", "Nb.Adultes", "Nb.Enfants")] -colnames(freq.ss.site.red.)[c(3,4,9:11)] <- c(paste0("FrSsSi_", names(freq.ss.site.red.[, c(3,4,9:11)]))) -freq.ss.site.red.$Site.date.sortie <- paste0(freq.ss.site.red.$Site, ".", freq.ss.site.red.$date.sortie) -freq.ss.site.red. <- dplyr::arrange(freq.ss.site.red., Site, date.sortie, FrSsSi_Nb.Total) -freq.ss.site.red. <- freq.ss.site.red. %>% +freq_ss_site_red_ <- freq_ss_site[, c("Code.Sous.Site", "Libellé.Sous.Site", "Code.Site", "Libellé.Site", "Site", "Site_bis", "site", "date.sortie", "Nb.Total", "Nb.Adultes", "Nb.Enfants")] +colnames(freq_ss_site_red_)[c(3, 4, 9:11)] <- c(paste0("frsssi_", names(freq_ss_site_red_[, c(3, 4, 9:11)]))) +freq_ss_site_red_$Site.date.sortie <- paste0(freq_ss_site_red_$Site, ".", freq_ss_site_red_$date.sortie) +freq_ss_site_red_ <- dplyr::arrange(freq_ss_site_red_, Site, date.sortie, frsssi_Nb.Total) +freq_ss_site_red_ <- freq_ss_site_red_ %>% dplyr::group_by(Site.date.sortie) %>% dplyr::mutate(repl. = dplyr::row_number(Site.date.sortie)) -freq.ss.site.red. <- data.frame(freq.ss.site.red.) +freq_ss_site_red_ <- data.frame(freq_ss_site_red_) -freq.site.red. <- freq.site[, c("Code.Site", "Libellé.Site", "Site", "Site_bis", "site", "date.sortie", "Nb.Total", "Nb.Adultes", "Nb.Enfants", "Nb.Pecheurs.Site", "Nb.Pecheurs.Arrivee", "Nb.Pecheurs.Départ", "Nb.Pecheurs.Zone.Interdite")] -colnames(freq.site.red.)[c(1,2,7:13)] <- c(paste0("FrSi_", names(freq.site.red.[, c(1,2,7:13)]))) -freq.site.red.$Site.date.sortie <- paste0(freq.site.red.$Site, ".", freq.site.red.$date.sortie) -freq.site.red. <- dplyr::arrange(freq.site.red., Site, date.sortie, FrSi_Nb.Total) -freq.site.red. <- freq.site.red. %>% +freq_site_red_ <- freq_site[, c("Code.Site", "Libellé.Site", "Site", "Site_bis", "site", "date.sortie", "Nb.Total", "Nb.Adultes", "Nb.Enfants", "Nb.Pecheurs.Site", "Nb.Pecheurs.Arrivee", "Nb.Pecheurs.Départ", "Nb.Pecheurs.Zone.Interdite")] +colnames(freq_site_red_)[c(1, 2, 7:13)] <- c(paste0("frsi_", names(freq_site_red_[, c(1, 2, 7:13)]))) +freq_site_red_$Site.date.sortie <- paste0(freq_site_red_$Site, ".", freq_site_red_$date.sortie) +freq_site_red_ <- dplyr::arrange(freq_site_red_, Site, date.sortie, frsi_Nb.Total) +freq_site_red_ <- freq_site_red_ %>% dplyr::group_by(Site.date.sortie) %>% dplyr::mutate(repl. = dplyr::row_number(Site.date.sortie)) -freq.site.red. <- data.frame(freq.site.red.) +freq_site_red_ <- data.frame(freq_site_red_) -df. <- dplyr::full_join(freq.site.red.#[, c("Site", "date.sortie", "Site.date.sortie", "repl.", "Nb.Total.freq.site")] - , freq.ss.site.red.#[, c("Site", "date.sortie", "Site.date.sortie", "repl.", "Nb.Total.freq.ss.site")] +df <- dplyr::full_join(freq_site_red_#[, c("Site", "date.sortie", "Site.date.sortie", "repl.", "Nb.Total.freq_site")] + , freq_ss_site_red_#[, c("Site", "date.sortie", "Site.date.sortie", "repl.", "Nb.Total.freq_ss_site")] , by = c("Site", "Site_bis", "site", "date.sortie", "Site.date.sortie", "repl.")) -df. <- dplyr::full_join(df., freq.hab.red.#[, c("Site", "date.sortie", "Site.date.sortie", "repl.", "Nb.Total.freq.hab.")] +df <- dplyr::full_join(df, freq_habred_#[, c("Site", "date.sortie", "Site.date.sortie", "repl.", "Nb.Total.freq_hab")] , by = c("Site", "Site_bis", "site", "date.sortie", "Site.date.sortie", "repl.")) -names(df.) -df. <- dplyr::arrange(df., Site, date.sortie, repl.) -freq. <- data.frame(df.) -freq. <- freq.[, c(3:6, 14, 15, 1, 2, 7:13, 16:28)] -rm(df.) +df <- dplyr::arrange(df, Site, date.sortie, repl.) +freq_ <- data.frame(df) +freq_ <- freq_[, c(3:6, 14, 15, 1, 2, 7:13, 16:28)] +rm(df) -freq. <- tidyr::separate(freq., date.sortie, into = c("Annee", "Mois", "Jour"), remove = FALSE) -freq.$Annee <- as.numeric(freq.$Annee) -freq.$Mois <- as.numeric(freq.$Mois) -freq.$Jour <- as.numeric(freq.$Jour) +freq_ <- tidyr::separate(freq_, date.sortie, into = c("Annee", "Mois", "Jour"), remove = FALSE) +freq_$Annee <- as.numeric(freq_$Annee) +freq_$Mois <- as.numeric(freq_$Mois) +freq_$Jour <- as.numeric(freq_$Jour) # prior saving, I will create 3 new variables in the freq dataset: one with freq data, one with the corresponding origin variable, one for the diff between frequ data -freq.$Fr_Nb <- NA +freq_$Fr_Nb <- NA -for (i in 1:nrow(freq.)){ -if (!is.na(freq.$FrSi_Nb.Total[i])) { - freq.$Fr_Nb[i] = freq.$FrSi_Nb.Total[i] - } else if (is.na(freq.$FrSi_Nb.Total[i]) & !is.na(freq.$FrSi_Nb.Pecheurs.Site[i])) { - freq.$Fr_Nb[i] = freq.$FrSi_Nb.Pecheurs.Site[i] - } else if (is.na(freq.$FrSi_Nb.Total[i]) & !is.na(freq.$FrSsSi_Nb.Total[i])) { - freq.$Fr_Nb[i] = freq.$FrSsSi_Nb.Total[i] +for (i in 1:nrow(freq_)){ +if (!is.na(freq_$frsi_Nb.Total[i])) { + freq_$Fr_Nb[i] = freq_$frsi_Nb.Total[i] + } else if (is.na(freq_$frsi_Nb.Total[i]) & !is.na(freq_$frsi_Nb.Pecheurs.Site[i])) { + freq_$Fr_Nb[i] = freq_$frsi_Nb.Pecheurs.Site[i] + } else if (is.na(freq_$frsi_Nb.Total[i]) & !is.na(freq_$frsssi_Nb.Total[i])) { + freq_$Fr_Nb[i] = freq_$frsssi_Nb.Total[i] } else { - freq.$Fr_Nb[i] = freq.$FrHa_Nb.Total[i] + freq_$Fr_Nb[i] = freq_$frha_Nb.Total[i] } } -freq.$Fr_var <- NA +freq_$Fr_var <- NA -for (i in 1:nrow(freq.)){ - if (!is.na(freq.$FrSi_Nb.Total[i])) { - freq.$Fr_var[i] = "FrSi_Nb.Total" - } else if (is.na(freq.$FrSi_Nb.Total[i]) & !is.na(freq.$FrSi_Nb.Pecheurs.Site[i])) { - freq.$Fr_var[i] = "FrSi_Nb.Pecheurs.Site" - } else if (is.na(freq.$FrSi_Nb.Total[i]) & !is.na(freq.$FrSsSi_Nb.Total[i])) { - freq.$Fr_var[i] = "FrSsSi_Nb.Total" +for (i in 1:nrow(freq_)){ + if (!is.na(freq_$frsi_Nb.Total[i])) { + freq_$Fr_var[i] = "frsi_Nb.Total" + } else if (is.na(freq_$frsi_Nb.Total[i]) & !is.na(freq_$frsi_Nb.Pecheurs.Site[i])) { + freq_$Fr_var[i] = "frsi_Nb.Pecheurs.Site" + } else if (is.na(freq_$frsi_Nb.Total[i]) & !is.na(freq_$frsssi_Nb.Total[i])) { + freq_$Fr_var[i] = "frsssi_Nb.Total" } else { - freq.$Fr_var[i] = "FrHa_Nb.Total" + freq_$Fr_var[i] = "frha_Nb.Total" } } -freq.$Fr_diff <- NA +freq_$Fr_diff <- NA -for (i in 1:nrow(freq.)){ -freq.$Fr_diff[i] <- abs(max(freq.[i, c("FrSi_Nb.Total", "FrSi_Nb.Pecheurs.Site", "FrSsSi_Nb.Total", "FrHa_Nb.Total")], na.rm = T) - freq.$Fr_Nb[i]) +for (i in 1:nrow(freq_)){ +freq_$Fr_diff[i] <- abs(max(freq_[i, c("frsi_Nb.Total", "frsi_Nb.Pecheurs.Site", "frsssi_Nb.Total", "frha_Nb.Total")], na.rm = TRUE) - freq_$Fr_Nb[i]) } - -#saveRDS(freq., "freq.RDS") - - # final plot # almost only unique value ! -freq. %>% +freq_ %>% dplyr::group_by(repl.) %>% dplyr::summarise(nb = dplyr::n()) -Ymin. <- min(freq.$Annee) -Ymax. <- max(freq.$Annee) +ymin_ <- min(freq_$Annee) +ymax_ <- max(freq_$Annee) -for (i in c(1:length(unique(freq.[, "Site"])))) { +for (i in c(1:length(unique(freq_[, "Site"])))) { #i <- 3 - freq.i <- dplyr::filter(freq., Site == unique(freq.[, "Site"])[i]) + freq_i <- dplyr::filter(freq_, Site == unique(freq_[, "Site"])[i]) - xmin. <- as.Date(ifelse(min(freq.i$Annee, na.rm = T) >= 2014, as.Date("2014-01-01", origin = "1970-01-01"), as.Date(paste0(Ymin., "-01-01"), origin = "1970-01-01")), origin = "1970-01-01") - xmax. <- as.Date(ifelse(max(freq.i$Annee, na.rm = T) <= 2017, as.Date("2018-01-01", origin = "1970-01-01"), as.Date(paste0((Ymax.+1), "-01-01"), origin = "1970-01-01")), origin = "1970-01-01") + xmin_ <- as.Date(ifelse(min(freq_i$Annee, na.rm = TRUE) >= 2014, as.Date("2014-01-01", origin = "1970-01-01"), as.Date(paste0(ymin_, "-01-01"), origin = "1970-01-01")), origin = "1970-01-01") + xmax_ <- as.Date(ifelse(max(freq_i$Annee, na.rm = TRUE) <= 2017, as.Date("2018-01-01", origin = "1970-01-01"), as.Date(paste0((ymax_ + 1), "-01-01"), origin = "1970-01-01")), origin = "1970-01-01") # I have to invert plot(x,y) by plot(y,x) in order to work ... ? - png(paste0("freq_1_", unique(freq.i$Site), ".png")) - plot(freq.i$date.sortie, rep(0, length(freq.i$date.sortie)), - xlim = c(xmin., xmax.), - ylim = c(-8, 165) -#round(0-((max(#freq.i[, "FrSi_Nb.Total"], - #freq.i[, "FrSsSi_Nb.Total"]#, freq.i[, "FrHa_Nb.Total"] -#, na.rm = T)*1.05)-max(#freq.i[, "FrSi_Nb.Total"], -#freq.i[, "FrSsSi_Nb.Total"]#, freq.i[, "FrHa_Nb.Total"] -#, na.rm = T))), #log10 - # round(max(#freq.i[, "FrSi_Nb.Total"], -#freq.i[, "FrSsSi_Nb.Total"]#, freq.i[, "FrHa_Nb.Total"] -#, na.rm = T)*1.05)), - , main = unique(freq.i$Site), xlab = "", ylab = "fréquentation", col = "white") - - #points(c(freq.$FrSi_Nb.Total, freq.$FrSsSi_Nb.Total, freq.$FrHa_Nb.Total), c(freq.$date.sortie, freq.$date.sortie, freq.$date.sortie), pch = 19, col = "lightgrey", cex = .5) - - points(freq.i$date.sortie, #log10 - (freq.i$FrSi_Nb.Total - #+1 - ), pch = 19, col = "red") - points(freq.i$date.sortie, #log10 - (freq.i$FrSsSi_Nb.Total - #+1 - ), pch = 19, col = "darkblue") - points(freq.i$date.sortie, #log10 - (freq.i$FrHa_Nb.Total - #+1 - ), pch = 19, col = "forestgreen") - - legend("bottom", inset = c(0,-0.4), legend = c("Si","SsSi", "Ha"), pch = c(19,19,19), col = c("red", "darkblue", "forestgreen"), horiz = TRUE, bty = "n", xpd = TRUE) - #inset = c(-0.45, 0) # You will need to fine-tune the first value depending on the windows size - # xpd = TRUE # You need to specify this graphical parameter to put the legend outside the plot - + freq_plot <- ggplot2::ggplot() + + #ggplot2::geom_point(ggplot2::aes(x = freq_$date.sortie, y = freq_$frsi_Nb.Total), col = "grey") + + #ggplot2::geom_point(ggplot2::aes(x = freq_$date.sortie, y = freq_$frsssi_Nb.Total), col = "grey") + + #ggplot2::geom_point(ggplot2::aes(x = freq_$date.sortie, y = freq_$frha_Nb.Total), col = "grey") + + ggplot2::geom_point(ggplot2::aes(x = freq_i$date.sortie, y = freq_i$frsi_Nb.Total), col = "red") + + ggplot2::geom_point(ggplot2::aes(x = freq_i$date.sortie, y = freq_i$frsssi_Nb.Total), col = "darkblue") + + ggplot2::geom_point(ggplot2::aes(x = freq_i$date.sortie, y = freq_i$frha_Nb.Total), col = "forestgreen") + + ggplot2::xlab("Date") + + ggplot2::ylab("fréquentation") + + ggplot2::ggtitle(unique(freq_i$Site)) + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1), legend.position = "bottom") + +ggplot2::ggsave(paste0("freq_", unique(freq_i$Site), ".png"), freq_plot, height = 3, width = 3.5) + } @@ -378,338 +342,249 @@ for (i in c(1:length(unique(freq.[, "Site"])))) { # remove data if any (not the case anymore here) without date.sortie value. -freq. <- freq. %>% dplyr::filter(!is.na(freq.$date.sortie)) +freq_ <- freq_ %>% dplyr::filter(!is.na(freq_$date.sortie)) -table(( na.omit(freq.[, c("FrSi_Code.Site", "Site.date.sortie", "FrSi_Nb.Total")]) %>% - dplyr::group_by(FrSi_Code.Site, Site.date.sortie) %>% +table((na.omit(freq_[, c("frsi_Code.Site", "Site.date.sortie", "frsi_Nb.Total")]) %>% + dplyr::group_by(frsi_Code.Site, Site.date.sortie) %>% dplyr::summarise(nb = dplyr::n()) )["nb"]) -table(( na.omit(freq.[, c("Code.Sous.Site", "Site.date.sortie", "FrSsSi_Nb.Total")]) %>% +table((na.omit(freq_[, c("Code.Sous.Site", "Site.date.sortie", "frsssi_Nb.Total")]) %>% dplyr::group_by(Code.Sous.Site, Site.date.sortie) %>% dplyr::summarise(nb = dplyr::n()) )["nb"]) -table(( na.omit(freq.[, c("Code.Habitat", "Site.date.sortie", "FrHa_Nb.Total")]) %>% +table((na.omit(freq_[, c("Code.Habitat", "Site.date.sortie", "frha_Nb.Total")]) %>% dplyr::group_by(Code.Habitat, Site.date.sortie) %>% dplyr::summarise(nb = dplyr::n()) )["nb"]) # median and mean of two values are identical, so dplyr::filter for nb > 2 -# and for remaining data, I'll consider the mean value cfr df. above +# and for remaining data, I'll consider the mean value cfr df above # first add a semester variable -freq.$Mois <- as.numeric(freq.$Mois) -ifelse(freq.$Mois %in% c(1:6), "s1", "s2") -> freq.$semester -freq.$semester <- as.factor(freq.$semester) -freq.$Annee_semester <- paste0(freq.$Annee, "_", freq.$semester) +freq_$Mois <- as.numeric(freq_$Mois) +freq_$semester <- ifelse(freq_$Mois %in% c(1:6), "s1", "s2") +freq_$semester <- as.factor(freq_$semester) +freq_$Annee_semester <- paste0(freq_$Annee, "_", freq_$semester) -FrSi <- na.omit(freq.[, c("Site", "Site_bis", "FrSi_Code.Site", "Site.date.sortie", "Annee_semester", "FrSi_Nb.Total")]) %>% - dplyr::group_by(Site, Site_bis, FrSi_Code.Site, Site.date.sortie, Annee_semester) %>% - dplyr::summarise(FrSi_Nb.Total = mean(FrSi_Nb.Total), nb = dplyr::n()) -FrSi <- data.frame(FrSi) +frsi <- na.omit(freq_[, c("Site", "Site_bis", "frsi_Code.Site", "Site.date.sortie", "Annee_semester", "frsi_Nb.Total")]) %>% + dplyr::group_by(Site, Site_bis, frsi_Code.Site, Site.date.sortie, Annee_semester) %>% + dplyr::summarise(frsi_Nb.Total = mean(frsi_Nb.Total), nb = dplyr::n()) +frsi <- data.frame(frsi) -FrSsSi <- na.omit(freq.[, c("Site", "Site_bis", "Code.Sous.Site", "Site.date.sortie", "Annee_semester", "FrSsSi_Nb.Total")]) %>% +frsssi <- na.omit(freq_[, c("Site", "Site_bis", "Code.Sous.Site", "Site.date.sortie", "Annee_semester", "frsssi_Nb.Total")]) %>% dplyr::group_by(Site, Site_bis, Code.Sous.Site, Site.date.sortie, Annee_semester) %>% - dplyr::summarise(FrSsSi_Nb.Total = mean(FrSsSi_Nb.Total), nb = dplyr::n()) -FrSsSi <- data.frame(FrSsSi) + dplyr::summarise(frsssi_Nb.Total = mean(frsssi_Nb.Total), nb = dplyr::n()) +frsssi <- data.frame(frsssi) -FrHa <- na.omit(freq.[, c("Site", "Site_bis", "Code.Habitat", "Site.date.sortie", "Annee_semester", "FrHa_Nb.Total")]) %>% +frha <- na.omit(freq_[, c("Site", "Site_bis", "Code.Habitat", "Site.date.sortie", "Annee_semester", "frha_Nb.Total")]) %>% dplyr::group_by(Site, Site_bis, Code.Habitat, Site.date.sortie, Annee_semester) %>% - dplyr::summarise(FrHa_Nb.Total = mean(FrHa_Nb.Total), nb = dplyr::n()) -FrHa <- data.frame(FrHa) + dplyr::summarise(frha_Nb.Total = mean(frha_Nb.Total), nb = dplyr::n()) +frha <- data.frame(frha) -FrSi.stat <- na.omit(FrSi[, c("Site", "Site_bis", "FrSi_Code.Site", "Site.date.sortie", "Annee_semester", "FrSi_Nb.Total")]) %>% - dplyr::group_by(Site, Site_bis, FrSi_Code.Site, Annee_semester) %>% - dplyr::summarise(mean.FrSi_Nb.Total = mean(FrSi_Nb.Total), median.FrSi_Nb.Total = median(FrSi_Nb.Total), min.FrSi_Nb.Total = min(FrSi_Nb.Total), max.FrSi_Nb.Total = max(FrSi_Nb.Total), nb.FrSi_Nb.Total = dplyr::n()) -FrSi.stat <- data.frame(FrSi.stat) +frsi_stat <- na.omit(frsi[, c("Site", "Site_bis", "frsi_Code.Site", "Site.date.sortie", "Annee_semester", "frsi_Nb.Total")]) %>% + dplyr::group_by(Site, Site_bis, frsi_Code.Site, Annee_semester) %>% + dplyr::summarise(mean.frsi_Nb.Total = mean(frsi_Nb.Total), median.frsi_Nb.Total = median(frsi_Nb.Total), min.frsi_Nb.Total = min(frsi_Nb.Total), max.frsi_Nb.Total = max(frsi_Nb.Total), nb.frsi_Nb.Total = dplyr::n()) +frsi_stat <- data.frame(frsi_stat) -FrSsSi.stat <- na.omit(FrSsSi[, c("Site", "Site_bis", "Code.Sous.Site", "Site.date.sortie", "Annee_semester", "FrSsSi_Nb.Total")]) %>% +frsssi_stat <- na.omit(frsssi[, c("Site", "Site_bis", "Code.Sous.Site", "Site.date.sortie", "Annee_semester", "frsssi_Nb.Total")]) %>% dplyr::group_by(Site, Site_bis, Code.Sous.Site, Annee_semester) %>% - dplyr::summarise(mean.FrSsSi_Nb.Total = mean(FrSsSi_Nb.Total), median.FrSsSi_Nb.Total = median(FrSsSi_Nb.Total), min.FrSsSi_Nb.Total = min(FrSsSi_Nb.Total), max.FrSsSi_Nb.Total = max(FrSsSi_Nb.Total), nb.FrSsSi_Nb.Total = dplyr::n()) -FrSsSi.stat <- data.frame(FrSsSi.stat) + dplyr::summarise(mean.frsssi_Nb.Total = mean(frsssi_Nb.Total), median.frsssi_Nb.Total = median(frsssi_Nb.Total), min.frsssi_Nb.Total = min(frsssi_Nb.Total), max.frsssi_Nb.Total = max(frsssi_Nb.Total), nb.frsssi_Nb.Total = dplyr::n()) +frsssi_stat <- data.frame(frsssi_stat) -FrHa.stat <- na.omit(FrHa[, c("Site", "Site_bis", "Code.Habitat", "Site.date.sortie", "Annee_semester", "FrHa_Nb.Total")]) %>% +frha_stat <- na.omit(frha[, c("Site", "Site_bis", "Code.Habitat", "Site.date.sortie", "Annee_semester", "frha_Nb.Total")]) %>% dplyr::group_by(Site, Site_bis, Code.Habitat, Annee_semester) %>% - dplyr::summarise(mean.FrHa_Nb.Total = mean(FrHa_Nb.Total), median.FrHa_Nb.Total = median(FrHa_Nb.Total), min.FrHa_Nb.Total = min(FrHa_Nb.Total), max.FrHa_Nb.Total = max(FrHa_Nb.Total), nb.FrHa_Nb.Total = dplyr::n()) -FrHa.stat <- data.frame(FrHa.stat) + dplyr::summarise(mean.frha_Nb.Total = mean(frha_Nb.Total), median.frha_Nb.Total = median(frha_Nb.Total), min.frha_Nb.Total = min(frha_Nb.Total), max.frha_Nb.Total = max(frha_Nb.Total), nb.frha_Nb.Total = dplyr::n()) +frha_stat <- data.frame(frha_stat) -freq.Nb.Total.stat <- dplyr::full_join(FrSi.stat, FrSsSi.stat, by = c("Site", "Site_bis", "Annee_semester")) -freq.Nb.Total.stat <- dplyr::full_join(freq.Nb.Total.stat, FrHa.stat, by = c("Site", "Site_bis", "Annee_semester")) +freq_nb_total_stat <- dplyr::full_join(frsi_stat, frsssi_stat, by = c("Site", "Site_bis", "Annee_semester")) +freq_nb_total_stat <- dplyr::full_join(freq_nb_total_stat, frha_stat, by = c("Site", "Site_bis", "Annee_semester")) #missing Annee_semester value and site value; add dummy observations to fill in gaps -freq.Nb.Total.stat[nrow(freq.Nb.Total.stat)+2, ] <- NA -freq.Nb.Total.stat[nrow(freq.Nb.Total.stat)-1, "Annee_semester"] <- "2012_s1" -freq.Nb.Total.stat[nrow(freq.Nb.Total.stat), "Annee_semester"] <- "2017_s1" - -freq.Nb.Total.stat$Annee_semester_to.nb <- freq.Nb.Total.stat$Annee_semester -freq.Nb.Total.stat$Annee_semester_to.nb <- as.factor(freq.Nb.Total.stat$Annee_semester_to.nb) +freq_nb_total_stat[nrow(freq_nb_total_stat) + 2, ] <- NA +freq_nb_total_stat[nrow(freq_nb_total_stat) - 1, "Annee_semester"] <- "2012_s1" +freq_nb_total_stat[nrow(freq_nb_total_stat), "Annee_semester"] <- "2017_s1" -levels(freq.Nb.Total.stat$Annee_semester_to.nb) <- c(1:length(sort(unique(freq.Nb.Total.stat$Annee_semester_to.nb)))) -freq.Nb.Total.stat$Annee_semester_to.nb <- as.numeric(freq.Nb.Total.stat$Annee_semester_to.nb) +freq_nb_total_stat$Annee_semester_to.nb <- freq_nb_total_stat$Annee_semester +freq_nb_total_stat$Annee_semester_to.nb <- as.factor(freq_nb_total_stat$Annee_semester_to.nb) -#saveRDS(freq.Nb.Total.stat, "Freq. & cptmt/freq.Nb.Total.stat.RDS") +levels(freq_nb_total_stat$Annee_semester_to.nb) <- c(1:length(sort(unique(freq_nb_total_stat$Annee_semester_to.nb)))) +freq_nb_total_stat$Annee_semester_to.nb <- as.numeric(freq_nb_total_stat$Annee_semester_to.nb) -op <- par(mar = c(8,4,4,2) + 0.1) ## default is c(5,4,4,2) + 0.1 +#saveRDS(freq_nb_total_stat, "freq_ & cptmt/freq_nb_total_stat.RDS") -for (i in c(1:length(na.omit(unique(freq.Nb.Total.stat[, "Site"]))))) { +op <- par(mar = c(8, 4, 4, 2) + 0.1) ## default is c(5,4,4,2) + 0.1 - # add na.omit cfr I had to add empty observations with missing factor levels - #i <- 1 - - dplyr::filter(freq.Nb.Total.stat, Site == unique(freq.Nb.Total.stat[, "Site"])[i]) -> freq.i - - # I have to invert plot(x,y) by plot(y,x) in order to work ... ? - png(paste0("freq_2_", unique(freq.i$Site), ".png")) - plot(freq.i$Annee_semester_to.nb, rep(0, length(freq.i$Annee_semester_to.nb)), - xlim = c(1, length(unique(freq.Nb.Total.stat$Annee_semester_to.nb))), - ylim = c(-5, 200)#c(0-( - #(max(c(freq.i[, "max.FrSi_Nb.Total"], freq.i[, "max.FrSsSi_Nb.Total"], freq.i[, "max.FrHa_Nb.Total"]), na.rm = T)*1.05) - - # (max(c(freq.i[, "max.FrSi_Nb.Total"], freq.i[, "max.FrSsSi_Nb.Total"], freq.i[, "max.FrHa_Nb.Total"]), na.rm = T))), - #(max(c(freq.i[, "max.FrSi_Nb.Total"], freq.i[, "max.FrSsSi_Nb.Total"], freq.i[, "max.FrHa_Nb.Total"]), na.rm = T)*1.05) ) - - , main = unique(freq.i$Site), xlab = "", ylab = "fréquentation", xaxt = "n", col = "white") - - axis(1, at = 1:length(unique(freq.Nb.Total.stat$Annee_semester_to.nb)), labels = sort(unique(freq.Nb.Total.stat$Annee_semester)), las = 2) - - #title(xlab = "Annee_semester", cex.lab = 1, line = 5) - - dplyr::filter(freq.i, nb.FrSi_Nb.Total >= 2) -> freq.i.red - if (nrow(freq.i.red) >= 1) { - arrows(as.numeric(freq.i.red$Annee_semester_to.nb), freq.i.red$median.FrSi_Nb.Total, as.numeric(freq.i.red$Annee_semester_to.nb), freq.i.red$max.FrSi_Nb.Total, code = 3, angle = 90, length = 0.00) - } else { - points(0, 0, col = "white") - } - - if (nrow(freq.i.red) >= 1) { - arrows(as.numeric(freq.i.red$Annee_semester_to.nb), freq.i.red$median.FrSi_Nb.Total, as.numeric(freq.i.red$Annee_semester_to.nb), freq.i.red$min.FrSi_Nb.Total, code = 3, angle = 90, length = 0.00) - } else { - points(0, 0, col = "white") - } - - points(freq.i$Annee_semester_to.nb, freq.i$mean.FrSi_Nb.Total, pch = 19, col = "orange") - points(freq.i$Annee_semester_to.nb, freq.i$median.FrSi_Nb.Total, pch = 19, col = "red") - - rm(freq.i.red) - - dplyr::filter(freq.i, nb.FrSsSi_Nb.Total >= 2) -> freq.i.red - - if (nrow(freq.i.red) >= 1) { - arrows(as.numeric(freq.i.red$Annee_semester_to.nb)+.25, freq.i.red$median.FrSsSi_Nb.Total, as.numeric(freq.i.red$Annee_semester_to.nb)+.25, freq.i.red$max.FrSsSi_Nb.Total, code = 3, angle = 90, length = 0.00) - } else { - points(0, 0, col = "white") - } - - if (nrow(freq.i.red) >= 1) { - arrows(as.numeric(freq.i.red$Annee_semester_to.nb)+.25, freq.i.red$median.FrSsSi_Nb.Total, as.numeric(freq.i.red$Annee_semester_to.nb)+.25, freq.i.red$min.FrSsSi_Nb.Total, code = 3, angle = 90, length = 0.00) - } else { - points(0, 0, col = "white") - } - - points(as.numeric(freq.i$Annee_semester_to.nb) + .25, freq.i$mean.FrSsSi_Nb.Total, pch = 19, col = "lightblue") - points(as.numeric(freq.i$Annee_semester_to.nb) + .25, freq.i$median.FrSsSi_Nb.Total, pch = 19, col = "darkblue") - - rm(freq.i.red) - - dplyr::filter(freq.i, nb.FrHa_Nb.Total >= 2) -> freq.i.red - - if (nrow(freq.i.red) >= 1) { - arrows(as.numeric(freq.i.red$Annee_semester_to.nb)-.25, freq.i.red$median.FrHa_Nb.Total, as.numeric(freq.i.red$Annee_semester_to.nb)-.25, freq.i.red$max.FrHa_Nb.Total, code = 3, angle = 90, length = 0.00) - } else { - points(0, 0, col = "white") - } - - if (nrow(freq.i.red) >= 1) { - arrows(as.numeric(freq.i.red$Annee_semester_to.nb)-.25, freq.i.red$median.FrHa_Nb.Total, as.numeric(freq.i.red$Annee_semester_to.nb)-.25, freq.i.red$min.FrHa_Nb.Total, code = 3, angle = 90, length = 0.00) - } else { - points(0, 0, col = "white") - } - - points(as.numeric(freq.i$Annee_semester_to.nb) - .25, freq.i$mean.FrHa_Nb.Total, pch = 19, col = "lightgreen") - points(as.numeric(freq.i$Annee_semester_to.nb) - .25, freq.i$median.FrHa_Nb.Total, pch = 19, col = "forestgreen") - - rm(freq.i.red) - - legend("bottom", - #inset = c(0,-0.75), - inset = c(0,-0.5), - legend = c("méd. Si", "moy. Si","méd. SsSi", "moy. SsSi", "méd. Ha", "moy. Ha"), pch = c(19,19,19,19,19,19), col = c("red", "orange", "darkblue", "lightblue", "forestgreen", "lightgreen")#, horiz = TRUE - , bty = "n", xpd = TRUE, - ncol = 3) - #inset = c(-0.45, 0) # You will need to fine-tune the first value depending on the windows size - # xpd = TRUE # You need to specify this graphical parameter to put the legend outside the plot - -} # par(mfrow = c(1,1)) par(op) # to be able to run the Rmarkdown loop -freq.$Site <- as.factor(freq.$Site) -freq.$Site <- ordered(freq.$Site) -freq.$Site.nb <- freq.$Site -levels(freq.$Site.nb) <- 1:length(unique(freq.$Site)) -freq.$Site.nb <- as.integer(freq.$Site.nb) - +freq_$Site <- as.factor(freq_$Site) +freq_$Site <- ordered(freq_$Site) +freq_$Site.nb <- freq_$Site +levels(freq_$Site.nb) <- 1:length(unique(freq_$Site)) +freq_$Site.nb <- as.integer(freq_$Site.nb) -freq.$Site.nb <- ifelse(freq.$Site %in% c("ARMO_Piegu", "ARMO_Verdelet") == TRUE, - unique(dplyr::filter(freq., Site == "ARMO_Piégu / Verdelet")[,"Site.nb"]), - freq.$Site.nb) -freq.$Site.nb <- ifelse(freq.$Site %in% c("BASQ_FlotsBleusZF", "BASQ_FlotsBleusZP"), - unique(dplyr::filter(freq., Site == "BASQ_FlotsBleus")[, "Site.nb"]), - freq.$Site.nb) +freq_$Site.nb <- ifelse(freq_$Site %in% c("ARMO_Piegu", "ARMO_Verdelet") == TRUE, + unique(dplyr::filter(freq_, Site == "ARMO_Piégu / Verdelet")[,"Site.nb"]), + freq_$Site.nb) +freq_$Site.nb <- ifelse(freq_$Site %in% c("BASQ_FlotsBleusZF", "BASQ_FlotsBleusZP"), + unique(dplyr::filter(freq_, Site == "BASQ_FlotsBleus")[, "Site.nb"]), + freq_$Site.nb) -saveRDS(freq., "freq.RDS") +saveRDS(freq_, "freq.RDS") -freq.Nb.Total.stat$Site <- as.factor(freq.Nb.Total.stat$Site) -freq.Nb.Total.stat$Site <- ordered(freq.Nb.Total.stat$Site) -freq.Nb.Total.stat$Site.nb <- freq.Nb.Total.stat$Site -levels(freq.Nb.Total.stat$Site.nb) <- 1:length(unique(freq.Nb.Total.stat$Site)) -freq.Nb.Total.stat$Site.nb <- as.integer(freq.Nb.Total.stat$Site.nb) -#freq.Nb.Total.stat$Site <- as.character(freq.Nb.Total.stat$Site) -freq.Nb.Total.stat$Site.nb <- ifelse(freq.Nb.Total.stat$Site %in% c("ARMO_Piegu", "ARMO_Verdelet") == TRUE, - unique(dplyr::filter(freq.Nb.Total.stat, Site == "ARMO_Piégu / Verdelet")[,"Site.nb"]), - freq.Nb.Total.stat$Site.nb) -freq.Nb.Total.stat$Site.nb <- ifelse(freq.Nb.Total.stat$Site %in% c("BASQ_FlotsBleusZF", "BASQ_FlotsBleusZP"), - unique(dplyr::filter(freq.Nb.Total.stat, Site == "BASQ_FlotsBleus")[, "Site.nb"]), - freq.Nb.Total.stat$Site.nb) +freq_nb_total_stat$Site <- as.factor(freq_nb_total_stat$Site) +freq_nb_total_stat$Site <- ordered(freq_nb_total_stat$Site) +freq_nb_total_stat$Site.nb <- freq_nb_total_stat$Site +levels(freq_nb_total_stat$Site.nb) <- 1:length(unique(freq_nb_total_stat$Site)) +freq_nb_total_stat$Site.nb <- as.integer(freq_nb_total_stat$Site.nb) +freq_nb_total_stat$Site.nb <- ifelse(freq_nb_total_stat$Site %in% c("ARMO_Piegu", "ARMO_Verdelet") == TRUE, + unique(dplyr::filter(freq_nb_total_stat, Site == "ARMO_Piégu / Verdelet")[,"Site.nb"]), + freq_nb_total_stat$Site.nb) +freq_nb_total_stat$Site.nb <- ifelse(freq_nb_total_stat$Site %in% c("BASQ_FlotsBleusZF", "BASQ_FlotsBleusZP"), + unique(dplyr::filter(freq_nb_total_stat, Site == "BASQ_FlotsBleus")[, "Site.nb"]), + freq_nb_total_stat$Site.nb) -#saveRDS(freq.Nb.Total.stat, "freq.Nb.Total.stat.RDS") ## Work on comportment data now -obs. <- read.csv2(input_obs, fileEncoding = "Latin1") +obs_ <- read.csv2(input_obs, fileEncoding = "Latin1") -fiche <- read.csv2(fiche_val, fileEncoding = "Latin1") +fiche <- read.csv2(fiche_2, fileEncoding = "Latin1") -obs. <- dplyr::left_join(obs., fiche[, c("ID.Fiche", "date.sortie", "libellé.campagne", "libellé.sortie", "territoire", "code.site", "site", "sous.site", "zone.habitat", "type.protocole", "version.protocole")], by = "ID.Fiche") +obs_ <- dplyr::left_join(obs_, fiche) -obs.$date.sortie <- as.Date(obs.$date.sortie, origin = "1970-01-01") +obs_$date.sortie <- as.Date(obs_$date.sortie, origin = "1970-01-01") -obs.$Heure.Debut <- ifelse(obs.$Heure.Debut == "", NA, obs.$Heure.Debut) -obs.$Heure.Fin <- ifelse(obs.$Heure.Debut == "", NA, obs.$Heure.Fin) -obs. <- tibble::add_column(obs., Tps.obs = difftime(strptime(obs.$Heure.Fin, format = "%H:%M"), strptime(obs.$Heure.Debut, format = "%H:%M"), units = "mins"), .after = "Heure.Fin") -table(obs.$Tps.obs)["15"] -obs.$Tps.obs <- as.integer(obs.$Tps.obs) +obs_$Heure.Debut <- ifelse(obs_$Heure.Debut == "", NA, obs_$Heure.Debut) +obs_$Heure.Fin <- ifelse(obs_$Heure.Debut == "", NA, obs_$Heure.Fin) +obs_ <- tibble::add_column(obs_, Tps.obs = difftime(strptime(obs_$Heure.Fin, format = "%H:%M"), strptime(obs_$Heure.Debut, format = "%H:%M"), units = "mins"), .after = "Heure.Fin") +table(obs_$Tps.obs)["15"] +obs_$Tps.obs <- as.integer(obs_$Tps.obs) -qecbNew <- readRDS(input_qecb) +qecbnew <- readRDS(input_qecb) -# insert site name based on qecb df. +# insert site name based on qecb df -(df. <- unique(qecbNew[, c("code.site", "Site", "Site_bis")])) +(df <- unique(qecbnew[, c("code.site", "Site", "Site_bis")])) -(df.join <- unique(obs.[, c("code.site", "zone.habitat" +(df_join <- unique(obs_[, c("code.site", "zone.habitat" )])) -(df.join <- dplyr::left_join(df.join, df., by = "code.site")) -df.join <- df.join[!(df.join$code.site == "ARMO_042-043" & df.join$Site == "ARMO_Piegu"),] -df.join <- df.join[!(df.join$zone.habitat == "Les Flots Bleus (champ de blocs) - zone pècheurs" & df.join$Site == "BASQ_FlotsBleusZF"),] -df.join <- df.join[!(df.join$zone.habitat == "Les Flots Bleus (champ de blocs) - zone famille" & df.join$Site == "BASQ_FlotsBleusZP"),] +(df_join <- dplyr::left_join(df_join, df, by = "code.site")) +df_join <- df_join[!(df_join$code.site == "ARMO_042-043" & df_join$Site == "ARMO_Piegu"),] +df_join <- df_join[!(df_join$zone.habitat == "Les Flots Bleus (champ de blocs) - zone pècheurs" & df_join$Site == "BASQ_FlotsBleusZF"),] +df_join <- df_join[!(df_join$zone.habitat == "Les Flots Bleus (champ de blocs) - zone famille" & df_join$Site == "BASQ_FlotsBleusZP"),] library(data.table) -for (i in c(1:nrow(df.join))) { +for (i in c(1:nrow(df_join))) { #i <- 1 - setDT(obs.)[code.site == df.join[i, "code.site"], Site := df.join[i,"Site"]] + setDT(obs_)[code.site == df_join[i, "code.site"], Site := df_join[i,"Site"]] } -obs. <- data.frame(obs.) +obs_ <- data.frame(obs_) -obs.$Site <- ifelse(obs.$code.site == "BASQ_01" & obs.$zone.habitat == "Les Flots Bleus (champ de blocs) - zone pècheurs", "BASQ_FlotsBleusZP", obs.$Site) +obs_$Site <- ifelse(obs_$code.site == "BASQ_01" & obs_$zone.habitat == "Les Flots Bleus (champ de blocs) - zone pècheurs", "BASQ_FlotsBleusZP", obs_$Site) -rm(df.join) -rm(df.) +rm(df_join) +rm(df) -obs. <- tibble::add_column(obs., Site = obs.$Site, .after = "code.site") -obs. <- obs.[ , !names(obs.) %in% c("Site")] -obs. <- dplyr::rename(obs., Site = Site.1) +obs_ <- tibble::add_column(obs_, Site = obs_$Site, .after = "code.site") +obs_ <- obs_[ , !names(obs_) %in% c("Site")] +obs_ <- dplyr::rename(obs_, Site = Site.1) # There are some missing time data, either because there was no fishers, either because data were not encoded. Has to be completed. -obs.$Tps.obs <- ifelse(!is.na(obs.$Heure.Debut) & is.na(obs.$Tps.obs), 15, obs.$Tps.obs) +obs_$Tps.obs <- ifelse(!is.na(obs_$Heure.Debut) & is.na(obs_$Tps.obs), 15, obs_$Tps.obs) -obs. <- tibble::add_column(obs., Nb.Blocs.Retournes.remis.norm15min = obs.$Nb.Blocs.Retournes.remis/obs.$Tps.obs*15, .before = "date.sortie") -obs. <- tibble::add_column(obs., Nb.Blocs.Deplaces.non.remis.norm15min = obs.$Nb.Blocs.Deplaces.non.remis/obs.$Tps.obs*15, .before = "date.sortie") -obs. <- tibble::add_column(obs., Nb.Blocs.Retournes.non.remis.norm15min = obs.$Nb.Blocs.Retournes.non.remis/obs.$Tps.obs*15, .before = "date.sortie") +obs_ <- tibble::add_column(obs_, Nb.Blocs.Retournes.remis.norm15min = obs_$Nb.Blocs.Retournes.remis/obs_$Tps.obs * 15, .before = "date.sortie") +obs_ <- tibble::add_column(obs_, Nb.Blocs.Deplaces.non.remis.norm15min = obs_$Nb.Blocs.Deplaces.non.remis/obs_$Tps.obs * 15, .before = "date.sortie") +obs_ <- tibble::add_column(obs_, Nb.Blocs.Retournes.non.remis.norm15min = obs_$Nb.Blocs.Retournes.non.remis/obs_$Tps.obs * 15, .before = "date.sortie") # some issues with encoded data -dplyr::filter(obs., is.na(obs.$Nb.Blocs.Retournes.remis) | is.na(obs.$Nb.Blocs.Retournes.remis.norm15min) | is.na(obs.$Nb.Blocs.Retournes.non.remis) | is.na(obs.$Nb.Blocs.Retournes.non.remis.norm15min) | is.na(obs.$Nb.Blocs.Deplaces.non.remis) | is.na(obs.$Nb.Blocs.Deplaces.non.remis.norm15min)) -> check -obs.$Nb.Blocs.Retournes.non.remis <- ifelse(obs.$ID.Fiche == 373166, 0, obs.$Nb.Blocs.Retournes.non.remis) -obs.$Nb.Blocs.Retournes.non.remis.norm15min <- ifelse(obs.$ID.Fiche == 373166, 0, obs.$Nb.Blocs.Retournes.non.remis.norm15min) -obs.$Nb.Blocs.Deplaces.non.remis <- ifelse(obs.$ID.Fiche == 373166, 0, obs.$Nb.Blocs.Deplaces.non.remis) -obs.$Nb.Blocs.Deplaces.non.remis.norm15min <- ifelse(obs.$ID.Fiche == 373166, 0, obs.$Nb.Blocs.Deplaces.non.remis.norm15min) -dplyr::filter(obs., is.na(obs.$Nb.Blocs.Retournes.remis) | is.na(obs.$Nb.Blocs.Retournes.remis.norm15min) | is.na(obs.$Nb.Blocs.Retournes.non.remis) | is.na(obs.$Nb.Blocs.Retournes.non.remis.norm15min) | is.na(obs.$Nb.Blocs.Deplaces.non.remis) | is.na(obs.$Nb.Blocs.Deplaces.non.remis.norm15min)) -> check -obs. <- dplyr::filter(obs., ID.Fiche %notin% c(check$ID.Fiche)) +check <- dplyr::filter(obs_, is.na(obs_$Nb.Blocs.Retournes.remis) | is.na(obs_$Nb.Blocs.Retournes.remis.norm15min) | is.na(obs_$Nb.Blocs.Retournes.non.remis) | is.na(obs_$Nb.Blocs.Retournes.non.remis.norm15min) | is.na(obs_$Nb.Blocs.Deplaces.non.remis) | is.na(obs_$Nb.Blocs.Deplaces.non.remis.norm15min)) +obs_$Nb.Blocs.Retournes.non.remis <- ifelse(obs_$ID.Fiche == 373166, 0, obs_$Nb.Blocs.Retournes.non.remis) +obs_$Nb.Blocs.Retournes.non.remis.norm15min <- ifelse(obs_$ID.Fiche == 373166, 0, obs_$Nb.Blocs.Retournes.non.remis.norm15min) +obs_$Nb.Blocs.Deplaces.non.remis <- ifelse(obs_$ID.Fiche == 373166, 0, obs_$Nb.Blocs.Deplaces.non.remis) +obs_$Nb.Blocs.Deplaces.non.remis.norm15min <- ifelse(obs_$ID.Fiche == 373166, 0, obs_$Nb.Blocs.Deplaces.non.remis.norm15min) +check <- dplyr::filter(obs_, is.na(obs_$Nb.Blocs.Retournes.remis) | is.na(obs_$Nb.Blocs.Retournes.remis.norm15min) | is.na(obs_$Nb.Blocs.Retournes.non.remis) | is.na(obs_$Nb.Blocs.Retournes.non.remis.norm15min) | is.na(obs_$Nb.Blocs.Deplaces.non.remis) | is.na(obs_$Nb.Blocs.Deplaces.non.remis.norm15min)) +obs_ <- dplyr::filter(obs_, ID.Fiche %notin% c(check$ID.Fiche)) # add replicate numbers for fishers observed during the same survey -obs. <- tibble::add_column(obs., Site.date.sortie = paste0(obs.$Site, ".", obs.$date.sortie), .after = "Site") +obs_ <- tibble::add_column(obs_, Site.date.sortie = paste0(obs_$Site, ".", obs_$date.sortie), .after = "Site") -obs. <- obs. %>% +obs_ <- obs_ %>% dplyr::group_by(Site.date.sortie) %>% dplyr::mutate(repl. = dplyr::row_number(Site.date.sortie)) -obs. <- data.frame(obs.) +obs_ <- data.frame(obs_) -# add a unique variable for the three obs. cptmt., with ponderation. +# add a unique variable for the three obs_ cptmt., with ponderation. -ratio. <- ((obs.$Nb.Blocs.Retournes.remis*0)+(obs.$Nb.Blocs.Deplaces.non.remis*0.5)+(obs.$Nb.Blocs.Retournes.non.remis*1))/(obs.$Nb.Blocs.Retournes.remis+obs.$Nb.Blocs.Deplaces.non.remis+obs.$Nb.Blocs.Retournes.non.remis) -obs. <- tibble::add_column(obs., ratio.obs = ratio., .after = "Nb.Blocs.Retournes.non.remis") -obs.$ratio.obs <- ifelse(obs.$Nb.Blocs.Retournes.remis == 0 & obs.$Nb.Blocs.Deplaces.non.remis == 0 & obs.$Nb.Blocs.Retournes.non.remis == 0, 0, obs.$ratio.obs) +ratio_ <- ((obs_$Nb.Blocs.Retournes.remis * 0) + (obs_$Nb.Blocs.Deplaces.non.remis * 0.5) + (obs_$Nb.Blocs.Retournes.non.remis * 1)) /(obs_$Nb.Blocs.Retournes.remis+obs_$Nb.Blocs.Deplaces.non.remis+obs_$Nb.Blocs.Retournes.non.remis) +obs_ <- tibble::add_column(obs_, ratio_obs = ratio_, .after = "Nb.Blocs.Retournes.non.remis") +obs_$ratio_obs <- ifelse(obs_$Nb.Blocs.Retournes.remis == 0 & obs_$Nb.Blocs.Deplaces.non.remis == 0 & obs_$Nb.Blocs.Retournes.non.remis == 0, 0, obs_$ratio_obs) -ratio.norm15min <- ((obs.$Nb.Blocs.Retournes.remis.norm15min*0)+(obs.$Nb.Blocs.Deplaces.non.remis.norm15min*0.5)+(obs.$Nb.Blocs.Retournes.non.remis.norm15min*1))/(obs.$Nb.Blocs.Retournes.remis.norm15min+obs.$Nb.Blocs.Deplaces.non.remis.norm15min+obs.$Nb.Blocs.Retournes.non.remis.norm15min) -obs. <- tibble::add_column(obs., ratio.norm15min.obs = ratio.norm15min, .after = "Nb.Blocs.Retournes.non.remis.norm15min") -obs.$ratio.norm15min.obs <- ifelse(obs.$Nb.Blocs.Retournes.remis.norm15min == 0 & obs.$Nb.Blocs.Deplaces.non.remis.norm15min == 0 & obs.$Nb.Blocs.Retournes.non.remis.norm15min == 0, 0, obs.$ratio.norm15min.obs) +ratio_norm15min <- ((obs_$Nb.Blocs.Retournes.remis.norm15min * 0) + (obs_$Nb.Blocs.Deplaces.non.remis.norm15min * 0.5) + (obs_$Nb.Blocs.Retournes.non.remis.norm15min * 1)) /(obs_$Nb.Blocs.Retournes.remis.norm15min+obs_$Nb.Blocs.Deplaces.non.remis.norm15min+obs_$Nb.Blocs.Retournes.non.remis.norm15min) +obs_ <- tibble::add_column(obs_, ratio_norm15min.obs = ratio_norm15min, .after = "Nb.Blocs.Retournes.non.remis.norm15min") +obs_$ratio_norm15min.obs <- ifelse(obs_$Nb.Blocs.Retournes.remis.norm15min == 0 & obs_$Nb.Blocs.Deplaces.non.remis.norm15min == 0 & obs_$Nb.Blocs.Retournes.non.remis.norm15min == 0, 0, obs_$ratio_norm15min.obs) -saveRDS(obs., "obs.RDS") +saveRDS(obs_, "obs.RDS") -xmin. <- as.Date("2013-01-01", origin = "1970-01-01") -xmax. <- as.Date("2022-01-01", origin = "1970-01-01") +xmin_ <- as.Date("2013-01-01", origin = "1970-01-01") +xmax_ <- as.Date("2022-01-01", origin = "1970-01-01") -ymax. = max(c(obs.$Nb.Blocs.Deplaces.non.remis, obs.$Nb.Blocs.Deplaces.non.remis.norm15min, obs.$Nb.Blocs.Retournes.non.remis, obs.$Nb.Blocs.Retournes.non.remis.norm15min, obs.$Nb.Blocs.Retournes.remis, obs.$Nb.Blocs.Retournes.remis.norm15min), na.rm = T) +ymax_ = max(c(obs_$Nb.Blocs.Deplaces.non.remis, obs_$Nb.Blocs.Deplaces.non.remis.norm15min, obs_$Nb.Blocs.Retournes.non.remis, obs_$Nb.Blocs.Retournes.non.remis.norm15min, obs_$Nb.Blocs.Retournes.remis, obs_$Nb.Blocs.Retournes.remis.norm15min), na.rm = TRUE) -par(mfrow = c(1,1)) +par(mfrow = c(1, 1)) -for (i in c(1:length(unique(obs.$Site)))) { +for (i in c(1:length(unique(obs_$Site)))) { #i <- 5 - dplyr::filter(obs., Site == unique(obs.$Site)[i]) -> obs.i + obs_i <- dplyr::filter(obs_, Site == unique(obs_$Site)[i]) - ymax. = max(c(obs.i$Nb.Blocs.Deplaces.non.remis, obs.i$Nb.Blocs.Deplaces.non.remis.norm15min, obs.i$Nb.Blocs.Retournes.non.remis, obs.i$Nb.Blocs.Retournes.non.remis.norm15min, obs.i$Nb.Blocs.Retournes.remis, obs.i$Nb.Blocs.Retournes.remis.norm15min), na.rm = T) + ymax_ <- max(c(obs_i$Nb.Blocs.Deplaces.non.remis, obs_i$Nb.Blocs.Deplaces.non.remis.norm15min, obs_i$Nb.Blocs.Retournes.non.remis, obs_i$Nb.Blocs.Retournes.non.remis.norm15min, obs_i$Nb.Blocs.Retournes.remis, obs_i$Nb.Blocs.Retournes.remis.norm15min), na.rm = TRUE) - png(paste0("Observations_1", unique(obs.i$Site), ".png")) - plot(obs.i$date.sortie, obs.i$Nb.Blocs.Retournes.remis.norm15min, xlim = c(xmin.,xmax.), ylim = c(0,105)#ymax.*1.05) -, xlab = "", ylab = "Nbr de blocs (normalisé pour 15 min d'obs.)", main = unique(obs.i$Site), col = "darkgreen") - points(obs.i$date.sortie, obs.i$Nb.Blocs.Retournes.non.remis.norm15min, col = "red") - points(obs.i$date.sortie, obs.i$Nb.Blocs.Deplaces.non.remis.norm15min, col = "orange") - - legend("bottom", inset = c(0,-0.4), legend = c("ret.re.","ret.n.re.", "dépl.n.re."), pch = c(19,19,19), col = c("darkgreen", "red", "orange"), horiz = TRUE, bty = "n", xpd = TRUE) - -} + obs_plot <- ggplot2::ggplot() + + #ggplot2::geom_point(ggplot2::aes(x = obs_$date.sortie, y = obs_$$Nb.Blocs.Retournes.remis.norm15min), col = "Other sites good behaviour", alpha = 0.2) + + ggplot2::geom_point(ggplot2::aes(obs_i$date.sortie, obs_i$Nb.Blocs.Retournes.remis.norm15min, colour = "Good : put the bloc back")) + + ggplot2::geom_point(ggplot2::aes(obs_i$date.sortie, obs_i$Nb.Blocs.Retournes.non.remis.norm15min, colour = "Bad : turned bloc")) + + ggplot2::geom_point(ggplot2::aes(obs_i$date.sortie, obs_i$Nb.Blocs.Deplaces.non.remis.norm15min, colour = "Bad : moved bloc")) + + #ggplot2::geom_point(ggplot2::aes(x = obs_$date.sortie, y = obs_$Nb.Blocs.Retournes.remis.norm15min), col = "grey") + + #ggplot2::geom_point(ggplot2::aes(x = obs_$date.sortie, y = obs_$Nb.Blocs.Retournes.non.remis.norm15min), col = "grey") + + #ggplot2::geom_point(ggplot2::aes(x = obs_$date.sortie, y = obs_$Nb.Blocs.Deplaces.non.remis.norm15min), col = "grey") + + ggplot2::scale_fill_manual(values = c("#35B117", "#E12D07", "#F9A806"), name = "Behaviour", breaks = c("Good : put the bloc back", "Bad : turned bloc", "Bad : moved bloc")) + + ggplot2::ylab("Nbr de blocs (normalisé pour 15 min d'obs_)") + + ggplot2::xlab("Date") + + ggplot2::ggtitle(unique(obs_i$Site)) + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1), legend.title = ggplot2::element_text("Behaviour")) -for (i in c(1:length(unique(obs.$Site)))) { - - #i <- 5 - - dplyr::filter(obs., Site == unique(obs.$Site)[i]) -> obs.i - - png(paste0("Observations_2", unique(obs.i$Site), ".png")) - plot(obs.i$date.sortie, obs.i$ratio.norm15min.obs*100, xlim = c(xmin.,xmax.), ylim = c(0,105), xlab = "", ylab = "comportement relatif (normalisé pour 15 min d'obs.)", main = unique(obs.i$Site), col = "black") + #plot(obs_i$date.sortie, obs_i$Nb.Blocs.Retournes.remis.norm15min, xlim = c(xmin_,xmax_), ylim = c(0,105)#ymax_*1.05) +#, xlab = "", ylab = "Nbr de blocs (normalisé pour 15 min d'obs_)", main = unique(obs_i$Site), col = "darkgreen") + #points(obs_i$date.sortie, obs_i$Nb.Blocs.Retournes.non.remis.norm15min, col = "red") + #points(obs_i$date.sortie, obs_i$Nb.Blocs.Deplaces.non.remis.norm15min, col = "orange") + #legend("bottom", inset = c(0, -0.4), legend = c("ret.re.","ret.n.re.", "dépl.n.re."), pch = c(19, 19, 19), col = c("darkgreen", "red", "orange"), horiz = TRUE, bty = "n", xpd = TRUE) + ggplot2::ggsave(paste0("obs_", unique(obs_i$Site), ".png"), obs_plot, height = 4, width = 6) } diff --git a/cb_freq.xml b/cb_freq.xml index e20b386..ae4dfe6 100644 --- a/cb_freq.xml +++ b/cb_freq.xml @@ -9,6 +9,14 @@ r-tibble r-ggplot2 + + + + + + + + - - - - - - + + + + + + + - - + + + + + - - - - - + + + + + + + - + + + + + + - + + - + diff --git a/freq_graph b/freq_graph new file mode 100644 index 0000000..00a1c66 --- /dev/null +++ b/freq_graph @@ -0,0 +1,181 @@ +png(paste0("freq_1_", unique(freq_i$Site), ".png")) + plot(freq_i$date.sortie, rep(0, length(freq_i$date.sortie)), + xlim = c(xmin_, xmax_), + ylim = c(-8, 165) +#round(0-((max(#freq_i[, "frsi_Nb.Total"], + #freq_i[, "frsssi_Nb.Total"]#, freq_i[, "frha_Nb.Total"] +#, na.rm = TRUE)*1.05)-max(#freq_i[, "frsi_Nb.Total"], +#freq_i[, "frsssi_Nb.Total"]#, freq_i[, "frha_Nb.Total"] +#, na.rm = TRUE))), #log10 + # round(max(#freq_i[, "frsi_Nb.Total"], +#freq_i[, "frsssi_Nb.Total"]#, freq_i[, "frha_Nb.Total"] +#, na.rm = TRUE)*1.05)), + , main = unique(freq_i$Site), xlab = "", ylab = "fréquentation", col = "white") + + #points(c(freq_$frsi_Nb.Total, freq_$frsssi_Nb.Total, freq_$frha_Nb.Total), c(freq_$date.sortie, freq_$date.sortie, freq_$date.sortie), pch = 19, col = "lightgrey", cex = .5) + + points(freq_i$date.sortie, #log10 + (freq_i$frsi_Nb.Total + #+1 + ), pch = 19, col = "red") + points(freq_i$date.sortie, #log10 + (freq_i$frsssi_Nb.Total + #+1 + ), pch = 19, col = "darkblue") + points(freq_i$date.sortie, #log10 + (freq_i$frha_Nb.Total + #+1 + ), pch = 19, col = "forestgreen") + + legend("bottom", inset = c(0, -0.4), legend = c("Si", "SsSi", "Ha"), pch = c(19, 19, 19), col = c("red", "darkblue", "forestgreen"), horiz = TRUE, bty = "n", xpd = TRUE) + #inset = c(-0.45, 0) # You will need to fine-tune the first value depending on the windows size + # xpd = TRUE # You need to specify this graphical parameter to put the legend outside the plot + + + + + + + + + + + + +# I have to invert plot(x,y) by plot(y,x) in order to work ... ? + png(paste0("freq_2_", unique(freq_i$Site), ".png")) + plot(freq_i$Annee_semester_to.nb, rep(0, length(freq_i$Annee_semester_to.nb)), + xlim = c(1, length(unique(freq_nb_total_stat$Annee_semester_to.nb))), + ylim = c(-5, 200)#c(0-( + #(max(c(freq_i[, "max.frsi_Nb.Total"], freq_i[, "max.frsssi_Nb.Total"], freq_i[, "max.frha_Nb.Total"]), na.rm = TRUE)*1.05) - + # (max(c(freq_i[, "max.frsi_Nb.Total"], freq_i[, "max.frsssi_Nb.Total"], freq_i[, "max.frha_Nb.Total"]), na.rm = TRUE))), + #(max(c(freq_i[, "max.frsi_Nb.Total"], freq_i[, "max.frsssi_Nb.Total"], freq_i[, "max.frha_Nb.Total"]), na.rm = TRUE)*1.05) ) + + , main = unique(freq_i$Site), xlab = "", ylab = "fréquentation", xaxt = "n", col = "white") + + axis(1, at = 1:length(unique(freq_nb_total_stat$Annee_semester_to.nb)), labels = sort(unique(freq_nb_total_stat$Annee_semester)), las = 2) + + #title(xlab = "Annee_semester", cex.lab = 1, line = 5) + + + + + + + + +for (i in c(1:length(na.omit(unique(freq_nb_total_stat[, "Site"]))))) { + + # add na.omit cfr I had to add empty observations with missing factor levels + #i <- 1 + + freq_i <- dplyr::filter(freq_nb_total_stat, Site == unique(freq_nb_total_stat[, "Site"])[i]) + + freq_i_red <- dplyr::filter(freq_i, nb.frsi_Nb.Total >= 2) + + freq_plot <- ggplot2::ggplot() + + ggplot2::ylab("fréquentation") + + ggplot2::xlab("semestres") + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1), legend.position = "bottom") + + + + if (nrow(freq_i_red) >= 1) { + freq_plot <- freq_plot + ggplot2::geom_pointrange(ggplot2::aes(as.numeric(freq_i_red$Annee_semester_to.nb), freq_i_red$median.frsi_Nb.Total, as.numeric(freq_i_red$Annee_semester_to.nb), freq_i_red$max.frsi_Nb.Total), code = 3, angle = 90, length = 0.00) +#arrows(as.numeric(freq_i_red$Annee_semester_to.nb), freq_i_red$median.frsi_Nb.Total, as.numeric(freq_i_red$Annee_semester_to.nb), freq_i_red$max.frsi_Nb.Total, code = 3, angle = 90, length = 0.00) + } else { + freq_plot <- freq_plot + ggplot2::geom_point(ggplot2::aes(0, 0), col = "white") + #points(0, 0, col = "white") + } + + if (nrow(freq_i_red) >= 1) { + freq_plot <- freq_plot + ggplot2::geom_pointrange(ggplot2::aes(as.numeric(freq_i_red$Annee_semester_to.nb), freq_i_red$median.frsi_Nb.Total, as.numeric(freq_i_red$Annee_semester_to.nb), freq_i_red$min.frsi_Nb.Total), code = 3, angle = 90, length = 0.00) + #arrows(as.numeric(freq_i_red$Annee_semester_to.nb), freq_i_red$median.frsi_Nb.Total, as.numeric(freq_i_red$Annee_semester_to.nb), freq_i_red$min.frsi_Nb.Total, code = 3, angle = 90, length = 0.00) + } else { + freq_plot <- freq_plot + ggplot2::geom_point(ggplot2::aes(0, 0), col = "white") + } + + freq_plot <- freq_plot + ggplot2::geom_point(ggplot2::aes(freq_i$Annee_semester_to.nb, freq_i$mean.frsi_Nb.Total), pch = 19, col = "orange") + + ggplot2::geom_point(ggplot2::aes(freq_i$Annee_semester_to.nb, freq_i$median.frsi_Nb.Total), pch = 19, col = "red") + #points(freq_i$Annee_semester_to.nb, freq_i$mean.frsi_Nb.Total, pch = 19, col = "orange") + #points(freq_i$Annee_semester_to.nb, freq_i$median.frsi_Nb.Total, pch = 19, col = "red") + + + freq_i_red <- dplyr::filter(freq_i, nb.frsssi_Nb.Total >= 2) + + if (nrow(freq_i_red) >= 1) { + freq_plot <- freq_plot + ggplot2::geom_pointrange(ggplot2::aes(as.numeric(freq_i_red$Annee_semester_to.nb) + .25, freq_i_red$median.frsssi_Nb.Total, as.numeric(freq_i_red$Annee_semester_to.nb) + .25, freq_i_red$max.frsssi_Nb.Total), code = 3, angle = 90, length = 0.00) + #arrows(as.numeric(freq_i_red$Annee_semester_to.nb) + .25, freq_i_red$median.frsssi_Nb.Total, as.numeric(freq_i_red$Annee_semester_to.nb) + .25, freq_i_red$max.frsssi_Nb.Total, code = 3, angle = 90, length = 0.00) + } else { + freq_plot <- freq_plot + ggplot2::geom_point(ggplot2::aes(0, 0), col = "white") + #points(0, 0, col = "white") + } + + if (nrow(freq_i_red) >= 1) { + freq_plot <- freq_plot + ggplot2::geom_pointrange(ggplot2::aes(as.numeric(freq_i_red$Annee_semester_to.nb) + .25, freq_i_red$median.frsssi_Nb.Total, as.numeric(freq_i_red$Annee_semester_to.nb) + .25, freq_i_red$min.frsssi_Nb.Total), code = 3, angle = 90, length = 0.00) + #arrows(as.numeric(freq_i_red$Annee_semester_to.nb) + .25, freq_i_red$median.frsssi_Nb.Total, as.numeric(freq_i_red$Annee_semester_to.nb) + .25, freq_i_red$min.frsssi_Nb.Total, code = 3, angle = 90, length = 0.00) + } else { + freq_plot <- freq_plot + ggplot2::geom_point(ggplot2::aes(0, 0), col = "white") + #points(0, 0, col = "white") + } + + freq_plot <- freq_plot + ggplot2::geom_point(ggplot2::aes(as.numeric(freq_i$Annee_semester_to.nb) + .25, freq_i$mean.frsssi_Nb.Total), pch = 19, col = "lightblue") + + ggplot2::geom_point(ggplot2::aes(as.numeric(freq_i$Annee_semester_to.nb) + .25, freq_i$median.frsssi_Nb.Total), pch = 19, col = "darkblue") + #points(as.numeric(freq_i$Annee_semester_to.nb) + .25, freq_i$mean.frsssi_Nb.Total, pch = 19, col = "lightblue") + #points(as.numeric(freq_i$Annee_semester_to.nb) + .25, freq_i$median.frsssi_Nb.Total, pch = 19, col = "darkblue") + + + freq_i_red <- dplyr::filter(freq_i, nb.frha_Nb.Total >= 2) + + if (nrow(freq_i_red) >= 1) { + freq_plot <- freq_plot + ggplot2::geom_pointrange(ggplot2::aes(as.numeric(freq_i_red$Annee_semester_to.nb) - .25, freq_i_red$median.frha_Nb.Total, as.numeric(freq_i_red$Annee_semester_to.nb) - .25, freq_i_red$max.frha_Nb.Total), code = 3, angle = 90, length = 0.00) + #arrows(as.numeric(freq_i_red$Annee_semester_to.nb) - .25, freq_i_red$median.frha_Nb.Total, as.numeric(freq_i_red$Annee_semester_to.nb) - .25, freq_i_red$max.frha_Nb.Total, code = 3, angle = 90, length = 0.00) + } else { + freq_plot <- freq_plot + ggplot2::geom_point(ggplot2::aes(0, 0), col = "white") + #points(0, 0, col = "white") + } + + if (nrow(freq_i_red) >= 1) { + freq_plot <- freq_plot + ggplot2::geom_pointrange(ggplot2::aes(as.numeric(freq_i_red$Annee_semester_to.nb) - .25, freq_i_red$median.frha_Nb.Total, as.numeric(freq_i_red$Annee_semester_to.nb) - .25, freq_i_red$min.frha_Nb.Total), code = 3, angle = 90, length = 0.00) + #arrows(as.numeric(freq_i_red$Annee_semester_to.nb) - .25, freq_i_red$median.frha_Nb.Total, as.numeric(freq_i_red$Annee_semester_to.nb) - .25, freq_i_red$min.frha_Nb.Total, code = 3, angle = 90, length = 0.00) + } else { + freq_plot <- freq_plot + ggplot2::geom_point(ggplot2::aes(0, 0), col = "white") + #points(0, 0, col = "white") + } + + freq_plot <- freq_plot + ggplot2::geom_point(ggplot2::aes(as.numeric(freq_i$Annee_semester_to.nb) - .25, freq_i$mean.frha_Nb.Total), pch = 19, col = "lightgreen") + + ggplot2::geom_point(ggplot2::aes(as.numeric(freq_i$Annee_semester_to.nb) - .25, freq_i$median.frha_Nb.Total), pch = 19, col = "forestgreen") + #points(as.numeric(freq_i$Annee_semester_to.nb) - .25, freq_i$mean.frha_Nb.Total, pch = 19, col = "lightgreen") + #points(as.numeric(freq_i$Annee_semester_to.nb) - .25, freq_i$median.frha_Nb.Total, pch = 19, col = "forestgreen") + + + #legend("bottom", + #inset = c(0,-0.75), + # inset = c(0, -0.5), + # legend = c("méd. Si", "moy. Si","méd. SsSi", "moy. SsSi", "méd. Ha", "moy. Ha"), pch = c(19, 19, 19, 19, 19, 19), col = c("red", "orange", "darkblue", "lightblue", "forestgreen", "lightgreen")#, horiz = TRUE + # , bty = "n", xpd = TRUE, + # ncol = 3) + #inset = c(-0.45, 0) # You will need to fine-tune the first value depending on the windows size + # xpd = TRUE # You need to specify this graphical parameter to put the legend outside the plot + ggplot2::ggsave(paste0("freq_2_", unique(freq_i$Site), ".png"), freq_plot, height = 3, width = 3.5) +dev.off() +} + + + + + + + +for (i in c(1:length(unique(obs_$Site)))) { + + #i <- 5 + + obs_i <- dplyr::filter(obs_, Site == unique(obs_$Site)[i]) + + png(paste0("obs_2_", unique(obs_i$Site), ".png")) + plot(obs_i$date.sortie, obs_i$ratio_norm15min.obs * 100, xlim = c(xmin_, xmax_), ylim = c(0, 105), xlab = "", ylab = "comportement relatif (normalisé pour 15 min d'obs_)", main = unique(obs_i$Site), col = "black") + dev.off() + +} + diff --git a/functions.r b/functions.r new file mode 100644 index 0000000..09d5d6c --- /dev/null +++ b/functions.r @@ -0,0 +1,79 @@ +#Rscript + +######################### +## Functions ## +######################### + +#####Packages : raster +# sp +# ggplot2 + +####Set paramaters for all tools using BiodivMapR + +# path for the Mask raster corresponding to image to process +# expected to be in ENVI HDR format, 1 band, integer 8bits +# expected values in the raster: 0 = masked, 1 = selected +# set to FALSE if no mask available +input_mask_file <- FALSE + +# relative or absolute path for the Directory where results will be stored +# For each image processed, a subdirectory will be created after its name +output_dir <- "RESULTS" + +# SPATIAL RESOLUTION +# resolution of spatial units for alpha and beta diversity maps (in pixels), relative to original image +# if Res.Map = 10 for images with 10 m spatial resolution, then spatial units will be 10 pixels x 10m = 100m x 100m surfaces +# rule of thumb: spatial units between 0.25 and 4 ha usually match with ground data +# too small window_size results in low number of pixels per spatial unit, hence limited range of variation of diversity in the image +window_size <- 10 + +# PCA FILTERING: Set to TRUE if you want second filtering based on PCA outliers to be processed. Slower +filterpca <- TRUE + +################################################################################ +## DEFINE PARAMETERS FOR METHOD ## +################################################################################ +nbcpu <- 4 +maxram <- 0.5 +nbclusters <- 50 + +################################################################################ +## PROCESS IMAGE ## +################################################################################ +# 1- Filter data in order to discard non vegetated / shaded / cloudy pixels +ndvi_thresh <- 0.5 +blue_thresh <- 500 +nir_thresh <- 1500 +continuum_removal <- TRUE + + + +#### Convert raster to dataframe + +# Convert raster to SpatialPointsDataFrame +convert_raster <- function(data_raster) { +r_pts <- raster::rasterToPoints(data_raster, spatial = TRUE) + +# reproject sp object +geo_prj <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" +r_pts <- sp::spTransform(r_pts, sp::CRS(geo_prj)) + + +# Assign coordinates to @data slot, display first 6 rows of data.frame +r_pts@data <- data.frame(r_pts@data, longitude = sp::coordinates(r_pts)[, 1], + latitude = sp::coordinates(r_pts)[, 2]) + +return(r_pts@data) +} + + +#### Potting indices + +plot_indices <- function(data, titre) { + graph_indices <- ggplot2::ggplot(data) + + ggplot2::geom_point(ggplot2::aes_string(x = data[, 2], y = data[, 3], color = data[, titre]), shape = "square", size = 2) + ggplot2::scale_colour_gradient(low = "blue", high = "orange", na.value = "grey50") + + ggplot2::xlab("Longitude") + ggplot2::ylab("Latitude") + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1), plot.title = ggplot2::element_text(color = "black", size = 12, face = "bold")) + ggplot2::ggtitle(titre) + +ggplot2::ggsave(paste0(titre, ".png"), graph_indices, width = 12, height = 10, units = "cm") +} diff --git a/intersite.r b/intersite.r index 09bcabf..326bac3 100644 --- a/intersite.r +++ b/intersite.r @@ -156,7 +156,7 @@ p_freq <- ggplot2::ggplot(freq_new, ggplot2::aes(x = Site, y = Fr_Nb)) + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1)) + ggplot2::geom_boxplot(data = freq_min5, ggplot2::aes(x = Site, y = Fr_Nb)) + ggplot2::geom_point(data = freq_inf5, ggplot2::aes(x = Site, y = Fr_Nb)) -ggplot2::ggsave("5_Frequentation_2021.png", p_freq) +ggplot2::ggsave("5_Frequentation.png", p_freq) # plot frequentation data with missing site @@ -198,11 +198,11 @@ p_freq <- ggplot2::ggplot(freq_new, ggplot2::aes(x = Site, y = Fr_Nb)) + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1)) + ggplot2::geom_boxplot(data = freq_min5, ggplot2::aes(x = Site, y = Fr_Nb)) + ggplot2::geom_point(data = freq_inf5, ggplot2::aes(x = Site, y = Fr_Nb)) -ggplot2::ggsave("5_Frequentation_na_2021.png", p_freq) +ggplot2::ggsave("5_Frequentation_na.png", p_freq) freq_new_stats_site_fr_nb <- na.omit(freq_new[, c("Site", "Fr_Nb")]) %>% dplyr::group_by(Site) %>% dplyr::summarise(min. = min(Fr_Nb), max. = max(Fr_Nb), mean. = mean(Fr_Nb), median. = median(Fr_Nb), nb = dplyr::n()) freq_new_stats_site_fr_nb %>% print(n = nrow(freq_new_stats_site_fr_nb)) -saveRDS(freq_new_stats_site_fr_nb, "freq_new_stats_site_fr_nb.RDS") +#saveRDS(freq_new_stats_site_fr_nb, "freq_new_stats_site_fr_nb.RDS") freq_new_stats_effort <- na.omit(freq_new[, c("Site", "Annee_semester")]) %>% dplyr::group_by(Site, Annee_semester) %>% dplyr::summarize(nb. = dplyr::n()) freq_new_stats_effort <- freq_new_stats_effort %>% tidyr::spread(Annee_semester, nb.) @@ -218,16 +218,16 @@ freq_new <- dplyr::filter(freq_new, Site %notin% c("FINS_Quemenes", "FINS_SeinGo stat_obs_ <- obs_ %>% dplyr::group_by(Site) %>% dplyr::tally() -obs_ <- dplyr::filter(obs_, Nb.Blocs.Retournes.remis.norm15min >= 0) +obs_ <- dplyr::filter(obs_, all(obs_$Nb.Blocs.Retournes.remis.norm15min) >= 0) -p_obs_cpt <- ggplot2::ggplot(obs_, ggplot2::aes(x = Site, y = ratio.norm15min.obs)) + +p_obs_cpt <- ggplot2::ggplot(obs_, ggplot2::aes(x = Site, y = ratio_norm15min.obs)) + ggplot2::geom_boxplot() + ggplot2::xlab("") + ggplot2::ylab("Comportement des pêcheurs") + ggplot2::labs(subtitle = "1 : Bad behavior don't put the block back \n0 : Good behavior put the block back") + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1)) -ggplot2::ggsave("6_Comportement_pecheurs_2021.png", p_obs_cpt) +ggplot2::ggsave("6_Comportement_pecheurs.png", p_obs_cpt) # plot comportment data with missing site @@ -248,25 +248,25 @@ obs_ <- dplyr::left_join(obs_, obs_sitebis, by = "Site") obs_ <- dplyr::arrange(obs_, Site) obs_$Site <- as.factor(obs_$Site) -obs_$Site <- factor(obs_$Site, levels = c("GONB_IlotStMichel", "ARMO_Piegu", "ARMO_Verdelet", "ARMO_Bilfot", "ARMO_IlePlate", "PDMO_Perharidy", "FINS_Quemenes", "BRES_Keraliou", "FINS_SeinGoulenez", "FINS_SeinKilaourou", "FINS_Mousterlin", "FINS_StNicolasGlenan", "GDMO_Locmariaquer", "GDMO_BegLann", "FOUR_PlateauFour", "EGMP_GroinCou", "EGMP_PasEmsembert", "EGMP_PerreAntiochat", "EGMP_Chassiron", "EGMP_BreeBains", "BASQ_FlotsBleusZP", "BASQ_FlotsBleusZF")) +obs_$Site <- factor(obs_$Site) stat_obs_ <- obs_ %>% dplyr::group_by(Site) %>% dplyr::tally() dplyr::arrange(stat_obs_, n) obs_min5 <- dplyr::filter(obs_, Site %notin% c(as.vector(unlist((dplyr::filter(stat_obs_, n < 5))["Site"])))) obs_inf5 <- dplyr::filter(obs_, Site %notin% c(as.vector(unlist((dplyr::filter(stat_obs_, n >= 5))["Site"])))) -p_obs_cpt <- ggplot2::ggplot(obs_, ggplot2::aes(x = Site, y = ratio.norm15min.obs)) + +p_obs_cpt <- ggplot2::ggplot(obs_, ggplot2::aes(x = Site, y = ratio_norm15min.obs)) + ggplot2::geom_blank() + #ggplot2::geom_jitter(shape = 16, position = position_jitter(0.2), color = "gray") + ggplot2::xlab("") + ggplot2::ylab("comportement (échelle relative)") + ggplot2::labs(subtitle = "1 : Bad behavior don't put the block back \n0 : Good behavior put the block back") + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1)) + - ggplot2::geom_boxplot(data = obs_min5, ggplot2::aes(x = Site, y = ratio.norm15min.obs)) + - ggplot2::geom_point(data = obs_inf5, ggplot2::aes(x = Site, y = ratio.norm15min.obs)) -ggplot2::ggsave("6_Comportement_pecheurs_na_2021.png", p_obs_cpt) + ggplot2::geom_boxplot(data = obs_min5, ggplot2::aes(x = Site, y = ratio_norm15min.obs)) + + ggplot2::geom_point(data = obs_inf5, ggplot2::aes(x = Site, y = ratio_norm15min.obs)) +ggplot2::ggsave("6_Comportement_pecheurs_na.png", p_obs_cpt) -obs_stats_site_ratio <- na.omit(obs_[, c("Site", "ratio.norm15min.obs")]) %>% dplyr::group_by(Site) %>% dplyr::summarise(min. = min(ratio.norm15min.obs), max. = max(ratio.norm15min.obs), mean. = mean(ratio.norm15min.obs), median. = median(ratio.norm15min.obs), nb = dplyr::n()) +obs_stats_site_ratio <- na.omit(obs_[, c("Site", "ratio_norm15min.obs")]) %>% dplyr::group_by(Site) %>% dplyr::summarise(min. = min(obs_$ratio_norm15min.obs), max. = max(obs_$ratio_norm15min.obs), mean. = mean(obs_$ratio_norm15min.obs), median. = median(obs_$ratio_norm15min.obs), nb = dplyr::n()) saveRDS(obs_stats_site_ratio, "obs_stats_site_ratio.RDS") obs_ <- tidyr::separate(obs_, date.sortie, into = c("Year", "Month", "Day"), remove = FALSE) diff --git a/intersite.xml b/intersite.xml index 353a2e9..070832c 100644 --- a/intersite.xml +++ b/intersite.xml @@ -8,12 +8,12 @@ r-readxl - - - - + + + + + - @@ -25,18 +25,19 @@ '$input_dissim' '$input_div' '$input_freq' - '$__tool_directory__/obs.RDS' + '$input_obs' '$__tool_directory__/list_date_mean.RDS' '$output_table' '$plots' ]]> - + - + + @@ -46,13 +47,14 @@ - - - - + + + + + - + @@ -66,7 +68,7 @@ Intersites comparison **What it does** -Compare indicators between sites. ⚠️ Warning ! The frequentation, the fishermen behavior and the oceano graphs are plotted for data stopping on april 2021 ! ⚠️ +Compare indicators between sites. ⚠️ Warning ! The oceano graphs are plotted for data stopping on april 2021 ! ⚠️ **Input description** diff --git a/macrobis.xml b/macrobis.xml new file mode 100644 index 0000000..c359550 --- /dev/null +++ b/macrobis.xml @@ -0,0 +1,101 @@ + + 0.0.1 + + + r-base + r-r.utils + r-rgdal + r-rgeos + r-stars + r-stringr + r-biodivmapr + r-xfun + r-rastervis + r-ggplot2 + r-mapview + r-gridextra + r-emstreer + r-filesstrings + + + + + + + + + + + + @Manual{, + title = {obisindicators: Develop marine biodiversity indicators from OBIS data}, + author = {Ben Ben and Pieter Provoost and Tylar Murray}, + year = {2022}, + note = {https://marinebon.org/obisindicators, + https://github.com/marinebon/obisindicators}, + } + + + + + + + @Manual{, + title = {preprocS2: preprocessing of Sentinel-2 data downloaded from various data hubs / produced from various atmospheric correction methods}, + author = {Jean-Baptiste Feret}, + year = {2022}, + note = {R package version 1.1.3}, + url = {https://gitlab.com/jbferet/preprocS2}, + } + + + + + + + @Manual{, + title = {biodivMapR: biodivMapR: an R package for a- and ß-diversity mapping using remotely-sensed images}, + author = {Jean-Baptiste Feret and Florian {de Boissieu}}, + year = {2022}, + note = {R package version 1.9.4}, + url = {https://github.com/jbferet/biodivMapR}, + } + + + + + + doi:10.1016/j.ecolind.2016.07.039 + doi:10.1101/2021.01.23.427872 + + + + + + @Manual{, + title = {prosail: PROSAIL leaf and canopy radiative transfer model and inversion routines}, + author = {Jean-Baptiste Feret and Florian {de Boissieu}}, + year = {2022}, + note = {R package version 1.1.1}, + url = {https://gitlab.com/jbferet/prosail}, + } + + + + + + doi.org/10.5281/zenodo.5907920 + + + + + doi.org/10.1016/j.rse.2013.09.022 + + + + + topic_0610 + topic_3050 + + +