diff --git a/plotting/datasheet.R b/plotting/datasheet.R index eca2270d031ab85fd687de122896abefc5c714d1..53d670d9753c848890db0e5ebe3ea19f3ef42007 100644 --- a/plotting/datasheet.R +++ b/plotting/datasheet.R @@ -35,7 +35,7 @@ source('plotting/shortcut.R', encoding='UTF-8') # Manages datasheets creations for all stations. Makes the call to # the different headers, trend analysis graphs and realises arranging # 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 # and max limits for y axis @@ -174,13 +174,19 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, mean_period, co if (!is.null(time_header)) { # Extracts the data serie corresponding to the 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), - max(time_header_code$Date)) + + if (is.null(axis_xlim)) { + # 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 Htime = time_panel(time_header_code, df_trend_code=NULL, 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é', grid=TRUE, ymin_lim=0, NspaceMax=NspaceMax[k], @@ -272,8 +278,6 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, mean_period, co color = append(color, colortmp) grid = FALSE } - } else { - axis_xlim = NULL } if (var != 'sqrt(Q)' & var != 'Q') { @@ -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, type=type, alpha=alpha, colorForce=colorForce, 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, ymin_lim=ymin_lim, color=color, NspaceMax=NspaceMax[k], first=first, @@ -467,12 +472,12 @@ time_panel = function (df_data_code, df_trend_code, var, type, alpha=0.1, colorF i = i + 1 } maxQ_lim = maxQtmp - - + # If x axis limits are specified if (!is.null(axis_xlim)) { - minor_minDatetmp_lim = as.Date(axis_xlim[1]) - minor_maxDatetmp_lim = as.Date(axis_xlim[2]) + axis_xlim = as.Date(axis_xlim) + minor_minDatetmp_lim = axis_xlim[1] + minor_maxDatetmp_lim = axis_xlim[2] # Otherwise } else { 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 listDate[which.min(abs(listDate - minor_minDatetmp_lim))] minor_maxDate_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='')) - minor_maxDate_lim = as.Date(paste(minor_maxDate_lim, + minor_maxDate_lim = as.Date(paste(round(minor_maxDate_lim), '-01-01', sep='')) - + # Open new plot p = ggplot() + theme_ash @@ -1158,21 +1163,29 @@ time_panel = function (df_data_code, df_trend_code, var, type, alpha=0.1, colorF 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 p = p + - scale_x_date(breaks=seq(minDate_lim, maxDate_lim, - by=paste(breakDate, 'years')), - minor_breaks=seq(minor_minDate_lim, - minor_maxDate_lim, - by=paste(minor_breakDate, - 'years')), + scale_x_date(breaks=breaks, + minor_breaks=minor_breaks, guide='axis_minor', - date_labels="%Y", + date_labels=date_labels, limits=limits, position=position, - expand=c(0, 0)) - - + expand=c(0, 0)) # If it is a date variable if (type == 'saisonnalité') { diff --git a/plotting/layout.R b/plotting/layout.R index 4242719dab5c8a7d8d0e2b00a1606412632e01bb..b1da9cf972f30815f331e7a7f281e06c0b48919e 100644 --- a/plotting/layout.R +++ b/plotting/layout.R @@ -125,24 +125,24 @@ contour = void + ## 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, - toplot=c('datasheet', 'matrix', 'map'), - figdir='', filedir_opt='', filename_opt='', - variable='', df_trend=NULL, - alpha=0.1, unit2day=365.25, var='', - type='', glose=NULL, trend_period=NULL, - mean_period=NULL, colorForce=FALSE, - axis_xlim=NULL, - missRect=TRUE, time_header=NULL, - info_header=NULL, foot_note=TRUE, - info_height=2.8, time_ratio=2, - var_ratio=3, foot_height=1.25, - df_shapefile=NULL, - resources_path=NULL, - logo_dir=NULL, - AEAGlogo_file=NULL, - INRAElogo_file=NULL, - FRlogo_file=NULL) { +layout_panel = function (df_data, df_meta, layout_matrix, + toplot=c('datasheet', 'matrix', 'map'), + figdir='', filedir_opt='', filename_opt='', + variable='', df_trend=NULL, + alpha=0.1, unit2day=365.25, var='', + type='', glose=NULL, trend_period=NULL, + mean_period=NULL, colorForce=FALSE, + axis_xlim=NULL, + missRect=TRUE, time_header=NULL, + info_header=NULL, foot_note=TRUE, + info_height=2.8, time_ratio=2, + var_ratio=3, foot_height=1.25, + df_shapefile=NULL, + resources_path=NULL, + logo_dir=NULL, + AEAGlogo_file=NULL, + INRAElogo_file=NULL, + FRlogo_file=NULL) { # Name of the document outfile = "Panels" @@ -291,6 +291,7 @@ datasheet_layout = function (df_data, df_meta, layout_matrix, df_meta, trend_period=trend_period, mean_period=mean_period, + axis_xlim=axis_xlim, colorForce=colorForce, info_header=info_header, time_header=time_header, diff --git a/processing/analyse.R b/processing/analyse.R index 87b1d3d96a86654b31824f57408014dca28501df..42540eee35309f58ff4a47e32964aeeac1d98696 100644 --- a/processing/analyse.R +++ b/processing/analyse.R @@ -700,7 +700,7 @@ get_hydrograph = function (df_data, period=NULL, df_meta=NULL) { 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) + regime = which.min(distance) distancemin = distance[which.min(distance)] if (regime < 7) { diff --git a/processing/format.R b/processing/format.R index 11ca6fe6a1cfecd436fa954ce4d49350c8567194..88e390b6081094257c0a7d998f5c69f3cf8372f8 100644 --- a/processing/format.R +++ b/processing/format.R @@ -95,6 +95,15 @@ missing_data = function (df_data, df_meta, yearLac_day=3, yearStart='01-01', Cod for (code in Code) { # Extracts the data corresponding to the 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) idyearStart = which(DateMD == yearStart) @@ -316,19 +325,22 @@ prepare_date = function(df_XEx, df_Xlist, per.start="01-01") { origin = as.Date(paste(format(df_dateStart$Date, "%Y"), '-', per.start, sep='')) - originHydro = as.Date(paste(format(df_dateStart$Date, "%Y"), - '-01-01', sep='')) + # originHydro = as.Date(paste(format(df_dateStart$Date, "%Y"), + # '-01-01', sep='')) for (i in 1:nrow(df_dateStart)) { dateJultmp = julian(date[i], origin=origin[i]) 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") for (group in df_dateStart$group) { + Ok_dateStart = df_dateStart$group == group Shift = df_dateStart$Date_julian[Ok_dateStart] year = df_dateStart$Year[Ok_dateStart] @@ -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] + Shift - OkXEx_code = df_XEx$group1 == group - - ShiftHydro = df_dateStart$DateHydro_julian[Ok_dateStart] - df_XEx$values[OkXEx_code] = df_XEx$values[OkXEx_code] + ShiftHydro + # OkXEx_code = df_XEx$group1 == group + # 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 # XEx_code = df_XEx$values[OkXEx_code] # meanXEx_code = mean(XEx_code, na.rm=TRUE) diff --git a/script.R b/script.R index 8ce0cc05483a5dd404e9fec5ad5db9689aca861d..18d5b49edacd6083148266b23c62ad842a68c082 100644 --- a/script.R +++ b/script.R @@ -55,22 +55,22 @@ 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 = - "" + # "" # "all" - # c( + c( # "S2235610_HYDRO_QJM.txt", - # "P1712910_HYDRO_QJM.txt", - # "P0885010_HYDRO_QJM.txt", - # "O5055010_HYDRO_QJM.txt", + # "P0885010_HYDRO_QJM.txt" + # "O3006710_HYDRO_QJM.txt", + # "O4704030_HYDRO_QJM.txt" # "O0384010_HYDRO_QJM.txt", - # "S4214010_HYDRO_QJM.txt", - # "Q7002910_HYDRO_QJM.txt", + "O0362510_HYDRO_QJM.txt" + # "Q7002910_HYDRO_QJM.txt" # "Q0214010_HYDRO_QJM.txt", # "O3035210_HYDRO_QJM.txt", # "O0554010_HYDRO_QJM.txt", # "Q6332510_HYDRO_QJM.txt" # "O0362510_HYDRO_QJM.txt" - # ) + ) ## AGENCE EAU ADOUR GARONNE SELECTION @@ -80,8 +80,8 @@ AEAGlistdir = "" AEAGlistname = - # "" - "Liste-station_RRSE.docx" + "" + # "Liste-station_RRSE.docx" ## NIVALE SELECTION @@ -96,7 +96,17 @@ INRAElistname = # "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 periodAll = c("1800-01-01", "2020-12-31") periodSub = c("1968-01-01", "2020-12-31") @@ -119,7 +129,7 @@ sampleSpan = c('05-01', '11-30') ## MAP # Is the hydrological network needs to be plot -is_river = TRUE +is_river = FALSE ############### END OF REGION TO MODIFY (without risk) ############### @@ -274,56 +284,58 @@ df_meta = get_hydrograph(df_data, df_meta, period=mean_period[[1]])$meta df_sqrt = compute_sqrt(df_data) ### 3.2. Trend analysis ______________________________________________ -# # QA trend -# res = get_QAtrend(df_data, df_meta, -# period=trend_period, -# alpha=alpha, -# yearLac_day=yearLac_day) -# df_QAdata = res$data -# df_QAmod = res$mod -# res_QAtrend = res$analyse - -# # QMNA tend -# res = get_QMNAtrend(df_data, df_meta, -# period=trend_period, -# alpha=alpha, -# sampleSpan=sampleSpan, -# yearLac_day=yearLac_day) -# df_QMNAdata = res$data -# df_QMNAmod = res$mod -# res_QMNAtrend = res$analyse - -# # VCN10 trend -# res = get_VCN10trend(df_data, df_meta, -# period=trend_period, -# alpha=alpha, -# sampleSpan=sampleSpan, -# yearLac_day=yearLac_day) -# df_VCN10data = res$data -# df_VCN10mod = res$mod -# res_VCN10trend = res$analyse - -# # Start date for low water trend -# res = get_tDEBtrend(df_data, df_meta, -# period=trend_period, -# alpha=alpha, -# sampleSpan=sampleSpan, -# thresold_type='VCN10', -# select_longest=TRUE, -# yearLac_day=yearLac_day) -# df_tDEBdata = res$data -# df_tDEBmod = res$mod -# res_tDEBtrend = res$analyse - -# # Center date for low water trend -# res = get_tCENtrend(df_data, df_meta, -# period=trend_period, -# alpha=alpha, -# sampleSpan=sampleSpan, -# yearLac_day=yearLac_day) -# df_tCENdata = res$data -# df_tCENmod = res$mod -# res_tCENtrend = res$analyse +if ('analyse' %in% which_layout) { + # QA trend + res = get_QAtrend(df_data, df_meta, + period=trend_period, + alpha=alpha, + yearLac_day=yearLac_day) + df_QAdata = res$data + df_QAmod = res$mod + res_QAtrend = res$analyse + + # QMNA tend + res = get_QMNAtrend(df_data, df_meta, + period=trend_period, + alpha=alpha, + sampleSpan=sampleSpan, + yearLac_day=yearLac_day) + df_QMNAdata = res$data + df_QMNAmod = res$mod + res_QMNAtrend = res$analyse + + # VCN10 trend + res = get_VCN10trend(df_data, df_meta, + period=trend_period, + alpha=alpha, + sampleSpan=sampleSpan, + yearLac_day=yearLac_day) + df_VCN10data = res$data + df_VCN10mod = res$mod + res_VCN10trend = res$analyse + + # Start date for low water trend + res = get_tDEBtrend(df_data, df_meta, + period=trend_period, + alpha=alpha, + sampleSpan=sampleSpan, + thresold_type='VCN10', + select_longest=TRUE, + yearLac_day=yearLac_day) + df_tDEBdata = res$data + df_tDEBmod = res$mod + res_tDEBtrend = res$analyse + + # Center date for low water trend + res = get_tCENtrend(df_data, df_meta, + period=trend_period, + alpha=alpha, + sampleSpan=sampleSpan, + yearLac_day=yearLac_day) + df_tCENdata = res$data + df_tCENmod = res$mod + res_tCENtrend = res$analyse +} ### 3.3. Break analysis ______________________________________________ # df_break = get_break(res_QAtrend$data, df_meta) @@ -362,28 +374,31 @@ df_shapefile = ini_shapefile(resources_path, ### 5.1. Simple time panel to criticize station data _________________ # Plot time panel of debit by stations -# datasheet_layout(toplot=c('datasheet'), -# df_meta=df_meta, -# df_data=list(df_data, -# df_sqrt), -# var=list('Q', 'sqrt(Q)'), -# type=list('data', 'data'), -# layout_matrix=matrix(c(1, 2), ncol=1), -# info_header=df_data, -# df_shapefile=df_shapefile, -# figdir=figdir, -# resources_path=resources_path, -# logo_dir=logo_dir, -# AEAGlogo_file=AEAGlogo_file, -# INRAElogo_file=INRAElogo_file, -# FRlogo_file=FRlogo_file) - +if ('serie' %in% which_layout) { + layout_panel(toplot=c('datasheet'), + df_meta=df_meta, + df_data=list(df_data, + df_sqrt), + var=list('Q', 'sqrt(Q)'), + type=list('data', 'data'), + axis_xlim=axis_xlim, + layout_matrix=matrix(c(1, 2), ncol=1), + info_header=df_data, + df_shapefile=df_shapefile, + figdir=figdir, + resources_path=resources_path, + logo_dir=logo_dir, + AEAGlogo_file=AEAGlogo_file, + INRAElogo_file=INRAElogo_file, + FRlogo_file=FRlogo_file) +} ### 5.2. Analysis layout _____________________________________________ -datasheet_layout(toplot=c( - 'datasheet', - 'matrix', - 'map' +if ('analyse' %in% which_layout) { + layout_panel(toplot=c( + 'datasheet' + # 'matrix', + # 'map' ), df_meta=df_meta, @@ -428,3 +443,4 @@ datasheet_layout(toplot=c( AEAGlogo_file=AEAGlogo_file, INRAElogo_file=INRAElogo_file, FRlogo_file=FRlogo_file) +}