From 8248b35ceee52bd5716f093759145b6e32c977cb Mon Sep 17 00:00:00 2001
From: "louis.heraut" <louis.heraut@inrae.fr>
Date: Fri, 21 Jan 2022 16:40:39 +0100
Subject: [PATCH] INI variable

---
 processing/analyse.R | 50 +++++++++++++++++++++++++++++++++---
 processing/format.R  | 61 +++++++++++---------------------------------
 2 files changed, 62 insertions(+), 49 deletions(-)

diff --git a/processing/analyse.R b/processing/analyse.R
index f5409f7..53691f0 100644
--- a/processing/analyse.R
+++ b/processing/analyse.R
@@ -273,6 +273,12 @@ get_VCN10trend = function (df_data, df_meta, period, p_thresold) {
 }
 
 ### 1.4. tINI date
+
+which_underfirst = function (L, UpLim) {
+    id = which(L < UpLim)[1]
+    return (id)
+}
+
 get_tINItrend = function (df_data, df_meta, period, p_thresold) {
 
     # Get all different stations code
@@ -308,6 +314,7 @@ get_tINItrend = function (df_data, df_meta, period, p_thresold) {
     for (per in period) {
 
         df_tINIEx = tibble()
+        df_tINIlist = list(data=tibble(), info=tibble())
         
         # For all the code
         for (k in 1:nCode) {
@@ -317,27 +324,64 @@ get_tINItrend = function (df_data, df_meta, period, p_thresold) {
             # Get the data associated to the code
             df_data_roll_code = df_data_roll[df_data_roll$code == code,]
 
+
+            # Get the data associated to the code
+            df_data_code = df_data[df_data$code == code,]
+
+            # Prepare the data to fit the entry of extract.Var
+            df_QNAlist_code = prepare(df_data_code,
+                                      colnamegroup=c('code'))
+            
+            # Compute the yearly mean over the data
+            df_QNAEx_code = extract.Var(data.station=df_QNAlist_code,
+                                        funct=min,
+                                        timestep='year',
+                                        period=per,
+                                        pos.datetime=1,
+                                        na.rm=TRUE)
+
+            print(code)
+            
+            print(df_QNAEx_code)
+            QNAmax = max(df_QNAEx_code$values, na.rm=TRUE)
+            print(QNAmax)
+
+
             # Prepare the data to fit the entry of extract.Var
             df_tINIlist_code = prepare(df_data_roll_code,
                                        colnamegroup=c('code'))
 
             per.start = df_meta$start_year[df_meta$code == code]
             per.start = paste(sprintf("%02d", per.start), '-01', sep='')
+
+            
+            print(per.start)
             
             # Compute the yearly min over the averaged data
             df_tINIEx_code = extract.Var(data.station=df_tINIlist_code,
-                                         funct=which.min,
+                                         funct=which_underfirst,
                                          period=per,
                                          per.start=per.start,
                                          timestep='year',
-                                         pos.datetime=1)
-            #######
+                                         pos.datetime=1,
+                                         UpLim=QNAmax)
+
             df_tINIEx_code$group1 = k
+            df_tINIlist_code$data$group = k
+            df_tINIlist_code$info$group = k
             
             # Store the results
             df_tINIEx = bind_rows(df_tINIEx, df_tINIEx_code)
             
+            df_tINIlist$data = bind_rows(df_tINIlist$data,
+                                         df_tINIlist_code$data)
+            df_tINIlist$info = bind_rows(df_tINIlist$info,
+                                         df_tINIlist_code$info)
+            
         }
+
+        # print(df_tINIEx)
+        # print(df_tINIlist)
         
         # Converts index of the tINI to the julian date associated
         df_tINIEx = prepare_date(df_tINIEx, df_tINIlist)
diff --git a/processing/format.R b/processing/format.R
index d6398ec..a7a03db 100644
--- a/processing/format.R
+++ b/processing/format.R
@@ -112,40 +112,6 @@ 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
-
     dateStart_group = summarise(group_by(df_Xlist$data, group),
                                 Date = min(Date))
     # filter(group_by(df_Xlist$data, group), Date == min(Date))
@@ -173,18 +139,16 @@ prepare_date = function(df_XEx, df_Xlist) {
         OkXEx_code = df_XEx$group1 == group        
         XEx_code = df_XEx$values[OkXEx_code]
         
-        # dXExUp = abs(XEx_code - c(XEx_code[2:length(XEx_code)],
-                                  # XEx_code[length(XEx_code)]))
-        # dXExDown = abs(XEx_code - c(XEx_code[1],
-                                    # XEx_code[1:(length(XEx_code) - 1)]))
-        # dExMean =  apply(rbind(dXExUp, dXExDown), 2, mean, na.rm=TRUE)
-
         meanXEx_code = mean(XEx_code, na.rm=TRUE)
         dXEx_code = meanXEx_code - XEx_code
         stdXEx_code = sd(XEx_code, na.rm=TRUE)
+
+        OkOverStd = dXEx_code >= stdXEx_code*3
+        OkOverStd[is.na(OkOverStd)] = FALSE
         
-        XEx_code[dXEx_code >= stdXEx_code*3] =
-            XEx_code[dXEx_code >= stdXEx_code*3] + 365
+        XEx_code[OkOverStd] = XEx_code[OkOverStd] + 365
+
+        print(OkOverStd)
 
         # print(group)
         # print(df_XEx$datetime[df_XEx$group1 == group][dXEx_code >= stdXEx_code*3])
@@ -306,10 +270,15 @@ get_period = function (per, df_Xtrend, df_XEx, df_Xlist) {
         id = which(df_Xtrend$group1 == g)
 
         # Computes index of the nearest accessible start and end date
-        iStart = which.min(abs(df_XExtmp_code$Date
-                               - as.Date(per[1])))
-        iEnd = which.min(abs(df_XExtmp_code$Date 
-                             - as.Date(per[2])))
+        OkStart = df_XExtmp_code$Date >= as.Date(per[1])
+        OkEnd = df_XExtmp_code$Date <= as.Date(per[2])
+        
+        distStart = abs(df_XExtmp_code$Date[OkStart] - as.Date(per[1]))
+        distEnd = abs(df_XExtmp_code$Date[OkEnd] - as.Date(per[2]))
+        
+        iStart = which.min(distStart)
+        iEnd = which.min(distEnd)
+
 
         # Stores the start and end of the trend analysis
         df_Xtrend$period_start[id] =
-- 
GitLab