Skip to content
GitLab
Select Git revision
  • master default
Find file BlameHistoryPermalink
script.R 15.11 KiB
1 
# \\\
2 
# Copyright 2021-2022 Louis Héraut*1
3 
#
4 
# *1   INRAE, France
5 
#      louis.heraut@inrae.fr
6 
#
7 
# This file is part of ash R toolbox.
8 
#
9 
# ash R toolbox is free software: you can redistribute it and/or modify
10 
# it under the terms of the GNU General Public License as published by
11 
# the Free Software Foundation, either version 3 of the License, or (at
12 
# your option) any later version.
13 
#
14 
# ash R toolbox is distributed in the hope that it will be useful, but
15 
# WITHOUT ANY WARRANTY; without even the implied warranty of
16 
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17 
# General Public License for more details.
18 
#
19 
# You should have received a copy of the GNU General Public License
20 
# along with ash R toolbox.  If not, see <https://www.gnu.org/licenses/>.
21 
# ///
22 
#
23 
#
24 
# script.R
25 
#
26 
# Script file to manage the trend analysis of the Adour-Garonne basin.
27 
# Performs the necessary calls to processing and plotting functions in
28 
# order to realise the hydrologic trend analysis of stations according
29 
# to the input parameters. The nearest area belove is where you need to
30 
# write your prefer parameters for the analysis. See the 'README.txt'
31 
# file for more information.
32
33
34 
############## START OF REGION TO MODIFY (without risk) ##############
35 
# Path to the data
36 
computer_data_path =
37 
    "/home/louis/Documents/bouleau/INRAE/CDD_stationnarite/data"
38 
    # "C:\\Users\\louis.heraut\\Documents\\CDD_stationnarite\\data"
39
40 
# Work path (it needs to end with '/ASH' directory)
41 
computer_work_path =
42 
    "/home/louis/Documents/bouleau/INRAE/CDD_stationnarite/ASH"
43 
    # "C:\\Users\\louis.heraut\\Documents\\CDD_stationnarite\\ASH"
44
45
46 
## BANQUE HYDRO
47 
# Path to the directory where Banque Hydro (BH) data is stored
48 
# from the work path
49 
filedir =
50 
    # ""
51 
    "BanqueHydro_Export2021"
52
53
54 
## MANUAL SELECTION
55 
# Name of the file that will be analysed from the BH directory
56 
# (if 'all', all the file of the directory will be chosen)
57 
filename =
58 
    # ""
59 
    # "all"
