script.R

# \\\
# Copyright 2021-2022 Louis Héraut*1
#
# *1   INRAE, France
#      louis.heraut@inrae.fr
#
# This file is part of ash R toolbox.
#
# ash R toolbox is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# ash R toolbox is distributed in the hope that it will be useful, but 
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with ash R toolbox.  If not, see <https://www.gnu.org/licenses/>.
# ///
#
#
# script.R
#
# Script file to manage the trend analysis of the Adour-Garonne basin.
# Performs the necessary calls to processing and plotting functions in
# order to realise the hydrologic trend analysis of stations according
# to the input parameters. The nearest area belove is where you need to
# write your prefer parameters for the analysis. See the 'README.txt'
# file for more information.


############## START OF REGION TO MODIFY (without risk) ##############
# Path to the data
computer_data_path = 
    "/home/louis/Documents/bouleau/INRAE/CDD_stationnarite/data"
    # "C:\\Users\\louis.heraut\\Documents\\CDD_stationnarite\\data"
  
# Work path (it needs to end with '/ASH' directory)
computer_work_path = 
    "/home/louis/Documents/bouleau/INRAE/CDD_stationnarite/ASH"
    # "C:\\Users\\louis.heraut\\Documents\\CDD_stationnarite\\ASH"


## BANQUE HYDRO
# Path to the directory where Banque Hydro (BH) data is stored
# from the work path
filedir = 
    # ""
    "BanqueHydro_Export2021"


## MANUAL SELECTION
# 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(
        "S2235610_HYDRO_QJM.txt",
        "P0885010_HYDRO_QJM.txt",
        "P0364010_HYDRO_QJM.txt",
        "O7635010_HYDRO_QJM.txt",
        "O3141010_HYDRO_QJM.txt",
        "Q6332510_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
# 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',
           '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'
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


## 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, 
                        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)'),
                 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)
}