Mercurial > repos > ecology > pampa_communitymetrics

comparison FunctExeCalcCommIndexesGalaxy.r @ 0:ddd5b2e74b8b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/ColineRoyaux/PAMPA-Galaxy commit 07f1028cc764f920b1e6419c151f04ab4e3600fa"
author ecology
date Tue, 21 Jul 2020 06:00:10 -0400
parents
children 5bd7ddd7601f
comparison
equal deleted inserted replaced
-1:000000000000 0:ddd5b2e74b8b
1 #Rscript
2
3 #####################################################################################################################
4 #####################################################################################################################
5 ################################# Calculate community indexes from observation data #################################
6 #####################################################################################################################
7 #####################################################################################################################
8
9 ###################### Packages R
10
11 suppressMessages(library(tidyr))
12
13 ###################### Load arguments and declaring variables
14
15 args = commandArgs(trailingOnly=TRUE)
16 #options(encoding = "UTF-8")
17
18 if (length(args) < 4) {
19 stop("At least one argument must be supplied, an input dataset file (.tabular).", call.=FALSE) #si pas d'arguments -> affiche erreur et quitte / if no args -> error and exit1
20
21 } else {
22 Importdata<-args[1] ###### Nom du fichier importé avec son extension / file name imported with the file type ".filetype"
23 index <- args[2] ###### List of selected metrics to calculate
24 source(args[3]) ###### Import functions
25
26 }
27 #### Data must be a dataframe with at least 3 variables : unitobs representing location and year ("observation.unit"), species code ("species.code") and abundance ("number")
28
29
30 #Import des données / Import data
31 obs<- read.table(Importdata,sep="\t",dec=".",header=TRUE,encoding="UTF-8") #
32 obs[obs == -999] <- NA
33 factors <- fact.det.f(Obs=obs)
34 ObsType <- def.typeobs.f(Obs=obs)
35 obs <- create.unitobs(data=obs)
36
37 vars_data<-c("observation.unit","species.code","number")
38 err_msg_data<-"The input dataset doesn't have the right format. It need to have at least the following 3 variables :\n- observation.unit (or point and year)\n- species.code\n- number\n"
39 check_file(obs,err_msg_data,vars_data,3)
40
41
42
43 ####################################################################################################
44 ################## create community metrics table ## Function : calcBiodiv.f #######################
45 ####################################################################################################
46
47 ########################################################################################################################
48 calcBiodiv.f <- function(Data, MPA, unitobs="observation.unit", code.especes="species.code", nombres="number",
49 indices=index)
50 {
51 ## Purpose: calcul des indices de biodiversité
52 ## ----------------------------------------------------------------------
53 ## Arguments: Data : les données à partir desquelles calculer les
54 ## indices. Doivent comporter au minimum (colones) :
55 ## * unités d'observations/sites
56 ## * espèces présentes
57 ## * nombre d'individus /espèce/unitobs.
58 ## refesp : le référentiel espèces.
59 ## MPA : l'AMP (chaîne de charactères).
60 ## unitobs : nom de la colone d'unités d'observation.
61 ## code.especes : nom de la colone d'espèces.
62 ## nombres : nom de la colone de nombres.
63 ## indices : liste des indices à calculer
64 ## (vecteur de caractères)
65 ## ----------------------------------------------------------------------
66 ## Author: Yves Reecht, Date: 29 oct. 2010, 08:58
67
68 ## Supression de tout ce qui n'a pas d'espèce précisee (peut être du non biotique ou identification >= genre) :
69
70 notspline <- grep("(sp\\.)$|([1-9])$|^(Absencemacrofaune)$|^(NoID)$|^(Acrobranc)$|^(Acrodigit)$|^(Acroencr)$|^(Acrosubm)$|^(Acrotabu)$|^(Adredure)$|^(Adremoll)$|^(Algaturf)$|^(Balimona)$|^(Corablan)$|^(CoradurV)$|^(Coraenal)$|^(Coramor1)$|^(Coramor2)$|^(Coramou)$|^( Dallcora)$|^(Debrcora)$|^(Debris)$|^(Hare)$|^(HexaChar)$|^(MuraCong)$|^(Nacrbran)$|^(Nacrcham)$|^(Nacrencr)$|^(Nacrfoli)$|^(Nacrmass)$|^(Nacrsubm)$|^(Recrcora)$|^(Roche)$|^(Sable)$|^(Vase)$",Data[, code.especes], value=FALSE)
71 if (length(notspline) != 0)
72 {
73 Data <- Data[-notspline, ]
74 }else{}
75
76 ## Suppression des niveaux de facteur inutilisés :
77 Data <- dropLevels.f(df=Data)
78
79
80 ## Si les données ne sont pas encore agrégées /espèce/unitobs on le fait ici :
81 if (nrow(Data) > nrow(expand.grid(unique(Data[ , unitobs]), unique(Data[ , code.especes]))))
82 {
83 Data <- agregations.generic.f(Data=Data, metrics=nombres,
84 factors=c(unitobs, code.especes),
85 listFact=NULL)
86 }else{}
87
88 df.biodiv <- as.data.frame(as.table(tapply(Data[ , nombres],
89 Data[ , unitobs],
90 sum, na.rm=TRUE)))
91
92 colnames(df.biodiv) <- c(unitobs, nombres)
93
94 ## ##################################################
95 ## Richesse spécifique :
96 Data$pres.abs <- presAbs.f(nombres=Data[ , nombres], logical = FALSE)
97
98 df.biodiv$species.richness <- as.vector(tapply(Data$pres.abs,
99 Data[ , unitobs], sum, na.rm=TRUE),
100 "integer")
101 ## ... as.vector to avoid the class "array".
102
103 ## ##################################################
104 ## Indices de Simpson et Shannon et dérivés :
105
106 matNombres <- tapply(Data[ , nombres], # Matrice de nombres d'individus /espèce/unitobs.
107 list(Data[ , unitobs], Data[ , code.especes]),
108 sum, na.rm=TRUE)
109
110 matNombres[is.na(matNombres)] <- 0 # Vrais zéros
111
112 ## Proportion d'individus de chaque espèce dans l'unitobs :
113 propIndiv <- sweep(matNombres, 1,
114 apply(matNombres, 1, sum, na.rm = TRUE), # Nombre d'individus / unitobs ; équiv df.biodiv$nombre.
115 FUN="/")
116
117 ## Indices de Simpson :
118 df.biodiv$simpson <- apply(propIndiv^2, 1, sum, na.rm=TRUE)
119
120 if (any(is.element(c("all", "simpson.l"), indices)))
121 {
122 df.biodiv$simpson.l <- 1 - df.biodiv$simpson
123 }
124
125 ## calcul de l'indice de Shannon :
126 df.biodiv$shannon <- -1 * apply(propIndiv * log(propIndiv), 1, sum, na.rm=TRUE)
127
128 ## calcul de l'indice de Pielou :
129 if (any(is.element(c("all", "pielou"), indices)))
130 {
131 df.biodiv$pielou <- df.biodiv$shannon / log(df.biodiv$species.richness)
132 }
133
134 ## calcul de l'indice de Hill :
135 if (any(is.element(c("all", "hill"), indices)))
136 {
137 df.biodiv$hill <- (1 - df.biodiv$simpson) / exp(df.biodiv$shannon)
138 # équiv df.biodiv$l.simpson / exp(df.biodiv$shannon)
139 }
140
141
142 return(df.biodiv)
143 }
144
145 ################# Analysis
146
147 res <- calc.numbers.f(obs, ObsType=ObsType , factors=factors, nbName="number")
148
149 tableCommunityIndexes <- calcBiodiv.f(res, MPA, unitobs="observation.unit", code.especes="species.code", nombres="number",
150 indices=index)
151 tableCommunityIndexes <- create.year.point(tableCommunityIndexes)
152 #Save dataframe in a tabular format
153
154 filenameComm <- "TabCommunityIndexes.tabular"
155 write.table(tableCommunityIndexes, filenameComm, row.names=FALSE, sep="\t", dec=".",fileEncoding="UTF-8")
156 cat(paste("\nWrite table with Community indexes. \n--> \"",filenameComm,"\"\n",sep=""))
157