Cleaned and commented

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 

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

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),