60 
    c(
61 
        "S2235610_HYDRO_QJM.txt",
62 
        "P0885010_HYDRO_QJM.txt",
63 
        "P0364010_HYDRO_QJM.txt",
64 
        "O7635010_HYDRO_QJM.txt",
65 
        "O3141010_HYDRO_QJM.txt",
66 
        "Q6332510_HYDRO_QJM.txt",
67 
        "Q7002910_HYDRO_QJM.txt"
68 
        # "Q0214010_HYDRO_QJM.txt",
69 
        # "O3035210_HYDRO_QJM.txt",
70 
        # "O0554010_HYDRO_QJM.txt",
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
# "Q6332510_HYDRO_QJM.txt" # "O0362510_HYDRO_QJM.txt" ) ## AGENCE EAU ADOUR GARONNE SELECTION # Path to the 'docx' list file of station from the Agence de l'eau # Adour-Garonne that will be analysed AEAGlistdir = "" AEAGlistname = "" # "Liste-station_RRSE.docx" ## NIVALE SELECTION # Path to the 'txt' list file of station from INRAE that will be analysed # Generated with : # create_selection(computer_data_path, 'dirname', 'example.txt') INRAElistdir = "" INRAElistname = "" # "INRAE_selection.txt" which_layout = # ('serie') ('analyse') # ('serie', 'analyse') # Selection axis_xlim = NULL # c("1982-01-01", "1983-01-01") ## ANALYSIS # Time period to analyse 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", "1988-12-31") period2 = c("2000-01-01", "2020-12-31") mean_period = list(period1, period2) # alpha the risk alpha = 0.1 # Number of missing days per year before remove the year dayLac_lim = 3 # Maximum continuously missing years before removing everything # before it yearNA_lim = 10 # Local corrections of the data df_flag = tibble( code=c('O3141010', 'O7635010', 'O7635010', 'O7635010', 'O7635010' ), Date=c('1974-07-04', '1948-09-06', '1949-02-08',
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
'1950-07-20', '1953-07-22' ), newValue=c(9.5, 4, 3, 1, 3) # /!\ Unit ) # Sampling span of the data sampleSpan = c('05-01', '11-30') # Is the hydrological network needs to be plot show_river = FALSE # If results and data used in the analysis needs to be written saving = FALSE ############### END OF REGION TO MODIFY (without risk) ############### ## 1. FILE STRUCTURE _________________________________________________ # Set working directory setwd(computer_work_path) # Sourcing R file source('processing/extract.R', encoding='UTF-8') source('processing/format.R', encoding='UTF-8') source('processing/analyse.R', encoding='UTF-8') source('plotting/layout.R', encoding='UTF-8') source('processing/read_write.R', encoding='UTF-8') # Result directory resdir = file.path(computer_work_path, 'results') if (!(file.exists(resdir))) { dir.create(resdir) } print(paste('resdir :', resdir)) # Figure directory figdir = file.path(computer_work_path, 'figures') if (!(file.exists(figdir))) { dir.create(figdir) } print(paste('figdir :', figdir)) # Resources directory resources_path = file.path(computer_work_path, 'resources') if (!(file.exists(resources_path))) { dir.create(resources_path) } print(paste('resources_path :', resources_path)) # Logo filename logo_dir = 'logo' AEAGlogo_file = 'agence-de-leau-adour-garonne_logo.png' INRAElogo_file = 'Logo-INRAE_Transparent.png' FRlogo_file = 'Republique_Francaise_RVB.png' # Path to the shapefile for france contour from 'computer_data_path' fr_shpdir = 'map/france' fr_shpname = 'gadm36_FRA_0.shp' # Path to the shapefile for basin shape from 'computer_data_path' bs_shpdir = 'map/bassin' bs_shpname = 'BassinHydrographique.shp' # Path to the shapefile for sub-basin shape from 'computer_data_path' sbs_shpdir = 'map/sous_bassin'
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
sbs_shpname = 'SousBassinHydrographique.shp' # Path to the shapefile for river shape from 'computer_data_path' rv_shpdir = 'map/river' rv_shpname = 'CoursEau_FXX.shp' ## 2. SELECTION OF STATION ___________________________________________ # Initialization of null data frame if there is no data selected df_data_AEAG = NULL df_data_INRAE = NULL df_meta_AEAG = NULL df_meta_INRAE = NULL ### 2.1. Selection of the Agence de l'eau Adour-Garonne ______________ if (AEAGlistname != "") { # Get only the selected station from a list station file df_selec_AEAG = get_selection_AEAG(computer_data_path, AEAGlistdir, AEAGlistname, cnames=c('code', 'station', 'BV_km2', 'axe_principal_concerne', 'longueur_serie', 'commentaires', 'choix'), c_num=c('BV_km2', 'longueur_serie')) # Get filenames of the selection filename = df_selec_AEAG[df_selec_AEAG$ok,]$filename # Extract metadata about selected stations df_meta_AEAG = extract_meta(computer_data_path, filedir, filename) # Extract data about selected stations df_data_AEAG = extract_data(computer_data_path, filedir, filename) } ### 2.2. INRAE selection _____________________________________________ if (INRAElistname != ""){ # Get only the selected station from a list station file df_selec_INRAE = get_selection_INRAE(computer_data_path, INRAElistdir, INRAElistname) # Get filenames of the selection filename = df_selec_INRAE[df_selec_INRAE$ok,]$filename # Extract metadata about selected stations df_meta_INRAE = extract_meta(computer_data_path, filedir, filename) # Extract data about selected stations df_data_INRAE = extract_data(computer_data_path, filedir, filename) } ### 2.3. Manual selection ____________________________________________ if (AEAGlistname == "" & INRAElistname == "") { # Extract metadata about selected stations df_meta_AEAG = extract_meta(computer_data_path, filedir, filename) # Extract data about selected stations df_data_AEAG = extract_data(computer_data_path, filedir, filename) } ### 2.4. Data join ___________________________________________________ df_join = join_selection(df_data_AEAG, df_data_INRAE, df_meta_AEAG, df_meta_INRAE) df_data = df_join$data df_meta = df_join$meta
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
## 3. ANALYSE ________________________________________________________ var = list( 'QA', 'QMNA', 'VCN10', 'tDEB', 'tCEN' ) type = list( 'sévérité', 'sévérité', 'sévérité', 'saisonnalité', 'saisonnalité' ) glose = list( "Moyenne annuelle du débit journalier", "Minimum annuel de la moyenne mensuelle du débit journalier", "Minimum annuel de la moyenne sur 10 jours du débit journalier", "Début d'étiage (jour de l'année de la première moyenne sur 10 jours sous le maximum des VCN10)", "Centre d'étiage (jour de l'année du VCN10)" ) ### 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 # Square root df_sqrt = compute_sqrt(df_data) ### 3.2. Trend analysis ______________________________________________ if ('analyse' %in% which_layout) { # QA trend res = get_QAtrend(df_data, df_meta, period=trend_period, alpha=alpha, dayLac_lim=dayLac_lim, yearNA_lim=yearNA_lim, df_flag=df_flag) 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, dayLac_lim=dayLac_lim, yearNA_lim=yearNA_lim, df_flag=df_flag) 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, dayLac_lim=dayLac_lim, yearNA_lim=yearNA_lim, df_flag=df_flag) 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,
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
period=trend_period, alpha=alpha, sampleSpan=sampleSpan, thresold_type='VCN10', select_longest=TRUE, dayLac_lim=dayLac_lim, yearNA_lim=yearNA_lim, df_flag=df_flag) 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, dayLac_lim=dayLac_lim, yearNA_lim=yearNA_lim, df_flag=df_flag) 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) # df_break = get_break(res_QMNAtrend$data, df_meta) # df_break = get_break(res_VCN10trend$data, df_meta) # histogram(df_break$Date, df_meta, # figdir=figdir) # cumulative(df_break$Date, df_meta, dyear=8, # figdir=figdir) ## 4. SAVING _________________________________________________________ if (saving) { for (v in var) { df_datatmp = get(paste('df_', v, 'data', sep='')) df_modtmp = get(paste('df_', v, 'mod', sep='')) res_trendtmp = get(paste('res_', v, 'trend', sep='')) # Modified data saving write_data(df_datatmp, df_modtmp, resdir, optdir='modified_data', filedir=v) # Trend analysis saving write_analyse(res_trendtmp, resdir, optdir='trend_analyse', filedir=v) } } # res_tDEBtrend = read_listofdf(resdir, 'res_tDEBtrend') ## 5. PLOTTING _______________________________________________________ # Shapefile importation in order to it only once time df_shapefile = ini_shapefile(resources_path, fr_shpdir, fr_shpname, bs_shpdir, bs_shpname, sbs_shpdir, sbs_shpname, rv_shpdir, rv_shpname, show_river=show_river) ### 5.1. Simple time panel to criticize station data _________________ # Plot time panel of debit by stations if ('serie' %in% which_layout) { layout_panel(what_plot=c('datasheet'), df_meta=df_meta, df_data=list(df_data, df_sqrt), var=list('Q', 'sqrt(Q)'),
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485
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 _____________________________________________ if ('analyse' %in% which_layout) { layout_panel(what_plot=c( # 'datasheet' # 'matrix', 'map' ), df_meta=df_meta, df_data=list( res_QAtrend$data, res_QMNAtrend$data, res_VCN10trend$data, res_tDEBtrend$data, res_tCENtrend$data ), df_trend=list( res_QAtrend$trend, res_QMNAtrend$trend, res_VCN10trend$trend, res_tDEBtrend$trend, res_tCENtrend$trend ), var=var, type=type, glose=glose, layout_matrix=matrix(c(1, 2, 3, 4, 5), ncol=1), missRect=TRUE, trend_period=trend_period, mean_period=mean_period, colorForce=TRUE, info_header=df_data, time_header=df_data, foot_note=TRUE, info_height=2.8, time_ratio=2, var_ratio=3, foot_height=1.25, df_shapefile=df_shapefile, figdir=figdir, filename_opt='', resources_path=resources_path, logo_dir=logo_dir, AEAGlogo_file=AEAGlogo_file, INRAElogo_file=INRAElogo_file, FRlogo_file=FRlogo_file) }