diff --git a/plotting/layout.R b/plotting/layout.R
index 44a84f98694cbbeeb6bb3c9952d901a4ba0c005b..b046c4f491662d353ec3a15743fb2186eb5657e0 100644
--- a/plotting/layout.R
+++ b/plotting/layout.R
@@ -127,7 +127,7 @@ gg_circle = function(r, xc, yc, color="black", fill=NA, ...) {
## 3. LAYOUT
# Generates a PDF that gather datasheets, map and summarize matrix about the trend analyses realised on selected stations
datasheet_layout = function (df_data, df_meta, layout_matrix,
- isplot=c('datasheet', 'matrix', 'map'),
+ toplot=c('datasheet', 'matrix', 'map'),
figdir='', filedir_opt='', filename_opt='',
variable='', df_trend=NULL,
p_threshold=0.1, unit2day=365.25, var='',
@@ -231,18 +231,18 @@ datasheet_layout = function (df_data, df_meta, layout_matrix,
}
# If datasheets needs to be plot
- if ('datasheet' %in% isplot) {
+ if ('datasheet' %in% toplot) {
datasheet_panel(list_df2plot, df_meta, trend_period, info_header=info_header, time_header=time_header, layout_matrix=layout_matrix, info_ratio=info_ratio, time_ratio=time_ratio, var_ratio=var_ratio, outdirTmp=outdirTmp)
}
# If summarize matrix needs to be plot
- if ('matrix' %in% isplot) {
+ if ('matrix' %in% toplot) {
matrix_panel(list_df2plot, df_meta, trend_period, mean_period,
slice=12, outdirTmp=outdirTmp, A3=TRUE)
}
# If map needs to be plot
- if ('map' %in% isplot) {
+ if ('map' %in% toplot) {
map_panel(list_df2plot,
df_meta,
idPer=length(trend_period),
diff --git a/plotting/matrix.R b/plotting/matrix.R
index f8454d5175d045ddb94a539bef6b3f8c3bf1f434..4e940c43fc1dc9e3065c6a0f858f15f32c8b0987 100644
--- a/plotting/matrix.R
+++ b/plotting/matrix.R
@@ -515,9 +515,13 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice
Color_trend_per =
subColor_trend[subNPeriod_trend == j]
+ # Converts the variable list into levels for factor
+ levels = unlist(Var_trend_per)
# Converts the vector of hydrological variable to
# a vector of integer associated to those variable
- Xtmp = as.integer(factor(as.character(Var_trend_per)))
+ Xtmp = as.integer(factor(as.character(Var_trend_per),
+ levels=levels))
+
# Computes X position of the column for the period dates
Xc = j + (j - 1)*nbp*2
# Computes X positions of columns for the mean of
@@ -691,9 +695,12 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice
Color_mean_per =
subColor_mean[subNPeriod_mean == j]
+ # Converts the variable list into levels for factor
+ levels = unlist(Var_mean_per)
# Converts the vector of hydrological variable to
# a vector of integer associated to those variable
- Xtmp_mean = as.integer(factor(as.character(Var_mean_per)))
+ Xtmp_mean = as.integer(factor(as.character(Var_mean_per),
+ levels=levels))
# Computes X position of the column for the period dates
Xc_mean = j + (j - 1)*nbp + X[length(X)]
# Computes X positions of columns for the mean of
diff --git a/processing/analyse.R b/processing/analyse.R
index 9e8eb65ba0f4292107c600aeffc4e26d389bf05e..7ec789c5c1f7585ae824110ae9f457b681f12b89 100644
--- a/processing/analyse.R
+++ b/processing/analyse.R
@@ -275,10 +275,81 @@ get_VCN10trend = function (df_data, df_meta, period, p_thresold) {
}
-### 1.4. VCN10 date
+
+### 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)
+
+# # 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)
+# }
+
+
+
+### 1.5. tMID 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_tMIDtrend = function (df_data, df_meta, period, p_thresold) {
# Get all different stations code
Code = levels(factor(df_meta$code))
@@ -304,24 +375,24 @@ get_dateVCN10trend = function (df_data, df_meta, period, p_thresold) {
# Set the max interval period as the minimal possible
Imax = 0
# Blank tibble for data to return
- df_VCN10trendB = tibble()
+ df_tMIDtrendB = 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'))
+ df_tMIDlist = prepare(df_data_roll, colnamegroup=c('code'))
# Compute the yearly min over the averaged data
- df_VCN10Ex = extract.Var(data.station=df_VCN10list,
+ df_tMIDEx = extract.Var(data.station=df_tMIDlist,
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)
+ # Converts index of the tMID to the julian date associated
+ df_tMIDEx = prepare_date(df_tMIDEx, df_tMIDlist)
# Compute the trend analysis
- df_VCN10trend = Estimate.stats(data.extract=df_VCN10Ex,
+ df_tMIDtrend = Estimate.stats(data.extract=df_tMIDEx,
level=p_thresold)
# Get the associated time interval
@@ -330,24 +401,148 @@ get_dateVCN10trend = function (df_data, df_meta, period, p_thresold) {
if (I > Imax) {
# Store it and the associated data and info
Imax = I
- df_VCN10listB = df_VCN10list
- df_VCN10ExB = df_VCN10Ex
+ df_tMIDlistB = df_tMIDlist
+ df_tMIDExB = df_tMIDEx
}
# Specify the period of analyse
- df_VCN10trend = get_period(per, df_VCN10trend, df_VCN10Ex,
- df_VCN10list)
+ df_tMIDtrend = get_period(per, df_tMIDtrend, df_tMIDEx,
+ df_tMIDlist)
# Store the trend
- df_VCN10trendB = bind_rows(df_VCN10trendB, df_VCN10trend)
+ df_tMIDtrendB = bind_rows(df_tMIDtrendB, df_tMIDtrend)
}
# Clean results of trend analyse
- res_VCN10trend = clean(df_VCN10trendB, df_VCN10ExB, df_VCN10listB)
- return (res_VCN10trend)
+ res_tMIDtrend = clean(df_tMIDtrendB, df_tMIDExB, df_tMIDlistB)
+ return (res_tMIDtrend)
}
## 2. OTHER ANALYSES
-### 2.1. Break date
+### 2.1. Hydrograph
+xref = matrix(
+ c(0.033, 0.025, 0.033, 0.075, 0.167, 0.217, 0.142,
+ 0.092, 0.067, 0.058, 0.050, 0.042,
+ 0.050, 0.050, 0.058, 0.083, 0.150, 0.167, 0.117,
+ 0.083, 0.058, 0.058, 0.067, 0.058,
+ 0.058, 0.050, 0.100, 0.142, 0.158, 0.092, 0.067,
+ 0.050, 0.042, 0.058, 0.083, 0.100,
+ 0.112, 0.098, 0.101, 0.125, 0.122, 0.072, 0.036,
+ 0.024, 0.039, 0.067, 0.102, 0.102,
+ 0.139, 0.092, 0.082, 0.099, 0.087, 0.039, 0.015,
+ 0.012, 0.036, 0.108, 0.159, 0.131,
+ 0.136, 0.105, 0.091, 0.057, 0.042, 0.034, 0.036,
+ 0.052, 0.068, 0.111, 0.132, 0.136,
+ 0.156, 0.154, 0.117, 0.119, 0.086, 0.044, 0.025,
+ 0.015, 0.025, 0.044, 0.089, 0.127,
+ 0.204, 0.163, 0.118, 0.102, 0.060, 0.030, 0.018,
+ 0.012, 0.023, 0.041, 0.087, 0.143,
+ 0.157, 0.130, 0.119, 0.094, 0.062, 0.042, 0.028,
+ 0.021, 0.035, 0.062, 0.099, 0.150,
+ 0.128, 0.142, 0.122, 0.128, 0.105, 0.065, 0.035,
+ 0.024, 0.031, 0.044, 0.074, 0.101,
+ 0.133, 0.126, 0.111, 0.110, 0.081, 0.056, 0.038,
+ 0.027, 0.042, 0.063, 0.098, 0.117,
+ 0.099, 0.100, 0.101, 0.099, 0.088, 0.078, 0.072,
+ 0.064, 0.064, 0.069, 0.076, 0.089),
+ ncol=12, byrow=TRUE)
+
+colnames(xref) = seq(1, 12, 1)
+row.names(xref) = c('GROUP11', 'GROUP10', 'GROUP9', 'GROUP8',
+ 'GROUP7', 'GROUUK', 'GROUP6', 'GROUP5',
+ 'GROUP4', 'GROUP3', 'GROUP2', 'GROUP1')
+
+# Computes the hydrograph of a station
+get_hydrograph = function (df_data, period=NULL, df_meta=NULL) {
+
+ # If there is a specified period
+ if (!is.null(period)) {
+ # Extracts only the data of this period
+ df_data = df_data[df_data$Date >= as.Date(period[1])
+ & df_data$Date <= as.Date(period[2]),]
+ }
+
+ # If there is the metadata
+ if (!is.null(df_meta)) {
+ # New column in metadata for hydrological regime
+ df_meta$regime_hydro = NA
+
+ # Get all different stations code
+ Code = levels(factor(df_meta$code))
+ # Number of stations
+ nCode = length(Code)
+
+ # Otherwise it is just a list of flow from one station
+ } else {
+ # Only one code is present
+ nCode = 1
+ }
+
+ # Blank tibble to store data
+ df_QM = tibble()
+ # For all accessible code
+ for (k in 1:nCode) {
+ # If there is the metadata
+ if (!is.null(df_meta)) {
+ # Gets the code
+ code = Code[k]
+ # Get the associated data
+ df_data_code = df_data[df_data$code == code,]
+ } else {
+ # The data are the date for the current code
+ df_data_code = df_data
+ }
+
+ # Gets a list of the month of the data as numeric
+ monthData = as.numeric(format(df_data_code$Date, "%m"))
+ # Blank list to stock month mean
+ QM_code = c()
+ # For all months
+ for (i in 1:12) {
+ # Gets all the flow data associated to the current month
+ data = df_data_code$Value[monthData == i]
+ # Averages the data
+ QM_code[i] = mean(data, na.rm=TRUE)
+ }
+
+ regime = 0
+ classRegime = ""
+ distance = rep(0, length(xref[,1]))
+ distancemin = 0
+ for (j in 1:length(xref[,1])) {
+ distance[j] = sum((QM_code / mean(QM_code) - xref[j,])^2)
+ }
+ regime = 1 + length(xref[,1]) - which.min(distance)
+ distancemin = distance[which.min(distance)]
+
+ if (regime < 7) {
+ classRegime = "Pluvial"
+ } else if (regime >= 7 & regime < 10) {
+ classRegime = "Transition"
+ } else if (regime >= 10) {
+ classRegime = "Nival glaciaire"
+ }
+
+ # If there is the metadata
+ if (!is.null(df_meta)) {
+ # Creates a temporary tibble to store hydrograph results
+ df_QMtmp = tibble(QM=QM_code, code=code)
+ # Stores it
+ df_QM = bind_rows(df_QM, df_QMtmp)
+ # Stores result of the hydrological regime
+ df_meta$regime_hydro[df_meta$code == code] = classRegime
+
+ # Otherwise
+ } else {
+ # No tibble needed
+ df_QM = QM_code
+ df_meta = classRegime
+ }
+ }
+ # Returns the hydrograph and meta data
+ return (list(QM=df_QM, meta=df_meta))
+}
+
+### 2.2. Break date
# Compute the break date of the flow data by station
get_break = function (df_data, df_meta, p_thresold=0.05) {
@@ -394,7 +589,7 @@ get_break = function (df_data, df_meta, p_thresold=0.05) {
return (df_break)
}
-### 2.2. Time gap
+### 2.3. Time gap
# Compute the time gap by station
get_lacune = function (df_data, df_meta) {
diff --git a/processing/format.R b/processing/format.R
index ec84328e6788769167500aa3b3273098fc3984fc..4fb3ab1c5ce1fefe143f0ef626f90f17dba8d906 100644
--- a/processing/format.R
+++ b/processing/format.R
@@ -182,11 +182,10 @@ prepare_date = function(df_XEx, df_Xlist) {
meanXEx_code = mean(XEx_code, na.rm=TRUE)
dXEx_code = abs(XEx_code - meanXEx_code)
stdXEx_code = sd(XEx_code, na.rm=TRUE)
- print(stdXEx_code*2)
- print(dXEx_code)
- XEx_code[dXEx_code >= 365/3] = XEx_code[dXEx_code >= 365/3] + 365
- # df_XEx$values[OkXEx_code] = XEx_code
-
+
+ XEx_code[dXEx_code >= stdXEx_code*3] =
+ XEx_code[dXEx_code >= stdXEx_code*3] + 365
+ df_XEx$values[OkXEx_code] = XEx_code
}
return (df_XEx)
diff --git a/script.R b/script.R
index 037d0c2d9886bf4fc372eb4a888455720c1d376c..4ab6dd138bc5df78ffcd4368ae2be1e35b3ba78c 100644
--- a/script.R
+++ b/script.R
@@ -96,17 +96,17 @@ INlistname =
## TREND ANALYSIS
# Time period to analyse
-periodAll = c("1800-01-01", "2019-12-31")
-periodSub = c("1968-01-01", "2019-12-31")
+periodAll = c("1800-01-01", "2020-12-31")
+periodSub = c("1968-01-01", "2020-12-31")
trend_period = list(periodAll, periodSub)
# Time period to mean
-period1 = c("1968-01-01", "1994-12-31")
-period2 = c("1995-01-01", "2019-12-31")
+period1 = c("1968-01-01", "1988-12-31")
+period2 = c("2000-01-01", "2020-12-31")
mean_period = list(period1, period2)
# p value thresold
-p_thresold = 0.1 #c(0.01, 0.05, 0.1)
+p_thresold = 0.1
## MAP
@@ -226,8 +226,11 @@ df_meta = df_join$meta
## 3. ANALYSE
-### 3.1. Compute gap parameters for stations
+### 3.1. Compute other parameters for stations
+# Time gap
df_meta = get_lacune(df_data, df_meta)
+# Hydrograph
+df_meta = get_hydrograph(df_data, df_meta, period=mean_period[[1]])$meta
### 3.2. Trend analysis
# QA trend
@@ -243,10 +246,15 @@ 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)
+# Start date for low water trend
+# 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,
+ period=trend_period,
+ p_thresold=p_thresold)
### 3.3. Break analysis
# df_break = get_break(res_QAtrend$data, df_meta)
@@ -282,26 +290,40 @@ df_shapefile = ini_shapefile(computer_data_path,
# filename_opt='time')
### 4.2. Analysis layout
-datasheet_layout(isplot=c(
- 'datasheet'
- # 'matrix',
+datasheet_layout(toplot=c(
+ # 'datasheet'
+ 'matrix'
# 'map'
),
df_meta=df_meta,
+
df_data=list(res_QAtrend$data,
res_QMNAtrend$data,
res_VCN10trend$data,
- res_dateVCN10trend$data),
+ # res_tINItrend$data,
+ res_tMIDtrend$data),
+
df_trend=list(res_QAtrend$trend,
res_QMNAtrend$trend,
res_VCN10trend$trend,
- res_dateVCN10trend$trend),
- var=list('QA', 'QMNA', 'VCN10', 'date du VCN10'),
- type=list('flow', 'flow', 'flow', 'date'),
+ # res_tINItrend$trend,
+ res_tMIDtrend$trend),
+
+ var=list('QA',
+ 'QMNA',
+ 'VCN10',
+ # 'tINI',
+ 'tMID'),
+
+ type=list('flow',
+ 'flow',
+ 'flow',
+ # 'date',
+ 'date'),
layout_matrix=matrix(c(1, 2, 3, 4), ncol=1),
- missRect=list(TRUE, TRUE, TRUE, TRUE),
+ missRect=TRUE,
trend_period=trend_period,
mean_period=mean_period,
info_header=TRUE,