script.R

# \\\
# Copyright 2021-2022 Louis Hraut*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 =
    ""

    # c(
      # "S2235610_HYDRO_QJM.txt", 
      # "P1712910_HYDRO_QJM.txt", 
      # "P0885010_HYDRO_QJM.txt",
      # "O5055010_HYDRO_QJM.txt",
      # "O0384010_HYDRO_QJM.txt"
      # )

    # c("S4214010_HYDRO_QJM.txt",
      # "O0384010_HYDRO_QJM.txt",
      # "Q7002910_HYDRO_QJM.txt",
      # "O0485110_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
AGlistdir = 
    ""

AGlistname = 
    # ""
    "Liste-station_RRSE.docx" 


## NIVALE SELECTION
# Path to the 'txt' list file of station from INRAE that will be analysed
INlistdir =
    ""

INlistname = 
    ""
    # "INRAE_selection.txt"


## TREND ANALYSIS
# Time period to analyse
periodAll = c("1800-01-01", "2019-12-31")
periodSub = c("1968-01-01", "2019-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")
mean_period = list(period1, period2)

# p value thresold
p_thresold = 0.1 #c(0.01, 0.05, 0.1)


## MAP
# 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 river shape from 'computer_data_path' 
rv_shpdir = 'map/river'
rv_shpname = 'CoursEau_FXX.shp'
############### 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='latin1')
source('processing/format.R', encoding='latin1')
source('processing/analyse.R', encoding='latin1')
source('plotting/layout.R', encoding='latin1')

# 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))


## 2. SELECTION OF STATION
# Initialization of null data frame if there is no data selected
df_data_AG = NULL
df_data_IN = NULL
df_meta_AG = NULL
df_meta_IN = NULL

### 2.1. Selection of the Agence de l'eau Adour-Garonne 
if (AGlistname != ""){
    
    # Get only the selected station from a list station file
    df_selec_AG = get_selection_AG(computer_data_path, 
                             AGlistdir,
                             AGlistname,
                             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_AG[df_selec_AG$ok,]$filename

    #####
    # filename = filename[(1):(10)]
    #####

    # Extract metadata about selected stations
    df_meta_AG = extract_meta(computer_data_path, filedir, filename)
    # Extract data about selected stations
    df_data_AG = extract_data(computer_data_path, filedir, filename)
}

### 2.2. INRAE selection 
if (INlistname != ""){
    
    # Get only the selected station from a list station file
    df_selec_IN = get_selection_IN(computer_data_path, 
                                   INlistdir,
                                   INlistname)

    # Get filenames of the selection
    filename = df_selec_IN[df_selec_IN$ok,]$filename

    #####
    # filename = filename[(1+20):(16+20)]
    #####

    # Extract metadata about selected stations
    df_meta_IN = extract_meta(computer_data_path, filedir, filename)
    # Extract data about selected stations
    df_data_IN = extract_data(computer_data_path, filedir, filename)
} 

### 2.3. Manual selection 
if (AGlistname == "" & INlistname == "") {
    # Extract metadata about selected stations
    df_meta_AG = extract_meta(computer_data_path, filedir, filename)
    # Extract data about selected stations
    df_data_AG = extract_data(computer_data_path, filedir, filename)
}

### 2.4. Data join
df_join = join(df_data_AG, df_data_IN, df_meta_AG, df_meta_IN)
df_data = df_join$data
df_meta = df_join$meta


## 3. ANALYSE
### 3.1. Compute gap parameters for stations
df_meta = get_lacune(df_data, df_meta)

### 3.2. Trend analysis
# QA trend
res_QAtrend = get_QAtrend(df_data, period=trend_period,
                          p_thresold=p_thresold)

# QMNA tend
res_QMNAtrend = get_QMNAtrend(df_data, period=trend_period,
                              p_thresold=p_thresold)

# VCN10 trend
res_VCN10trend = get_VCN10trend(df_data, df_meta, 
                                period=trend_period,
                                p_thresold=p_thresold)

### 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. PLOTTING
# Shapefile importation in order to it only once time
df_shapefile = ini_shapefile(computer_data_path, fr_shpdir, fr_shpname, bs_shpdir, bs_shpname, rv_shpdir, rv_shpname, riv=FALSE)

### 4.1. Simple time panel to criticize station data
# Plot time panel of debit by stations
# datasheet_layout(list(df_data, df_data),
                 # layout_matrix=c(1, 2),
                 # df_meta=df_meta,
                 # missRect=list(TRUE, TRUE), 
                 # type=list('Q', 'sqrt(Q)'), 
                 # info_header=TRUE,
                 # time_header=NULL,
                 # var_ratio=3,
                 # figdir=figdir,
                 # filename_opt='time')

### 4.2. Analysis layout 
datasheet_layout(isplot=c(
                     'datasheet',
                     'matrix',
                     'map'
                 ),
                 df_data=list(res_QAtrend$data, res_QMNAtrend$data,
                              res_VCN10trend$data), 
                 layout_matrix=c(1, 2, 3),
                 df_meta=df_meta,
                 df_trend=list(res_QAtrend$trend, res_QMNAtrend$trend,
                               res_VCN10trend$trend), 
                 type=list('QA', 'QMNA', 'VCN10'),
                 missRect=list(TRUE, TRUE, TRUE),
                 trend_period=trend_period,
                 mean_period=mean_period,
                 info_header=TRUE,
                 time_header=df_data,
                 info_ratio=2, 
                 time_ratio=2, 
                 var_ratio=3,
                 df_shapefile=df_shapefile,
                 figdir=figdir,
                 filename_opt='')