Skip to content
GitLab
Commit 957260a9 authored by Heraut Louis's avatar Heraut Louis
Browse files

critique et debug regime

parent 5f04744b
master
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 176 additions and 134 deletions
+176 -134
plotting/datasheet.R
+ 38
− 25
View file @ 957260a9
...@@ -35,7 +35,7 @@ source('plotting/shortcut.R', encoding='UTF-8') ...@@ -35,7 +35,7 @@ source('plotting/shortcut.R', encoding='UTF-8')
# Manages datasheets creations for all stations. Makes the call to # Manages datasheets creations for all stations. Makes the call to
# the different headers, trend analysis graphs and realises arranging # the different headers, trend analysis graphs and realises arranging
# every plots. # every plots.
datasheet_panel = function (list_df2plot, df_meta, trend_period, mean_period, colorForce, info_header, time_header, foot_note, layout_matrix, info_height, time_ratio, var_ratio, foot_height, resources_path, logo_dir, AEAGlogo_file, INRAElogo_file, FRlogo_file, outdirTmp, df_page=NULL) { datasheet_panel = function (list_df2plot, df_meta, trend_period, mean_period, axis_xlim, colorForce, info_header, time_header, foot_note, layout_matrix, info_height, time_ratio, var_ratio, foot_height, resources_path, logo_dir, AEAGlogo_file, INRAElogo_file, FRlogo_file, outdirTmp, df_page=NULL) {
# The percentage of augmentation and diminution of the min # The percentage of augmentation and diminution of the min
# and max limits for y axis # and max limits for y axis
...@@ -174,13 +174,19 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, mean_period, co ...@@ -174,13 +174,19 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, mean_period, co
if (!is.null(time_header)) { if (!is.null(time_header)) {
# Extracts the data serie corresponding to the code # Extracts the data serie corresponding to the code
time_header_code = time_header[time_header$code == code,] time_header_code = time_header[time_header$code == code,]
# Gets the limits of the time serie
axis_xlim = c(min(time_header_code$Date), if (is.null(axis_xlim)) {
max(time_header_code$Date)) # Gets the limits of the time serie
axis_xlim_code = c(min(time_header_code$Date),
max(time_header_code$Date))
} else {
axis_xlim_code = axis_xlim
}
# Gets the time serie plot # Gets the time serie plot
Htime = time_panel(time_header_code, df_trend_code=NULL, Htime = time_panel(time_header_code, df_trend_code=NULL,
trend_period=trend_period, trend_period=trend_period,
axis_xlim=axis_xlim, missRect=TRUE, axis_xlim=axis_xlim_code, missRect=TRUE,
unit2day=365.25, var='Q', type='sévérité', unit2day=365.25, var='Q', type='sévérité',
grid=TRUE, ymin_lim=0, grid=TRUE, ymin_lim=0,
NspaceMax=NspaceMax[k], NspaceMax=NspaceMax[k],
...@@ -272,8 +278,6 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, mean_period, co ...@@ -272,8 +278,6 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, mean_period, co
color = append(color, colortmp) color = append(color, colortmp)
grid = FALSE grid = FALSE
} }
} else {
axis_xlim = NULL
} }
if (var != 'sqrt(Q)' & var != 'Q') { if (var != 'sqrt(Q)' & var != 'Q') {
...@@ -294,7 +298,8 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, mean_period, co ...@@ -294,7 +298,8 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, mean_period, co
p = time_panel(df_data_code, df_trend_code, var=var, p = time_panel(df_data_code, df_trend_code, var=var,
type=type, alpha=alpha, colorForce=colorForce, type=type, alpha=alpha, colorForce=colorForce,
missRect=missRect, trend_period=trend_period, missRect=missRect, trend_period=trend_period,
mean_period=mean_period, axis_xlim=axis_xlim, mean_period=mean_period,
axis_xlim=axis_xlim_code,
unit2day=unit2day, grid=grid, unit2day=unit2day, grid=grid,
ymin_lim=ymin_lim, color=color, ymin_lim=ymin_lim, color=color,
NspaceMax=NspaceMax[k], first=first, NspaceMax=NspaceMax[k], first=first,
...@@ -467,12 +472,12 @@ time_panel = function (df_data_code, df_trend_code, var, type, alpha=0.1, colorF ...@@ -467,12 +472,12 @@ time_panel = function (df_data_code, df_trend_code, var, type, alpha=0.1, colorF
i = i + 1 i = i + 1
} }
maxQ_lim = maxQtmp maxQ_lim = maxQtmp
# If x axis limits are specified # If x axis limits are specified
if (!is.null(axis_xlim)) { if (!is.null(axis_xlim)) {
minor_minDatetmp_lim = as.Date(axis_xlim[1]) axis_xlim = as.Date(axis_xlim)
minor_maxDatetmp_lim = as.Date(axis_xlim[2]) minor_minDatetmp_lim = axis_xlim[1]
minor_maxDatetmp_lim = axis_xlim[2]
# Otherwise # Otherwise
} else { } else {
minor_minDatetmp_lim = as.Date(df_data_code$Date[1]) minor_minDatetmp_lim = as.Date(df_data_code$Date[1])
...@@ -527,11 +532,11 @@ time_panel = function (df_data_code, df_trend_code, var, type, alpha=0.1, colorF ...@@ -527,11 +532,11 @@ time_panel = function (df_data_code, df_trend_code, var, type, alpha=0.1, colorF
listDate[which.min(abs(listDate - minor_minDatetmp_lim))] listDate[which.min(abs(listDate - minor_minDatetmp_lim))]
minor_maxDate_lim = minor_maxDate_lim =
listDate[which.min(abs(listDate - minor_maxDatetmp_lim))] listDate[which.min(abs(listDate - minor_maxDatetmp_lim))]
minor_minDate_lim = as.Date(paste(minor_minDate_lim, minor_minDate_lim = as.Date(paste(round(minor_minDate_lim),
'-01-01', sep='')) '-01-01', sep=''))
minor_maxDate_lim = as.Date(paste(minor_maxDate_lim, minor_maxDate_lim = as.Date(paste(round(minor_maxDate_lim),
'-01-01', sep='')) '-01-01', sep=''))
# Open new plot # Open new plot
p = ggplot() + theme_ash p = ggplot() + theme_ash
...@@ -1158,21 +1163,29 @@ time_panel = function (df_data_code, df_trend_code, var, type, alpha=0.1, colorF ...@@ -1158,21 +1163,29 @@ time_panel = function (df_data_code, df_trend_code, var, type, alpha=0.1, colorF
limits = axis_xlim limits = axis_xlim
} }
if (breakDate < 1) {
breaks = waiver()
minor_breaks = waiver()
date_labels = waiver()
} else {
breaks = seq(minDate_lim, maxDate_lim,
by=paste(breakDate, 'years'))
minor_breaks = seq(minor_minDate_lim,
minor_maxDate_lim,
by=paste(minor_breakDate,
'years'))
date_labels = "%Y"
}
# Parameters of the x axis contain the limit of the date data # Parameters of the x axis contain the limit of the date data
p = p + p = p +
scale_x_date(breaks=seq(minDate_lim, maxDate_lim, scale_x_date(breaks=breaks,
by=paste(breakDate, 'years')), minor_breaks=minor_breaks,
minor_breaks=seq(minor_minDate_lim,
minor_maxDate_lim,
by=paste(minor_breakDate,
'years')),
guide='axis_minor', guide='axis_minor',
date_labels="%Y", date_labels=date_labels,
limits=limits, limits=limits,
position=position, position=position,
expand=c(0, 0)) expand=c(0, 0))
# If it is a date variable # If it is a date variable
if (type == 'saisonnalité') { if (type == 'saisonnalité') {
......
plotting/layout.R
+ 19
− 18
View file @ 957260a9
...@@ -125,24 +125,24 @@ contour = void + ...@@ -125,24 +125,24 @@ contour = void +
## 3. LAYOUT _________________________________________________________ ## 3. LAYOUT _________________________________________________________
# Generates a PDF that gather datasheets, map and summarize matrix about the trend analyses realised on selected stations # 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, layout_panel = function (df_data, df_meta, layout_matrix,
toplot=c('datasheet', 'matrix', 'map'), toplot=c('datasheet', 'matrix', 'map'),
figdir='', filedir_opt='', filename_opt='', figdir='', filedir_opt='', filename_opt='',
variable='', df_trend=NULL, variable='', df_trend=NULL,
alpha=0.1, unit2day=365.25, var='', alpha=0.1, unit2day=365.25, var='',
type='', glose=NULL, trend_period=NULL, type='', glose=NULL, trend_period=NULL,
mean_period=NULL, colorForce=FALSE, mean_period=NULL, colorForce=FALSE,
axis_xlim=NULL, axis_xlim=NULL,
missRect=TRUE, time_header=NULL, missRect=TRUE, time_header=NULL,
info_header=NULL, foot_note=TRUE, info_header=NULL, foot_note=TRUE,
info_height=2.8, time_ratio=2, info_height=2.8, time_ratio=2,
var_ratio=3, foot_height=1.25, var_ratio=3, foot_height=1.25,
df_shapefile=NULL, df_shapefile=NULL,
resources_path=NULL, resources_path=NULL,
logo_dir=NULL, logo_dir=NULL,
AEAGlogo_file=NULL, AEAGlogo_file=NULL,
INRAElogo_file=NULL, INRAElogo_file=NULL,
FRlogo_file=NULL) { FRlogo_file=NULL) {
# Name of the document # Name of the document
outfile = "Panels" outfile = "Panels"
...@@ -291,6 +291,7 @@ datasheet_layout = function (df_data, df_meta, layout_matrix, ...@@ -291,6 +291,7 @@ datasheet_layout = function (df_data, df_meta, layout_matrix,
df_meta, df_meta,
trend_period=trend_period, trend_period=trend_period,
mean_period=mean_period, mean_period=mean_period,
axis_xlim=axis_xlim,
colorForce=colorForce, colorForce=colorForce,
info_header=info_header, info_header=info_header,
time_header=time_header, time_header=time_header,
......
processing/analyse.R
+ 1
− 1
View file @ 957260a9
...@@ -700,7 +700,7 @@ get_hydrograph = function (df_data, period=NULL, df_meta=NULL) { ...@@ -700,7 +700,7 @@ get_hydrograph = function (df_data, period=NULL, df_meta=NULL) {
for (j in 1:length(xref[,1])) { for (j in 1:length(xref[,1])) {
distance[j] = sum((QM_code / mean(QM_code) - xref[j,])^2) distance[j] = sum((QM_code / mean(QM_code) - xref[j,])^2)
} }
regime = 1 + length(xref[,1]) - which.min(distance) regime = which.min(distance)
distancemin = distance[which.min(distance)] distancemin = distance[which.min(distance)]
if (regime < 7) { if (regime < 7) {
......
processing/format.R
+ 20
− 8
View file @ 957260a9
...@@ -95,6 +95,15 @@ missing_data = function (df_data, df_meta, yearLac_day=3, yearStart='01-01', Cod ...@@ -95,6 +95,15 @@ missing_data = function (df_data, df_meta, yearLac_day=3, yearStart='01-01', Cod
for (code in Code) { for (code in Code) {
# Extracts the data corresponding to the code # Extracts the data corresponding to the code
df_data_code = df_data[df_data$code == code,] df_data_code = df_data[df_data$code == code,]
Date_NA = df_data_code$Date[is.na(df_data_code$Value)]
# print(Date_NA)
# Start_NA
DateMD = substr(df_data_code$Date, 6, 10) DateMD = substr(df_data_code$Date, 6, 10)
idyearStart = which(DateMD == yearStart) idyearStart = which(DateMD == yearStart)
...@@ -316,19 +325,22 @@ prepare_date = function(df_XEx, df_Xlist, per.start="01-01") { ...@@ -316,19 +325,22 @@ prepare_date = function(df_XEx, df_Xlist, per.start="01-01") {
origin = as.Date(paste(format(df_dateStart$Date, "%Y"), origin = as.Date(paste(format(df_dateStart$Date, "%Y"),
'-', per.start, sep='')) '-', per.start, sep=''))
originHydro = as.Date(paste(format(df_dateStart$Date, "%Y"), # originHydro = as.Date(paste(format(df_dateStart$Date, "%Y"),
'-01-01', sep='')) # '-01-01', sep=''))
for (i in 1:nrow(df_dateStart)) { for (i in 1:nrow(df_dateStart)) {
dateJultmp = julian(date[i], origin=origin[i]) dateJultmp = julian(date[i], origin=origin[i])
df_dateStart$Date_julian[i] = dateJultmp df_dateStart$Date_julian[i] = dateJultmp
dateJulHydrotmp = julian(date[i], origin=originHydro[i])
df_dateStart$DateHydro_julian[i] = dateJulHydrotmp # print(date[i])
# dateJulHydrotmp = julian(date[i], origin=originHydro[i])
# df_dateStart$DateHydro_julian[i] = dateJulHydrotmp
} }
df_dateStart$Year = format(df_dateStart$Date, "%Y") df_dateStart$Year = format(df_dateStart$Date, "%Y")
for (group in df_dateStart$group) { for (group in df_dateStart$group) {
Ok_dateStart = df_dateStart$group == group Ok_dateStart = df_dateStart$group == group
Shift = df_dateStart$Date_julian[Ok_dateStart] Shift = df_dateStart$Date_julian[Ok_dateStart]
year = df_dateStart$Year[Ok_dateStart] year = df_dateStart$Year[Ok_dateStart]
...@@ -336,11 +348,11 @@ prepare_date = function(df_XEx, df_Xlist, per.start="01-01") { ...@@ -336,11 +348,11 @@ prepare_date = function(df_XEx, df_Xlist, per.start="01-01") {
df_XEx$values[OkXEx_code_year] = df_XEx$values[OkXEx_code_year] =
df_XEx$values[OkXEx_code_year] + Shift df_XEx$values[OkXEx_code_year] + Shift
OkXEx_code = df_XEx$group1 == group # OkXEx_code = df_XEx$group1 == group
ShiftHydro = df_dateStart$DateHydro_julian[Ok_dateStart]
df_XEx$values[OkXEx_code] = df_XEx$values[OkXEx_code] + ShiftHydro
# ShiftHydro = df_dateStart$DateHydro_julian[Ok_dateStart]
# df_XEx$values[OkXEx_code] = df_XEx$values[OkXEx_code] + ShiftHydro
## Add 365 when the point is too remote ## Add 365 when the point is too remote
# XEx_code = df_XEx$values[OkXEx_code] # XEx_code = df_XEx$values[OkXEx_code]
# meanXEx_code = mean(XEx_code, na.rm=TRUE) # meanXEx_code = mean(XEx_code, na.rm=TRUE)
......
script.R
+ 98
− 82
View file @ 957260a9
...@@ -55,22 +55,22 @@ filedir = ...@@ -55,22 +55,22 @@ filedir =
# Name of the file that will be analysed from the BH directory # Name of the file that will be analysed from the BH directory
# (if 'all', all the file of the directory will be chosen) # (if 'all', all the file of the directory will be chosen)
filename = filename =
"" # ""
# "all" # "all"
# c( c(
# "S2235610_HYDRO_QJM.txt", # "S2235610_HYDRO_QJM.txt",
# "P1712910_HYDRO_QJM.txt", # "P0885010_HYDRO_QJM.txt"
# "P0885010_HYDRO_QJM.txt", # "O3006710_HYDRO_QJM.txt",
# "O5055010_HYDRO_QJM.txt", # "O4704030_HYDRO_QJM.txt"
# "O0384010_HYDRO_QJM.txt", # "O0384010_HYDRO_QJM.txt",
# "S4214010_HYDRO_QJM.txt", "O0362510_HYDRO_QJM.txt"
# "Q7002910_HYDRO_QJM.txt", # "Q7002910_HYDRO_QJM.txt"
# "Q0214010_HYDRO_QJM.txt", # "Q0214010_HYDRO_QJM.txt",
# "O3035210_HYDRO_QJM.txt", # "O3035210_HYDRO_QJM.txt",
# "O0554010_HYDRO_QJM.txt", # "O0554010_HYDRO_QJM.txt",
# "Q6332510_HYDRO_QJM.txt" # "Q6332510_HYDRO_QJM.txt"
# "O0362510_HYDRO_QJM.txt" # "O0362510_HYDRO_QJM.txt"
# ) )
## AGENCE EAU ADOUR GARONNE SELECTION ## AGENCE EAU ADOUR GARONNE SELECTION
...@@ -80,8 +80,8 @@ AEAGlistdir = ...@@ -80,8 +80,8 @@ AEAGlistdir =
"" ""
AEAGlistname = AEAGlistname =
# "" ""
"Liste-station_RRSE.docx" # "Liste-station_RRSE.docx"
## NIVALE SELECTION ## NIVALE SELECTION
...@@ -96,7 +96,17 @@ INRAElistname = ...@@ -96,7 +96,17 @@ INRAElistname =
# "INRAE_selection.txt" # "INRAE_selection.txt"
## TREND ANALYSIS which_layout =
# ('serie')
('analyse')
# ('serie', 'analyse')
# Selection
axis_xlim =
NULL
# c("2002-01-01", "2003-01-01")
## ANALYSIS
# Time period to analyse # Time period to analyse
periodAll = c("1800-01-01", "2020-12-31") periodAll = c("1800-01-01", "2020-12-31")
periodSub = c("1968-01-01", "2020-12-31") periodSub = c("1968-01-01", "2020-12-31")
...@@ -119,7 +129,7 @@ sampleSpan = c('05-01', '11-30') ...@@ -119,7 +129,7 @@ sampleSpan = c('05-01', '11-30')
## MAP ## MAP
# Is the hydrological network needs to be plot # Is the hydrological network needs to be plot
is_river = TRUE is_river = FALSE
############### END OF REGION TO MODIFY (without risk) ############### ############### END OF REGION TO MODIFY (without risk) ###############
...@@ -274,56 +284,58 @@ df_meta = get_hydrograph(df_data, df_meta, period=mean_period[[1]])$meta ...@@ -274,56 +284,58 @@ df_meta = get_hydrograph(df_data, df_meta, period=mean_period[[1]])$meta
df_sqrt = compute_sqrt(df_data) df_sqrt = compute_sqrt(df_data)
### 3.2. Trend analysis ______________________________________________ ### 3.2. Trend analysis ______________________________________________
# # QA trend if ('analyse' %in% which_layout) {
# res = get_QAtrend(df_data, df_meta, # QA trend
# period=trend_period, res = get_QAtrend(df_data, df_meta,
# alpha=alpha, period=trend_period,
# yearLac_day=yearLac_day) alpha=alpha,
# df_QAdata = res$data yearLac_day=yearLac_day)
# df_QAmod = res$mod df_QAdata = res$data
# res_QAtrend = res$analyse df_QAmod = res$mod
res_QAtrend = res$analyse
# # QMNA tend
# res = get_QMNAtrend(df_data, df_meta, # QMNA tend
# period=trend_period, res = get_QMNAtrend(df_data, df_meta,
# alpha=alpha, period=trend_period,
# sampleSpan=sampleSpan, alpha=alpha,
# yearLac_day=yearLac_day) sampleSpan=sampleSpan,
# df_QMNAdata = res$data yearLac_day=yearLac_day)
# df_QMNAmod = res$mod df_QMNAdata = res$data
# res_QMNAtrend = res$analyse df_QMNAmod = res$mod
res_QMNAtrend = res$analyse
# # VCN10 trend
# res = get_VCN10trend(df_data, df_meta, # VCN10 trend
# period=trend_period, res = get_VCN10trend(df_data, df_meta,
# alpha=alpha, period=trend_period,
# sampleSpan=sampleSpan, alpha=alpha,
# yearLac_day=yearLac_day) sampleSpan=sampleSpan,
# df_VCN10data = res$data yearLac_day=yearLac_day)
# df_VCN10mod = res$mod df_VCN10data = res$data
# res_VCN10trend = res$analyse df_VCN10mod = res$mod
res_VCN10trend = res$analyse
# # Start date for low water trend
# res = get_tDEBtrend(df_data, df_meta, # Start date for low water trend
# period=trend_period, res = get_tDEBtrend(df_data, df_meta,
# alpha=alpha, period=trend_period,
# sampleSpan=sampleSpan, alpha=alpha,
# thresold_type='VCN10', sampleSpan=sampleSpan,
# select_longest=TRUE, thresold_type='VCN10',
# yearLac_day=yearLac_day) select_longest=TRUE,
# df_tDEBdata = res$data yearLac_day=yearLac_day)
# df_tDEBmod = res$mod df_tDEBdata = res$data
# res_tDEBtrend = res$analyse df_tDEBmod = res$mod
res_tDEBtrend = res$analyse
# # Center date for low water trend
# res = get_tCENtrend(df_data, df_meta, # Center date for low water trend
# period=trend_period, res = get_tCENtrend(df_data, df_meta,
# alpha=alpha, period=trend_period,
# sampleSpan=sampleSpan, alpha=alpha,
# yearLac_day=yearLac_day) sampleSpan=sampleSpan,
# df_tCENdata = res$data yearLac_day=yearLac_day)
# df_tCENmod = res$mod df_tCENdata = res$data
# res_tCENtrend = res$analyse df_tCENmod = res$mod
res_tCENtrend = res$analyse
}
### 3.3. Break analysis ______________________________________________ ### 3.3. Break analysis ______________________________________________
# df_break = get_break(res_QAtrend$data, df_meta) # df_break = get_break(res_QAtrend$data, df_meta)
...@@ -362,28 +374,31 @@ df_shapefile = ini_shapefile(resources_path, ...@@ -362,28 +374,31 @@ df_shapefile = ini_shapefile(resources_path,
### 5.1. Simple time panel to criticize station data _________________ ### 5.1. Simple time panel to criticize station data _________________
# Plot time panel of debit by stations # Plot time panel of debit by stations
# datasheet_layout(toplot=c('datasheet'), if ('serie' %in% which_layout) {
# df_meta=df_meta, layout_panel(toplot=c('datasheet'),
# df_data=list(df_data, df_meta=df_meta,
# df_sqrt), df_data=list(df_data,
# var=list('Q', 'sqrt(Q)'), df_sqrt),
# type=list('data', 'data'), var=list('Q', 'sqrt(Q)'),
# layout_matrix=matrix(c(1, 2), ncol=1), type=list('data', 'data'),
# info_header=df_data, axis_xlim=axis_xlim,
# df_shapefile=df_shapefile, layout_matrix=matrix(c(1, 2), ncol=1),
# figdir=figdir, info_header=df_data,
# resources_path=resources_path, df_shapefile=df_shapefile,
# logo_dir=logo_dir, figdir=figdir,
# AEAGlogo_file=AEAGlogo_file, resources_path=resources_path,
# INRAElogo_file=INRAElogo_file, logo_dir=logo_dir,
# FRlogo_file=FRlogo_file) AEAGlogo_file=AEAGlogo_file,
INRAElogo_file=INRAElogo_file,
FRlogo_file=FRlogo_file)
}
### 5.2. Analysis layout _____________________________________________ ### 5.2. Analysis layout _____________________________________________
datasheet_layout(toplot=c( if ('analyse' %in% which_layout) {
'datasheet', layout_panel(toplot=c(
'matrix', 'datasheet'
'map' # 'matrix',
# 'map'
), ),
df_meta=df_meta, df_meta=df_meta,
...@@ -428,3 +443,4 @@ datasheet_layout(toplot=c( ...@@ -428,3 +443,4 @@ datasheet_layout(toplot=c(
AEAGlogo_file=AEAGlogo_file, AEAGlogo_file=AEAGlogo_file,
INRAElogo_file=INRAElogo_file, INRAElogo_file=INRAElogo_file,
FRlogo_file=FRlogo_file) FRlogo_file=FRlogo_file)
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment