diff --git a/processing/analyse.R b/processing/analyse.R index 35fea94202bad5e16a921e97a391183d26370962..f5409f7ae42aea3717c830046340ca4af9521753 100644 --- a/processing/analyse.R +++ b/processing/analyse.R @@ -216,19 +216,17 @@ get_VCN10trend = function (df_data, df_meta, period, p_thresold) { df_data_roll = tibble() # For all the code - for (c in Code) { + for (code in Code) { # Get the data associated to the code - df_data_code = df_data[df_data$code == c,] + df_data_code = df_data[df_data$code == code,] # Perform the roll mean of the flow over 10 days - df_data_code = tibble(Date=rollmean(df_data_code$Date, - 10, - fill=NA), - Value=rollmean(df_data_code$Value, - 10, - fill=NA), - code=c) + df_data_roll_code = tibble(Date=df_data_code$Date, + Value=rollmean(df_data_code$Value, + 10, + fill=NA), + code=code) # Store the results - df_data_roll = bind_rows(df_data_roll, df_data_code) + df_data_roll = bind_rows(df_data_roll, df_data_roll_code) } # Make sure to convert the period to a list @@ -274,77 +272,100 @@ get_VCN10trend = function (df_data, df_meta, period, p_thresold) { return (res_VCN10trend) } +### 1.4. tINI date +get_tINItrend = function (df_data, df_meta, period, p_thresold) { + # Get all different stations code + Code = levels(factor(df_meta$code)) + # Gets the number of station + nCode = length(Code) + + # Blank tibble to store the data averaged + df_data_roll = tibble() -### 1.4. tINI date -# get_tINItrend = 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, -# Value=rollmean(df_data_code$Value, -# 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_tMIDtrendB = tibble() - -# # For all periods -# for (per in period) { -# # Prepare the data to fit the entry of extract.Var -# df_tMIDlist = prepare(df_data_roll, colnamegroup=c('code')) -# # Compute the yearly min over the averaged data -# df_tMIDEx = extract.Var(data.station=df_tMIDlist, -# funct=which.min, -# period=per, -# timestep='year', -# pos.datetime=1) - -# # Converts index of the tMID to the julian date associated -# df_tMIDEx = prepare_date(df_tMIDEx, df_tMIDlist) + # For all the code + for (code in Code) { + # Get the data associated to the code + df_data_code = df_data[df_data$code == code,] + # Perform the roll mean of the flow over 10 days + df_data_roll_code = tibble(Date=df_data_code$Date, + Value=rollmean(df_data_code$Value, + 10, + fill=NA), + code=code) + # Store the results + df_data_roll = bind_rows(df_data_roll, df_data_roll_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_tINItrendB = tibble() + + # For all periods + for (per in period) { + + df_tINIEx = tibble() + + # For all the code + for (k in 1:nCode) { + # Gets the code + code = Code[k] + + # Get the data associated to the code + df_data_roll_code = df_data_roll[df_data_roll$code == code,] + + # 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='') + + # Compute the yearly min over the averaged data + df_tINIEx_code = extract.Var(data.station=df_tINIlist_code, + funct=which.min, + period=per, + per.start=per.start, + timestep='year', + pos.datetime=1) + ####### + df_tINIEx_code$group1 = k + + # Store the results + df_tINIEx = bind_rows(df_tINIEx, df_tINIEx_code) + + } + + # Converts index of the tINI to the julian date associated + df_tINIEx = prepare_date(df_tINIEx, df_tINIlist) -# # Compute the trend analysis -# df_tMIDtrend = Estimate.stats(data.extract=df_tMIDEx, -# 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_tMIDlistB = df_tMIDlist -# df_tMIDExB = df_tMIDEx -# } - -# # Specify the period of analyse -# df_tMIDtrend = get_period(per, df_tMIDtrend, df_tMIDEx, -# df_tMIDlist) -# # Store the trend -# df_tMIDtrendB = bind_rows(df_tMIDtrendB, df_tMIDtrend) -# } -# # Clean results of trend analyse -# res_tMIDtrend = clean(df_tMIDtrendB, df_tMIDExB, df_tMIDlistB) -# return (res_tMIDtrend) -# } + # Compute the trend analysis + df_tINItrend = Estimate.stats(data.extract=df_tINIEx, + 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_tINIlistB = df_tINIlist + df_tINIExB = df_tINIEx + } + # Specify the period of analyse + df_tINItrend = get_period(per, df_tINItrend, df_tINIEx, + df_tINIlist) + # Store the trend + df_tINItrendB = bind_rows(df_tINItrendB, df_tINItrend) + } + # Clean results of trend analyse + res_tINItrend = clean(df_tINItrendB, df_tINIExB, df_tINIlistB) + return (res_tINItrend) +} ### 1.5. tMID date # Realises the trend analysis of the date of the minimum 10 day @@ -357,17 +378,17 @@ get_tMIDtrend = function (df_data, df_meta, period, p_thresold) { df_data_roll = tibble() # For all the code - for (c in Code) { + for (code in Code) { # Get the data associated to the code - df_data_code = df_data[df_data$code == c,] + df_data_code = df_data[df_data$code == code,] # Perform the roll mean of the flow over 10 days - df_data_code = tibble(Date=df_data_code$Date, - Value=rollmean(df_data_code$Value, - 10, - fill=NA), - code=c) + df_data_roll_code = tibble(Date=df_data_code$Date, + Value=rollmean(df_data_code$Value, + 10, + fill=NA), + code=code) # Store the results - df_data_roll = bind_rows(df_data_roll, df_data_code) + df_data_roll = bind_rows(df_data_roll, df_data_roll_code) } # Make sure to convert the period to a list @@ -465,6 +486,8 @@ get_hydrograph = function (df_data, period=NULL, df_meta=NULL) { if (!is.null(df_meta)) { # New column in metadata for hydrological regime df_meta$regime_hydro = NA + # New column in metadata for the start of the hydrological year + df_meta$start_year = NA # Get all different stations code Code = levels(factor(df_meta$code)) @@ -530,7 +553,10 @@ get_hydrograph = function (df_data, period=NULL, df_meta=NULL) { df_QM = bind_rows(df_QM, df_QMtmp) # Stores result of the hydrological regime df_meta$regime_hydro[df_meta$code == code] = classRegime - + # Computes the month of the max QM + maxMonth = which.max(df_QM$QM) + # Stores it as the start of the hydrological year + df_meta$start_year[df_meta$code == code] = maxMonth # Otherwise } else { # No tibble needed diff --git a/script.R b/script.R index c659ff6f9167ffde6e7ba336a82009ad32361073..f478638708fdb5a5e885fad611550ecef3cd2927 100644 --- a/script.R +++ b/script.R @@ -244,9 +244,9 @@ res_VCN10trend = get_VCN10trend(df_data, df_meta, p_thresold=p_thresold) # Start date for low water trend -# res_tINItrend = get_tINItrend(df_data, df_meta, - # period=trend_period, - # p_thresold=p_thresold) +res_tINItrend = get_tINItrend(df_data, df_meta, + period=trend_period, + p_thresold=p_thresold) # Center date for low water trend res_tMIDtrend = get_tMIDtrend(df_data, df_meta, @@ -297,28 +297,28 @@ datasheet_layout(toplot=c( df_data=list(res_QAtrend$data, res_QMNAtrend$data, res_VCN10trend$data, - # res_tINItrend$data, + res_tINItrend$data, res_tMIDtrend$data), df_trend=list(res_QAtrend$trend, res_QMNAtrend$trend, res_VCN10trend$trend, - # res_tINItrend$trend, + res_tINItrend$trend, res_tMIDtrend$trend), var=list('QA', 'QMNA', 'VCN10', - # 'tINI', + 'tINI', 'tMID'), type=list('flow', 'flow', 'flow', - # 'date', + 'date', 'date'), - layout_matrix=matrix(c(1, 2, 3, 4), ncol=1), + layout_matrix=matrix(c(1, 2, 3, 4, 5), ncol=1), missRect=TRUE, trend_period=trend_period,