Skip to content
GitLab
Select Git revision
  • master default
Find file BlameHistoryPermalink
script.R 13.66 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 
        # "P1712910_HYDRO_QJM.txt",
63 
        # "P0885010_HYDRO_QJM.txt",
64 
        # "O5055010_HYDRO_QJM.txt",
65 
        # "O0384010_HYDRO_QJM.txt",
66 
        # "S4214010_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" # "O8255010_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" ## TREND 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 yearLac_day = 3 # Sampling span of the data sampleSpan = c('05-01', '11-30') ## MAP # Is the hydrological network needs to be plot is_river = TRUE ############### 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))) {
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
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
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
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 ______________________________________________ # QA trend res = get_QAtrend(df_data, df_meta, period=trend_period, alpha=alpha,
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
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) # 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 _________________________________________________________ # 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)
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
# } # 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, is_river=is_river) ### 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) ### 5.2. Analysis layout _____________________________________________ datasheet_layout(toplot=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,
421422423424425426427428429430431
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)