From c0b45a867014b10cb880b0cf50f99c10a8cd4cc7 Mon Sep 17 00:00:00 2001
From: "louis.heraut" <louis.heraut@inrae.fr>
Date: Wed, 12 Jan 2022 20:31:59 +0100
Subject: [PATCH] Cleaned and commented
---
processing/analyse.R | 74 +++++++++++++++++++++++++++++++++++++++++++-
processing/format.R | 47 ++++++++++++++++++++++++++++
script.R | 25 +++++++++------
3 files changed, 136 insertions(+), 10 deletions(-)
diff --git a/processing/analyse.R b/processing/analyse.R
index 74e606c..a56cf57 100644
--- a/processing/analyse.R
+++ b/processing/analyse.R
@@ -206,7 +206,8 @@ get_QMNAtrend = function (df_data, period, p_thresold) {
}
### 1.3. VCN10
-# Realise the trend analysis of the minimum 10 day average flow over the year (VCN10) hydrological variable
+# Realises the trend analysis of the minimum 10 day average flow
+# over the year (VCN10) hydrological variable
get_VCN10trend = function (df_data, df_meta, period, p_thresold) {
# Get all different stations code
@@ -274,6 +275,77 @@ get_VCN10trend = function (df_data, df_meta, period, p_thresold) {
}
+### 1.4. VCN10 date
+# Realises the trend analysis of the date of the minimum 10 day
+# average flow over the year (VCN10) hydrological variable
+get_dateVCN10trend = function (df_data, df_meta, period, p_thresold) {
+
+ # Get all different stations code
+ Code = levels(factor(df_meta$code))
+ # Blank tibble to store the data averaged
+ df_data_roll = tibble()
+
+ # For all the code
+ for (c in Code) {
+ # Get the data associated to the code
+ df_data_code = df_data[df_data$code == c,]
+ # Perform the roll mean of the flow over 10 days
+ df_data_code = tibble(Date=df_data_code$Date,
+ Qm3s=rollmean(df_data_code$Qm3s,
+ 10,
+ fill=NA),
+ code=c)
+ # Store the results
+ df_data_roll = bind_rows(df_data_roll, df_data_code)
+ }
+
+ # Make sure to convert the period to a list
+ period = as.list(period)
+ # Set the max interval period as the minimal possible
+ Imax = 0
+ # Blank tibble for data to return
+ df_VCN10trendB = tibble()
+
+ # For all periods
+ for (per in period) {
+ # Prepare the data to fit the entry of extract.Var
+ df_VCN10list = prepare(df_data_roll, colnamegroup=c('code'))
+ # Compute the yearly min over the averaged data
+ df_VCN10Ex = extract.Var(data.station=df_VCN10list,
+ funct=which.min,
+ period=per,
+ timestep='year',
+ pos.datetime=1)
+
+ # Converts index of the VCN10 to the julian date associated
+ df_VCN10Ex = prepare_date(df_VCN10Ex, df_VCN10list)
+
+ # Compute the trend analysis
+ df_VCN10trend = Estimate.stats(data.extract=df_VCN10Ex,
+ level=p_thresold)
+
+ # Get the associated time interval
+ I = interval(per[1], per[2])
+ # If it is the largest interval
+ if (I > Imax) {
+ # Store it and the associated data and info
+ Imax = I
+ df_VCN10listB = df_VCN10list
+ df_VCN10ExB = df_VCN10Ex
+ }
+
+ # Specify the period of analyse
+ df_VCN10trend = get_period(per, df_VCN10trend, df_VCN10Ex,
+ df_VCN10list)
+ # Store the trend
+ df_VCN10trendB = bind_rows(df_VCN10trendB, df_VCN10trend)
+ }
+ # Clean results of trend analyse
+ res_VCN10trend = clean(df_VCN10trendB, df_VCN10ExB, df_VCN10listB)
+ return (res_VCN10trend)
+}
+
+
## 2. OTHER ANALYSES
### 2.1. Break date
# Compute the break date of the flow data by station
diff --git a/processing/format.R b/processing/format.R
index 2123372..7b17395 100644
--- a/processing/format.R
+++ b/processing/format.R
@@ -107,6 +107,53 @@ reprepare = function(df_XEx, df_Xlist, colnamegroup=NULL) {
}
+
+
+
+prepare_date = function(df_XEx, df_Xlist) {
+
+ dateMD = substr(df_Xlist$data$Date, 6, 10)
+ idYearList = which(dateMD == "01-01")
+ YearList = df_Xlist$data$Date[idYearList]
+ groupList = df_Xlist$data$group[idYearList]
+
+ Group = levels(factor(df_XEx$group1))
+ for (group in Group) {
+
+ Year = YearList[groupList == group]
+
+ for (idyear in 1:length(Year)) {
+ year = Year[idyear]
+
+ idShift = which(df_Xlist$data$group == group
+ & df_Xlist$data$Date == year)
+
+ okEx = df_XEx$group1 == group & df_XEx$datetime == substr(year, 1, 4)
+
+ df_XEx$values[okEx] = df_XEx$values[okEx] + idShift
+ }
+ }
+
+ idDate = df_XEx$values
+ dateEx = as.Date(df_Xlist$data$Date[idDate])
+ origin = as.Date(paste(format(dateEx, "%Y"), '-01-01', sep=''))
+ ndate = length(origin)
+
+ dateExJul = rep(0, times=ndate)
+ for (i in 1:ndate) {
+ dateExJul[i] = julian(dateEx[i], origin=origin[i])
+ }
+
+ df_XEx$values = dateExJul
+
+ return (df_XEx)
+}
+
+
+
+
+
+
## 2. OUTPUT
# Cleans the trend results of the function 'Estimate.stats' in the
# 'StatsAnalysisTrend' package. It adds the station code and the
diff --git a/script.R b/script.R
index 02e9f9e..75400e2 100644
--- a/script.R
+++ b/script.R
@@ -55,7 +55,7 @@ filedir =
# Name of the file that will be analysed from the BH directory
# (if 'all', all the file of the directory will be chosen)
filename =
- ""
+ # ""
# c(
# "S2235610_HYDRO_QJM.txt",
@@ -65,10 +65,12 @@ filename =
# "O0384010_HYDRO_QJM.txt"
# )
- # c("S4214010_HYDRO_QJM.txt",
- # "O0384010_HYDRO_QJM.txt",
+ c(
+ # "S4214010_HYDRO_QJM.txt",
+ "O0384010_HYDRO_QJM.txt",
# "Q7002910_HYDRO_QJM.txt",
- # "O0485110_HYDRO_QJM.txt")
+ "O3006710_HYDRO_QJM.txt"
+ )
@@ -79,8 +81,8 @@ AGlistdir =
""
AGlistname =
- # ""
- "Liste-station_RRSE.docx"
+ ""
+ # "Liste-station_RRSE.docx"
## NIVALE SELECTION
@@ -238,6 +240,11 @@ res_VCN10trend = get_VCN10trend(df_data, df_meta,
period=trend_period,
p_thresold=p_thresold)
+# VCN10 date trend
+res_dateVCN10trend = get_dateVCN10trend(df_data, df_meta,
+ period=trend_period,
+ p_thresold=p_thresold)
+
### 3.3. Break analysis
# df_break = get_break(res_QAtrend$data, df_meta)
# df_break = get_break(res_QMNAtrend$data, df_meta)
@@ -269,9 +276,9 @@ df_shapefile = ini_shapefile(computer_data_path, fr_shpdir, fr_shpname, bs_shpdi
### 4.2. Analysis layout
datasheet_layout(isplot=c(
- 'datasheet',
- 'matrix',
- 'map'
+ 'datasheet'
+ # 'matrix',
+ # 'map'
),
df_data=list(res_QAtrend$data, res_QMNAtrend$data,
res_VCN10trend$data),
--
GitLab