script.R

###### A MODIFIER ######


# 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
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 BH data is stored
BHfiledir = 
    # ""
    "BanqueHydro_Export2021"

## Manual selection ##
# Name of the file that will be analysed from the BH directory
BHfilename =
    # ""
    c("H5920011_HYDRO_QJM.txt")#, "K4470010_HYDRO_QJM.txt")
    # "all"

## Or list selection ##
# Path to the list file of BH data that will be analysed
BHlistdir = 
    ""

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

### NIVALE ###
# Path to the directory where NV data is stored
NVfiledir = 
    ""
    # "France207"

# Name of the file that will be analysed from the NV directory
NVfilename = 
    ""
    # "all"


# Path to the list file of information about station that will be analysed
NVlistdir =
    ""

NVlistname = 
    ""
    # "liste_bv_principaux_global.txt"


### TREND ANALYSIS ###
# Time period to analyse
period = c("1960-01-01","2019-12-31")
# period = c("1960-01-01","2020-01-01")



########################


# FILE STRUCTURE #
# Set working directory
setwd(computer_work_path)

# Sourcing R file
source('processing/extractBH.R')
source('processing/extractNV.R')
source('processing/format.R')
source('processing/analyse.R')
source('plotting/panel.R')

# Usefull library
library(StatsAnalysisTrend)
library(dplyr)

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


# BANQUE HYDRO #
# Get only the selected station from a list station file
if (BHlistname != ""){
    df_selec = get_selection(computer_data_path, 
                             BHlistdir,
                             BHlistname,
                             cnames=c('code',
                                      'station', 
                                      'BV_km2',
                                      'axe_principal_concerne',
                                      'longueur_serie',
                                      'commentaires',
                                      'choix'), 
                             c_num=c('BV_km2',
                                      'longueur_serie'))

    BHfilename = df_selec[df_selec$ok,]$filename
}

# Extract information about selected stations
df_info_BH = extractBH_info(computer_data_path, BHfiledir, BHfilename)

# Extract data about selected stations
df_data_BH = extractBH_data(computer_data_path, BHfiledir, BHfilename)


# NIVALE #
# Extract information about selected stations
df_info_NV = extractNVlist_info(computer_data_path, NVfiledir, NVlistdir, NVlistname)
### /!\ missing station info ###

# Extract data about selected stations
df_data_NV = extractNV_data(computer_data_path, NVfiledir, NVfilename)


# JOIN #
df_join = join(df_data_BH, df_data_NV, df_info_BH, df_info_NV)
df_data = df_join$data
df_info = df_join$info


# TIME PANEL #
# Plot time panel of debit by stations
# panel(df_data, df_info, figdir, "")
# panel(df_data, df_info, figdir, "", is_sqrt=TRUE)
### /!\ Removed 185 row(s) containing missing values (geom_path) -> remove NA ###


# ANALYSE #
# Compute gap parameters for stations
# df_lac = get_lacune(df_data, df_info)


# AVERAGE ANNUAL FLOW : QA #
### /!\ verify order conservation ###
df_QAlist = prepare(df_data, colnamegroup=c('code'))

df_QAEx = extract.Var(data.station=df_QAlist,
                      funct=mean,
                      timestep='year',
                      period=period,
                      pos.datetime=1,
                      na.rm=TRUE)
df_QATrend = Estimate.stats(data.extract=df_QAEx)
df_QATrend = clean(df_QATrend, df_QAlist)


# MONTHLY MINIMUM FLOW IN THE YEAR : QMNA #
df_QMNAlist = prepare(df_data, colnamegroup=c('code'))

### /!\ PLUS RAPIDE ###
# fMNA = function (X) {
#     # prendre un paquet de 1 ans et faire la moyenne par mois et retourner le minimum des debit 
#     dpm = length(X)/12
#     # print(dpm)
#     # print(length(X))
#     monthmean = c()
#     for (i in 1:12) {
#         id = round(dpm*(i-1)+1, 0)
#         iu = round(i*dpm, 0)
#         monthmean = append(monthmean, mean(X[id:iu], na.rm=TRUE))
#         # print(paste('start', id))
#         # print(paste('end', iu))
#         # print('')
#     }
#     # print(monthmean)
#     return (min(monthmean, na.rm=TRUE))
# }
# df_QMNAEx = extract.Var(data.station=df_QMNAlist,
#                         funct=fMNA,
#                         period=period,
#                         pos.datetime=1)#,
                        # na.rm=TRUE) ### /!\ PAS COMPRIS ###

df_QMNAEx = extract.Var(data.station=df_QMNAlist,
                        funct=mean,
                        period=period,
                        timestep='month',
                        pos.datetime=1,
                        na.rm=TRUE)

### /!\ NOM DE COLONNE PAS CONSERVER ###
df_QMNAlist = reprepare(df_QMNAEx, df_QMNAlist, colnamegroup=c('code'))

df_QMNAEx = extract.Var(data.station=df_QMNAlist,
                        funct=min,
                        period=period,
                        timestep='year',
                        pos.datetime=1,
                        na.rm=TRUE)

df_QMNATrend = Estimate.stats(data.extract=df_QMNAEx)
df_QMNATrend = clean(df_QMNATrend, df_QMNAlist)