Mercurial > repos > siwaa > redelac_stics_a

comparison redelacComputeAclim.R @ 0:80842c5ecb58 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://forgemia.inra.fr/redelac/redelac-toolbox/-/tree/68f4e9db390070700af804f77468145abfc26bc0/tools/REDELACcomputeAclim commit 68f4e9db390070700af804f77468145abfc26bc0-dirty"
author siwaa
date Tue, 24 Feb 2026 08:53:06 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:80842c5ecb58
1 library(data.table) #version 1.14.2
2 library(dplyr) #version 1.0.10
3 library(ggplot2)
4 library(rlang)
5 library(stringr)
6 sessionInfo()
7
8 args = commandArgs(trailingOnly = TRUE)
9
10 climatPath = args[1]
11 horizons = args[2]
12 rdsOut = args[3]
13 csvOut = args[4]
14 png1Out = args[5]
15 png2Out = args[6]
16
17 liste_mailles_SAFRAN <- list.dirs(climatPath,full.names = F,recursive = F)
18
19 #agregation des fichiers climatiques dans un unique fichier
20 fichier_climat_tot<- NULL
21
22
23 for(maille_DRIAS in liste_mailles_SAFRAN){
24
25 setwd(paste0(climatPath,"/",maille_DRIAS))
26
27 filenames <- list.files()
28
29
30 fichier_climat_maille_tot <- NULL
31 f=1
32 for (filename in filenames)
33 {
34 file <- as.data.frame(fread(filename,sep="\t", dec=".", header=FALSE))
35 if (f==1) {fichier_climat_maille_tot <-file} else {fichier_climat_maille_tot <- bind_rows(fichier_climat_maille_tot,file)}
36 f=f+1
37 }
38 fichier_climat_tot<- bind_rows(fichier_climat_tot,fichier_climat_maille_tot)
39 setwd("../..")
40 }
41
42 #rm(file,fichier_climat_maille_tot)
43
44 #entete des colonnes
45 # column 1: name of weather file
46 # column 2: year
47 # column 3: month
48 # column 4: day in month
49 # column 5: Julian day
50 # column 6: minimum temperature (degrees C)
51 # column 7: maximum temperature (degrees C)
52 # column 8: global radiation (MJ m-2 j-1)
53 # column 9: Penman PET (mm j-1)
54 # column 10: rainfall (mm j-1)
55 # column 11: wind (m s-1)
56 # column 12: vapour pressure (mbars)
57 # column 13: CO2 content(ppm)
58
59 names(fichier_climat_tot)
60
61 names(fichier_climat_tot)<- c("maille_DRIAS","annee","mois","jcal","jjulien","TN","TX","RG","ETP","P","V","PV","CO2")
62
63
64 # sanityzing
65
66 trimedHorizon <-str_trim(horizons)
67 purgedHorizon <-str_sub(trimedHorizon, 1, nchar(trimedHorizon)-1)
68
69 listHory <- lapply(unlist(strsplit(purgedHorizon, ",")), function(x) {tmpdata <- str_split(str_trim(x), " ")[[1]]; paste0("(annee<", tmpdata[3], ")&(annee>",tmpdata[2], ")~\"", tmpdata[1], "\"" )})
70
71 listHory
72
73 paste(listHory, collapse = ",", sep = "")
74
75 horizonsFilter <- parse_expr(paste(listHory, collapse = ",", sep = ""))
76
77 listHory
78
79 fichier_climat_tot <- fichier_climat_tot %>%
80 mutate(horizon=case_when(!!horizonsFilter))
81
82 #fichier_climat_tot <- fichier_climat_tot %>%
83 # mutate(horizon=case_when((annee>1990)&(annee<2025)~"ref"))
84
85 #ecriture du fichier
86 saveRDS(fichier_climat_tot, rdsOut)
87
88 fichier_aclim_an <- fichier_climat_tot %>%
89 select(maille_DRIAS,horizon,annee,V)%>%
90 group_by(maille_DRIAS,horizon,annee)%>%
91 summarise(Vmed = median(V,na.rm=T))%>%
92 mutate(aclim=-6.2*log(Vmed)+17.5)
93
94 fichier_aclim_hor <- fichier_aclim_an %>%
95 select(maille_DRIAS,horizon,annee,aclim)%>%
96 group_by(maille_DRIAS,horizon)%>%
97 summarise(aclim = median(aclim,na.rm=T))%>%
98 filter(is.na(horizon)==F)
99
100
101 #ecriture des valeurs de aclim par maille SAFRAN et horizon
102 write.table(fichier_aclim_hor,csvOut,sep=";",dec=",",row.names = F)
103
104
105 #tracage des valeurs calculees de aclim
106 # ggplot()+geom_boxplot(data=fichier_aclim_an,aes(y=aclim))
107 # ggplot()+geom_boxplot(data=fichier_aclim_an,aes(y=aclim,x=horizon))
108 # ggplot()+geom_boxplot(data=fichier_aclim_an,aes(y=aclim,x=as.factor(maille_DRIAS)))
109
110 ggplot()+geom_line(data=fichier_aclim_an ,aes(x=annee,y=aclim,group=maille_DRIAS,color=as.factor(maille_DRIAS)))
111 ggsave(png1Out)
112 ggplot()+
113 geom_boxplot(data=fichier_aclim_an %>% filter(is.na(horizon)==F),aes(x=as.factor(maille_DRIAS),y=aclim,fill=horizon))+
114 xlab("Maille DRIAS")+
115 theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
116 ggsave(png2Out)