matrix.R 44.6 KB
Newer Older
Heraut Louis's avatar
Heraut Louis committed
# \\\
Heraut Louis's avatar
Heraut Louis committed
# Copyright 2021-2022 Louis Héraut*1
Heraut Louis's avatar
Heraut Louis committed
#
# *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/>.
# ///
#
#
# plotting/matrix.R
#
Heraut Louis's avatar
Heraut Louis committed
# Allows the creation of a summarizing matrix of trend and break analyses
Heraut Louis's avatar
Heraut Louis committed


Heraut Louis's avatar
Heraut Louis committed
## 1. MATRIX PANEL
Heraut Louis's avatar
Heraut Louis committed
# Generates a summarizing matrix of the trend analyses of all station for different hydrological variables and periods. Also shows difference of means between specific periods.
Heraut Louis's avatar
Heraut Louis committed
matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice=NULL, outdirTmp='', outnameTmp='matrix', title=NULL, A3=FALSE) {

Heraut Louis's avatar
Heraut Louis committed
    # Number of variable/plot
Heraut Louis's avatar
Heraut Louis committed
    nbp = length(list_df2plot)
    
    # Get all different stations code
Heraut Louis's avatar
Heraut Louis committed
    Code = levels(factor(df_meta$code))    
Heraut Louis's avatar
Heraut Louis committed
    nCode = length(Code)

Heraut Louis's avatar
Heraut Louis committed
    # Gets a trend example
Heraut Louis's avatar
Heraut Louis committed
    df_trend = list_df2plot[[1]]$trend

    # Convert 'trend_period' to list
    trend_period = as.list(trend_period)
    # Number of trend period
    nPeriod_trend = length(trend_period)
Heraut Louis's avatar
Heraut Louis committed

    # Fix the maximal number of period to the minimal possible
Heraut Louis's avatar
Heraut Louis committed
    nPeriod_max = 0
Heraut Louis's avatar
Heraut Louis committed
    # For all code 
Heraut Louis's avatar
Heraut Louis committed
    for (code in Code) {
Heraut Louis's avatar
Heraut Louis committed
        # Extracts the trend corresponding to the code
Heraut Louis's avatar
Heraut Louis committed
        df_trend_code = df_trend[df_trend$code == code,]
        
Heraut Louis's avatar
Heraut Louis committed
        # Extract start and end of trend periods
        Start = df_trend_code$period_start
Heraut Louis's avatar
Heraut Louis committed
        End = df_trend_code$period_end
Heraut Louis's avatar
Heraut Louis committed
        # Get the name of the different period
        UStart = levels(factor(Start))        
Heraut Louis's avatar
Heraut Louis committed
        UEnd = levels(factor(End))
Heraut Louis's avatar
Heraut Louis committed

        # Compute the max of different start and end
        # so the number of different period
Heraut Louis's avatar
Heraut Louis committed
        nPeriod = max(length(UStart), length(UEnd))

Heraut Louis's avatar
Heraut Louis committed
        # If the number of period for the trend is greater
        # than the current max period, stocks it
Heraut Louis's avatar
Heraut Louis committed
        if (nPeriod > nPeriod_max) {
            nPeriod_max = nPeriod
        }
    }

Heraut Louis's avatar
Heraut Louis committed
    tab_Start =  array(rep('', nCode*nbp*nPeriod_max),
                       dim=c(nCode, nbp, nPeriod_max))
    tab_End = array(rep('', nCode*nbp*nPeriod_max),
                    dim=c(nCode, nbp, nPeriod_max))
    tab_Code = array(rep('', nCode*nbp*nPeriod_max),
                     dim=c(nCode, nbp, nPeriod_max))
    tab_Periods = array(rep('', nCode*nbp*nPeriod_max),
                        dim=c(nCode, nbp, nPeriod_max))
    
    # For all code
    for (k in 1:nCode) {
Heraut Louis's avatar
Heraut Louis committed
        # Gets the code
Heraut Louis's avatar
Heraut Louis committed
        code = Code[k]

        for (i in 1:nbp) {
            df_trend = list_df2plot[[i]]$trend
            # Extracts the trend corresponding to the code
            df_trend_code = df_trend[df_trend$code == code,]
            
            # Extract start and end of trend periods
            Start = df_trend_code$period_start
            End = df_trend_code$period_end
            # Get the name of the different period
            UStart = levels(factor(Start))        
            UEnd = levels(factor(End))
            
            # Compute the max of different start and end
            # so the number of different period
            nPeriod = max(length(UStart), length(UEnd))

            # For all the period
            for (j in 1:nPeriod_max) {
                # Stocks period
                Periods = paste(Start[j],
                                End[j],
                                sep=' / ')
                
                tab_Start[k, i, j] = as.character(Start[j])
                tab_End[k, i, j] = as.character(End[j])
                tab_Code[k, i, j] = code
                tab_Periods[k, i, j] = Periods
                
            }
Heraut Louis's avatar
Heraut Louis committed
        }
    }
Heraut Louis's avatar
Heraut Louis committed

    
Heraut Louis's avatar
Heraut Louis committed
    
    # Blank array to store mean of the trend for each
    # station, perdiod and variable
Heraut Louis's avatar
Heraut Louis committed
    TrendMean_code = array(rep(1, nPeriod_trend*nbp*nCode),
                           dim=c(nPeriod_trend, nbp, nCode))
Heraut Louis's avatar
Heraut Louis committed
    # For all the trend period
Heraut Louis's avatar
Heraut Louis committed
    for (j in 1:nPeriod_trend) {
Heraut Louis's avatar
Heraut Louis committed
        # For all the code
Heraut Louis's avatar
Heraut Louis committed
        for (k in 1:nCode) {
Heraut Louis's avatar
Heraut Louis committed
            # Gets the code
Heraut Louis's avatar
Heraut Louis committed
            code = Code[k]
Heraut Louis's avatar
Heraut Louis committed
            # For all variable
Heraut Louis's avatar
Heraut Louis committed
            for (i in 1:nbp) {
Heraut Louis's avatar
Heraut Louis committed
                # Extracts the data corresponding to the
                # current variable
Heraut Louis's avatar
Heraut Louis committed
                df_data = list_df2plot[[i]]$data
Heraut Louis's avatar
Heraut Louis committed
                # Extracts the trend corresponding to the
                # current variable
Heraut Louis's avatar
Heraut Louis committed
                df_trend = list_df2plot[[i]]$trend
                p_threshold = list_df2plot[[i]]$p_threshold
Heraut Louis's avatar
Heraut Louis committed
                # Extracts the data corresponding to the code
                df_data_code = df_data[df_data$code == code,]
                # Extracts the trend corresponding to the code
Heraut Louis's avatar
Heraut Louis committed
                df_trend_code = df_trend[df_trend$code == code,]

Heraut Louis's avatar
Heraut Louis committed
                # Gets the associated time info
Heraut Louis's avatar
Heraut Louis committed
                Start = tab_Start[k, i, j]
                End = tab_End[k, i, j]
                Periods = tab_Periods[k, i, j]
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
                # Extracts the corresponding data for the period
Heraut Louis's avatar
Heraut Louis committed
                df_data_code_per =
                    df_data_code[df_data_code$Date >= Start 
                                 & df_data_code$Date <= End,]
Heraut Louis's avatar
Heraut Louis committed
                # Same for trend
Heraut Louis's avatar
Heraut Louis committed
                df_trend_code_per = 
                    df_trend_code[df_trend_code$period_start == Start 
                                  & df_trend_code$period_end == End,]

Heraut Louis's avatar
Heraut Louis committed
                # Computes the number of trend analysis selected
Heraut Louis's avatar
Heraut Louis committed
                Ntrend = nrow(df_trend_code_per)
Heraut Louis's avatar
Heraut Louis committed
                # If there is more than one trend on the same period
Heraut Louis's avatar
Heraut Louis committed
                if (Ntrend > 1) {
Heraut Louis's avatar
Heraut Louis committed
                    # Takes only the first because they are similar
Heraut Louis's avatar
Heraut Louis committed
                    df_trend_code_per = df_trend_code_per[1,]
                }
Heraut Louis's avatar
Heraut Louis committed

                # Computes the mean of the data on the period
Heraut Louis's avatar
Heraut Louis committed
                dataMean = mean(df_data_code_per$Value, na.rm=TRUE)
Heraut Louis's avatar
Heraut Louis committed
                # Normalises the trend value by the mean of the data
Heraut Louis's avatar
Heraut Louis committed
                trendMean = df_trend_code_per$trend / dataMean

Heraut Louis's avatar
Heraut Louis committed
                # If the p value is under the threshold
Heraut Louis's avatar
Heraut Louis committed
                if (df_trend_code_per$p <= p_threshold){
Heraut Louis's avatar
Heraut Louis committed
                    # Stores the averaged trend
Heraut Louis's avatar
Heraut Louis committed
                    TrendMean_code[j, i, k] = trendMean
Heraut Louis's avatar
Heraut Louis committed
                # Otherwise
Heraut Louis's avatar
Heraut Louis committed
                } else {
Heraut Louis's avatar
Heraut Louis committed
                    # Do not stocks it
Heraut Louis's avatar
Heraut Louis committed
                    TrendMean_code[j, i, k] = NA
                }
            }
        }
    }
Heraut Louis's avatar
Heraut Louis committed
    # Computes the min and the max of the mean trend for
    # all the station
Heraut Louis's avatar
Heraut Louis committed
    minTrendMean = apply(TrendMean_code, c(1, 2), min, na.rm=TRUE)
    maxTrendMean = apply(TrendMean_code, c(1, 2), max, na.rm=TRUE)

Heraut Louis's avatar
Heraut Louis committed
    # Blank vectors to store info about trend analyses
Heraut Louis's avatar
Heraut Louis committed
    Periods_trend = c()
    NPeriod_trend = c()
Heraut Louis's avatar
Heraut Louis committed
    Var_trend = c()
    Type_trend = c()
Heraut Louis's avatar
Heraut Louis committed
    Code_trend = c()
    Pthresold_trend = c()
    TrendMean_trend = c()
    DataMean_trend = c()
    Fill_trend = c()
    Color_trend = c()
Heraut Louis's avatar
Heraut Louis committed
    # For all the trend period
Heraut Louis's avatar
Heraut Louis committed
    for (j in 1:nPeriod_trend) {
Heraut Louis's avatar
Heraut Louis committed
        # For all code
Heraut Louis's avatar
Heraut Louis committed
        for (k in 1:nCode) {
            # Gets the code
            code = Code[k]
Heraut Louis's avatar
Heraut Louis committed
            # For all variable
Heraut Louis's avatar
Heraut Louis committed
            for (i in 1:nbp) {
Heraut Louis's avatar
Heraut Louis committed
                # Extracts the data corresponding to the current variable
Heraut Louis's avatar
Heraut Louis committed
                df_data = list_df2plot[[i]]$data
Heraut Louis's avatar
Heraut Louis committed
                # Extracts the trend corresponding to the
                # current variable
Heraut Louis's avatar
Heraut Louis committed
                df_trend = list_df2plot[[i]]$trend
                p_threshold = list_df2plot[[i]]$p_threshold
Heraut Louis's avatar
Heraut Louis committed
                # Extract the variable of the plot
Heraut Louis's avatar
Heraut Louis committed
                var = list_df2plot[[i]]$var
Heraut Louis's avatar
Heraut Louis committed
                # Extract the type of the variable to plot
                type = list_df2plot[[i]]$type
Heraut Louis's avatar
Heraut Louis committed
                # Extracts the data corresponding to the code
                df_data_code = df_data[df_data$code == code,]
                # Extracts the trend corresponding to the code
Heraut Louis's avatar
Heraut Louis committed
                df_trend_code = df_trend[df_trend$code == code,]

Heraut Louis's avatar
Heraut Louis committed
                # Gets the associated time info
Heraut Louis's avatar
Heraut Louis committed
                Start = tab_Start[k, i, j]
                End = tab_End[k, i, j]
                Periods = tab_Periods[k, i, j]
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
                # Extracts the corresponding data for the period
Heraut Louis's avatar
Heraut Louis committed
                df_data_code_per =
                    df_data_code[df_data_code$Date >= Start 
                                 & df_data_code$Date <= End,]
Heraut Louis's avatar
Heraut Louis committed
                # Same for trend
Heraut Louis's avatar
Heraut Louis committed
                df_trend_code_per = 
                    df_trend_code[df_trend_code$period_start == Start 
                                  & df_trend_code$period_end == End,]

Heraut Louis's avatar
Heraut Louis committed
                # Computes the number of trend analysis selected
Heraut Louis's avatar
Heraut Louis committed
                Ntrend = nrow(df_trend_code_per)
Heraut Louis's avatar
Heraut Louis committed
                # If there is more than one trend on the same period
Heraut Louis's avatar
Heraut Louis committed
                if (Ntrend > 1) {
Heraut Louis's avatar
Heraut Louis committed
                    # Takes only the first because they are similar
Heraut Louis's avatar
Heraut Louis committed
                    df_trend_code_per = df_trend_code_per[1,]
                }
Heraut Louis's avatar
Heraut Louis committed

                # Computes the mean of the data on the period
Heraut Louis's avatar
Heraut Louis committed
                dataMean = mean(df_data_code_per$Value, na.rm=TRUE)
Heraut Louis's avatar
Heraut Louis committed
                # Normalises the trend value by the mean of the data
Heraut Louis's avatar
Heraut Louis committed
                trendMean = df_trend_code_per$trend / dataMean

Heraut Louis's avatar
Heraut Louis committed
                # If the p value is under the threshold
Heraut Louis's avatar
Heraut Louis committed
                if (df_trend_code_per$p <= p_threshold){
Heraut Louis's avatar
Heraut Louis committed
                    # Gets the color associated to the averaged trend
Heraut Louis's avatar
Heraut Louis committed
                    color_res = get_color(trendMean, 
                                          minTrendMean[j, i],
                                          maxTrendMean[j, i],
                                          palette_name='perso',
                                          reverse=TRUE)
Heraut Louis's avatar
Heraut Louis committed
                    # Specifies the color fill and contour of
                    # table cells
Heraut Louis's avatar
Heraut Louis committed
                    fill = color_res
                    color = 'white'
                    Pthresold = p_thresold
Heraut Louis's avatar
Heraut Louis committed
                # Otherwise it is not significative
Heraut Louis's avatar
Heraut Louis committed
                } else { 
                    fill = 'white'
                    color = 'grey85'  
                    Pthresold = NA
                }

Heraut Louis's avatar
Heraut Louis committed
                # Stores info needed to plot
Heraut Louis's avatar
Heraut Louis committed
                Periods_trend = append(Periods_trend, Periods)
                NPeriod_trend = append(NPeriod_trend, j)
Heraut Louis's avatar
Heraut Louis committed
                Var_trend = append(Var_trend, var)
Heraut Louis's avatar
Heraut Louis committed
                Type_trend = append(Type_trend, type)
Heraut Louis's avatar
Heraut Louis committed
                Code_trend = append(Code_trend, code)
                Pthresold_trend = append(Pthresold_trend, Pthresold)
                TrendMean_trend = append(TrendMean_trend, trendMean)
                DataMean_trend = append(DataMean_trend, dataMean)
                Fill_trend = append(Fill_trend, fill)
                Color_trend = append(Color_trend, color)
            }
        }
    }

    # If there is a 'mean_period'
    if (!is.null(mean_period)) {
Heraut Louis's avatar
Heraut Louis committed
        # Blank vectors to store info about breaking analysis
Heraut Louis's avatar
Heraut Louis committed
        Periods_mean = c()
        NPeriod_mean = c()
Heraut Louis's avatar
Heraut Louis committed
        Var_mean = c()
        Type_mean = c()
Heraut Louis's avatar
Heraut Louis committed
        Code_mean = c()
        DataMean_mean = c()
        BreakMean_mean = c()
        
        # Convert 'mean_period' to list
        mean_period = as.list(mean_period)
        # Number of mean period
        nPeriod_mean = length(mean_period)

Heraut Louis's avatar
Heraut Louis committed
        # Blank array to store difference of mean between two periods
Heraut Louis's avatar
Heraut Louis committed
        BreakMean_code = array(rep(1, nPeriod_mean*nbp*nCode),
                               dim=c(nPeriod_mean, nbp, nCode))
Heraut Louis's avatar
Heraut Louis committed
        # Blank array to store mean for a temporary period in order
        # to compute the difference of mean with a second period
Heraut Louis's avatar
Heraut Louis committed
        dataMeantmp = array(rep(NA, nbp*nCode),
                            dim=c(nbp, nCode))
        
Heraut Louis's avatar
Heraut Louis committed
        # For all period of breaking analysis
Heraut Louis's avatar
Heraut Louis committed
        for (j in 1:nPeriod_mean) {
Heraut Louis's avatar
Heraut Louis committed
            # For all the code
Heraut Louis's avatar
Heraut Louis committed
            for (k in 1:nCode) {
Heraut Louis's avatar
Heraut Louis committed
                # Gets the code
Heraut Louis's avatar
Heraut Louis committed
                code = Code[k]
Heraut Louis's avatar
Heraut Louis committed
                # For all variable
Heraut Louis's avatar
Heraut Louis committed
                for (i in 1:nbp) {
Heraut Louis's avatar
Heraut Louis committed
                    # Extracts the data corresponding to
                    # the current variable
Heraut Louis's avatar
Heraut Louis committed
                    df_data = list_df2plot[[i]]$data
Heraut Louis's avatar
Heraut Louis committed
                    # Extract the variable of the plot
Heraut Louis's avatar
Heraut Louis committed
                    var = list_df2plot[[i]]$var
Heraut Louis's avatar
Heraut Louis committed
                    # Extract the type of the variable to plot
                    type = list_df2plot[[i]]$type
Heraut Louis's avatar
Heraut Louis committed
                    # Extracts the data corresponding to the code
Heraut Louis's avatar
Heraut Louis committed
                    df_data_code = df_data[df_data$code == code,] 
                    
                    # Get the current start and end of the sub period
                    Start_mean = mean_period[[j]][1]
                    End_mean = mean_period[[j]][2]
                    
                    # Extract the data corresponding to this sub period
                    df_data_code_per =
                        df_data_code[df_data_code$Date >= Start_mean 
                                     & df_data_code$Date <= End_mean,]
                    
Heraut Louis's avatar
Heraut Louis committed
                    # Min max for the sub period
Heraut Louis's avatar
Heraut Louis committed
                    Datemin = min(df_data_code_per$Date)
                    Datemax = max(df_data_code_per$Date)
Heraut Louis's avatar
Heraut Louis committed
                    # Creates a period name
Heraut Louis's avatar
Heraut Louis committed
                    Periods = paste(Datemin, Datemax,
                                    sep=' / ')

                    # Mean of the flow over the sub period
Heraut Louis's avatar
Heraut Louis committed
                    dataMean = mean(df_data_code_per$Value,
Heraut Louis's avatar
Heraut Louis committed
                                    na.rm=TRUE)

Heraut Louis's avatar
Heraut Louis committed
                    # If this in not the first period
Heraut Louis's avatar
Heraut Louis committed
                    if (j > 1) {
Heraut Louis's avatar
Heraut Louis committed
                        # Compute the difference of mean
Heraut Louis's avatar
Heraut Louis committed
                        Break = dataMean - dataMeantmp[i, k]
Heraut Louis's avatar
Heraut Louis committed
                    # Otherwise for the first period
Heraut Louis's avatar
Heraut Louis committed
                    } else {
Heraut Louis's avatar
Heraut Louis committed
                        # Stocks NA
Heraut Louis's avatar
Heraut Louis committed
                        Break = NA
                    }
Heraut Louis's avatar
Heraut Louis committed
                    # Normalises the break by the mean of the
                    # initial period
Heraut Louis's avatar
Heraut Louis committed
                    BreakMean = Break / dataMeantmp[i, k]
Heraut Louis's avatar
Heraut Louis committed
                    # Stores the result
Heraut Louis's avatar
Heraut Louis committed
                    BreakMean_code[j, i, k] = BreakMean
Heraut Louis's avatar
Heraut Louis committed
                    # Stores temporarily the mean of the current period
Heraut Louis's avatar
Heraut Louis committed
                    dataMeantmp[i, k] = dataMean
                    
Heraut Louis's avatar
Heraut Louis committed
                    # Stores info needed to plot
Heraut Louis's avatar
Heraut Louis committed
                    Periods_mean = append(Periods_mean, Periods)
                    NPeriod_mean = append(NPeriod_mean, j)
Heraut Louis's avatar
Heraut Louis committed
                    Var_mean = append(Var_mean, var)
Heraut Louis's avatar
Heraut Louis committed
                    Type_mean = append(Type_mean, type)
Heraut Louis's avatar
Heraut Louis committed
                    Code_mean = append(Code_mean, code)
                    DataMean_mean = append(DataMean_mean, dataMean)
                    BreakMean_mean = append(BreakMean_mean,
                                            BreakMean)
                }
            }
        }
Heraut Louis's avatar
Heraut Louis committed
        # Computes the min and the max of the averaged trend for
Heraut Louis's avatar
Heraut Louis committed
        # all the station
Heraut Louis's avatar
Heraut Louis committed
        minBreakMean = apply(BreakMean_code, c(1, 2),
                             min, na.rm=TRUE)
        maxBreakMean = apply(BreakMean_code, c(1, 2),
                             max, na.rm=TRUE)
Heraut Louis's avatar
Heraut Louis committed
        # Blanks vector to store color info
Heraut Louis's avatar
Heraut Louis committed
        Fill_mean = c()
        Color_mean = c()
Heraut Louis's avatar
Heraut Louis committed
        # Index to count over all break computed
Heraut Louis's avatar
Heraut Louis committed
        ii = 1
        for (j in 1:nPeriod_mean) {
Heraut Louis's avatar
Heraut Louis committed
            # For all the code
Heraut Louis's avatar
Heraut Louis committed
            for (k in 1:nCode) {
Heraut Louis's avatar
Heraut Louis committed
                # Gets the code
Heraut Louis's avatar
Heraut Louis committed
                code = Code[k]
Heraut Louis's avatar
Heraut Louis committed
                # For all variable
Heraut Louis's avatar
Heraut Louis committed
                for (i in 1:nbp) {
Heraut Louis's avatar
Heraut Louis committed
                    # Extracts averaged breaking
Heraut Louis's avatar
Heraut Louis committed
                    BreakMean = BreakMean_mean[ii]
Heraut Louis's avatar
Heraut Louis committed
                    # Gets the color associated
Heraut Louis's avatar
Heraut Louis committed
                    color_res = get_color(BreakMean, 
                                          minBreakMean[j, i],
                                          maxBreakMean[j, i],
                                          palette_name='perso',
                                          reverse=TRUE)
Heraut Louis's avatar
Heraut Louis committed
                    # Gets the fill and contour color
Heraut Louis's avatar
Heraut Louis committed
                    fill = color_res
                    color = 'white'
Heraut Louis's avatar
Heraut Louis committed
                    # Stores it
Heraut Louis's avatar
Heraut Louis committed
                    Fill_mean = append(Fill_mean, fill)
                    Color_mean = append(Color_mean, color)
Heraut Louis's avatar
Heraut Louis committed
                    # Passes to the next index
Heraut Louis's avatar
Heraut Louis committed
                    ii = ii + 1
                }
            }
        }
    }

Heraut Louis's avatar
Heraut Louis committed
    # If the slice option is not specified, the info for all
    # stations will be draw on the same page 
Heraut Louis's avatar
Heraut Louis committed
    if (is.null(slice)) {
        slice = nCode
    }

Heraut Louis's avatar
Heraut Louis committed
    # Gets all the different type of plots
    Type = levels(factor(Type_trend))
    nbType = length(Type)
    # For all the type of plots
    for (itype in 1:nbType) {
        # Gets the type
        type = Type[itype]
        
        # Extracts each possibilities of first letter of station code
        firstLetter = levels(factor(substr(Code, 1, 1)))
        # Number of different first letters
        nfL = length(firstLetter)
        # For all the available first letter
        for (ifL in 1:nfL) {
            # Gets the first letter
            fL = firstLetter[ifL]

            # Get only station code with the same first letter 
            subCodefL = Code[substr(Code, 1, 1) == fL]
            # Counts the number of station in it
            nsubCodefL = length(subCodefL)
            # Computes the number of pages needed to plot all stations
            nMat = as.integer(nsubCodefL/slice) + 1
            # For all the pages
            for (iMat in 1:nMat) {
                # Print the matrix name
                print(paste('Matrix ', iMat, '/', nMat,
                            ' of ', type,
                            ' for region : ', fL,
                            "   (",
                            round((ifL + nfL*(itype-1)) / (nfL*2) * 100,
                                  1),
                            " %)", 
                            sep=''))
Heraut Louis's avatar
Heraut Louis committed
                # Extracts the station for the current page
                subCode = subCodefL[(slice*(iMat-1)+1):(slice*iMat)]
                # Removes NA stations
                subCode = subCode[!is.na(subCode)]
                # Reverses verticale order of stations
                subCode = rev(subCode)
                # Gets the number of station for the page
                nsubCode = length(subCode)

                # Creates logical vector to select only info about
                # stations that will be plot on the page
                CodefL_trend =
                    Code_trend %in% subCode & Type_trend == type
                # Extracts those info
                subPeriods_trend = Periods_trend[CodefL_trend]
                subNPeriod_trend = NPeriod_trend[CodefL_trend]
                subVar_trend = Var_trend[CodefL_trend]
                subType_trend = Type_trend[CodefL_trend]
                subCode_trend = Code_trend[CodefL_trend]
                subPthresold_trend = Pthresold_trend[CodefL_trend]
                subTrendMean_trend = TrendMean_trend[CodefL_trend]
                subDataMean_trend = DataMean_trend[CodefL_trend]
                subFill_trend = Fill_trend[CodefL_trend]
                subColor_trend = Color_trend[CodefL_trend]

                # Same for breaking analysis
                CodefL_mean =
                    Code_mean %in% subCode & Type_mean == type
                # Extracts right info
                subPeriods_mean = Periods_mean[CodefL_mean]
                subNPeriod_mean = NPeriod_mean[CodefL_mean]
                subVar_mean = Var_mean[CodefL_mean]
                subType_mean = Type_mean[CodefL_mean]
                subCode_mean = Code_mean[CodefL_mean]
                subDataMean_mean = DataMean_mean[CodefL_mean]
                subBreakMean_mean = BreakMean_mean[CodefL_mean]
                subFill_mean = Fill_mean[CodefL_mean]
                subColor_mean = Color_mean[CodefL_mean]
                
                # Gets the number of variable to plot in
                # function of the current type
                nbpMod =
                    length(levels(factor(subVar_trend)))                
                
                ### Plot ###
                # Fixes the height and width of the table according to
                # the number of station and the number of column to draw
                height = nsubCode
                width = nbpMod * 2 * nPeriod_trend + nPeriod_trend + nPeriod_mean * nbpMod + nPeriod_mean + nbpMod        

                # Fixes the size of the plot area to keep proportion right
                options(repr.plot.width=width, repr.plot.height=height)

                # Open a new plot with a personalise theme
                mat = ggplot() + theme_ash + 
                    # Modification of theme in order to remove axis
                    theme(
                        panel.border=element_blank(),
                        axis.text.x=element_blank(),
                        axis.text.y=element_blank(),
                        axis.ticks.y=element_blank(),
                        axis.ticks.x=element_blank(),
                        axis.title.y=element_blank(),
                        plot.margin=margin(t=5, r=5, b=5, l=5, unit="mm")
                    )

                # Extracts the name of the currently hydrological
                # region plotted
                title = df_meta[df_meta$code == subCode[1],]$region_hydro

                subtitle = paste(' ', iMat, '/', nMat,' ', type,
                                 sep='')    
                # Postion and name of the title
                xt = 1 - 6
                yt = height + 2
                Title = bquote(bold(.(title))[.(subtitle)])
                # Writes the title
Heraut Louis's avatar
Heraut Louis committed
                mat = mat +
Heraut Louis's avatar
Heraut Louis committed
                    annotate("text", x=xt, y=yt,
                             label=Title,
                             hjust=0, vjust=1, 
                             size=6, color="#00A3A8")

                ### Trend ###
                # For all the trend period
                for (j in 1:nPeriod_trend) {
                    # Extracts the info to plot associated to the
                    # right period
                    Periods_trend_per =
                        subPeriods_trend[subNPeriod_trend == j]
                    NPeriods_trend_per =
                        subNPeriod_trend[subNPeriod_trend == j]
                    Var_trend_per =
                        subVar_trend[subNPeriod_trend == j]
                    Type_trend_per =
                        subType_trend[subNPeriod_trend == j]
                    Code_trend_per =
                        subCode_trend[subNPeriod_trend == j]
                    Pthresold_trend_per =
                        subPthresold_trend[subNPeriod_trend == j]
                    TrendMean_trend_per =
                        subTrendMean_trend[subNPeriod_trend == j]
                    DataMean_trend_per =
                        subDataMean_trend[subNPeriod_trend == j]
                    Fill_trend_per =
                        subFill_trend[subNPeriod_trend == j]
                    Color_trend_per =
                        subColor_trend[subNPeriod_trend == j]

                    # Converts the variable list into levels for factor
                    levels = unlist(Var_trend_per)
                    # Converts the vector of hydrological variable to
                    # a vector of integer associated to those variable
                    Xtmp = as.integer(factor(as.character(Var_trend_per),
                                             levels=levels))

                    # Computes X position of the column for
                    # the period dates
                    Xc = j + (j - 1)*nbpMod*2
                    # Computes X positions of columns for
                    # the mean of variables
                    Xm = Xtmp + (j - 1)*nbpMod*2 + j
                    # Computes X positions of columns for
                    # the averaged trend
                    X = Xtmp + (j - 1)*nbpMod*2 + nbpMod + j
Heraut Louis's avatar
Heraut Louis committed
                    
Heraut Louis's avatar
Heraut Louis committed
                    # Computes Y positions of each line for each station
                    Y = as.integer(factor(Code_trend_per))
                    # Reverses vertical order of stations
                    Y = rev(Y)

                    # Position of a line to delimite periods
                    x = Xc - 0.4
                    xend = X[length(X)] + 0.4
                    y = height + 1.1
                    yend = height + 1.1
                    # Drawing of the line
                    mat = mat +
                        annotate("segment",
                                 x=x, xend=xend,
                                 y=y, yend=yend, 
                                 color="grey40", size=0.35)
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
                    # Position of the name of the current period
                    yt = y + 0.15
                    Start = trend_period[[j]][1]
                    End = trend_period[[j]][2]
                    # Name of the period
                    periodName =
                        bquote(bold('Période')~bold(.(as.character(j))))
                    # Naming the period
Heraut Louis's avatar
Heraut Louis committed
                    mat = mat +
Heraut Louis's avatar
Heraut Louis committed
                        annotate("text", x=x, y=yt,
                                 label=periodName,
                                 hjust=0, vjust=0.5, 
Heraut Louis's avatar
Heraut Louis committed
                                 size=3, color='grey40')

Heraut Louis's avatar
Heraut Louis committed
                    # For all the variable
                    for (i in 1:length(X)) {
                        mat = mat +
                            # Plots circles for averaged trends
                            gg_circle(r=0.45, xc=X[i], yc=Y[i],
                                      fill=Fill_trend_per[i],
                                      color=Color_trend_per[i]) +
                            # Plots circles for averaged of variables
                            gg_circle(r=0.45, xc=Xm[i], yc=Y[i],
                                      fill='white', color='grey40') +
                            # Plots circles for the column of period dates
                            gg_circle(r=0.45, xc=Xc, yc=Y[i],
                                      fill='white', color='grey40') 
                    }
                    # For all averaged trends on this periods
                    for (i in 1:length(TrendMean_trend_per)) {
                        # Extracts the value of the averaged trend
                        trendMean = TrendMean_trend_per[i]
                        # Converts it to the right format with two
                        # significant figures
                        trendMeanC = signif(trendMean*100, 2)

                        # If it is significative
                        if (!is.na(Pthresold_trend_per[i])) {
                            # The text color is white
                            Tcolor = 'white'
                            # Otherwise
                        } else {
                            # The text is grey
                            Tcolor = 'grey85'
                        }
                        
                        # Same for averaged variables over
                        # the current period
                        dataMean = DataMean_trend_per[i]
                        dataMeanC = signif(dataMean, 2)

                        mat = mat +
                            # Writes the mean trend
                            annotate('text', x=X[i], y=Y[i],
                                     label=trendMeanC,
                                     hjust=0.5, vjust=0.5, 
                                     size=3, color=Tcolor) + 
                            # Writes the mean of the associated variable
                            annotate('text', x=Xm[i], y=Y[i],
                                     label=dataMeanC,
                                     hjust=0.5, vjust=0.5, 
                                     size=3, color='grey40')
                    }
Heraut Louis's avatar
Heraut Louis committed
                    
Heraut Louis's avatar
Heraut Louis committed
                    # Writes a name for the period dates column
Heraut Louis's avatar
Heraut Louis committed
                    mat = mat +
Heraut Louis's avatar
Heraut Louis committed
                        annotate('text', x=Xc, y=max(Y) + 0.9,
                                 label=bquote(bold('Début')),
Heraut Louis's avatar
Heraut Louis committed
                                 hjust=0.5, vjust=0.5, 
Heraut Louis's avatar
Heraut Louis committed
                                 size=3, color='grey20') + 
                        annotate('text', x=Xc, y=max(Y) + 0.63,
                                 label=bquote(bold('Fin')),
Heraut Louis's avatar
Heraut Louis committed
                                 hjust=0.5, vjust=0.5, 
Heraut Louis's avatar
Heraut Louis committed
                                 size=3, color='grey20')
                    
                    # For all variable
                    for (i in 1:nbpMod) {
                        # Extract the variable of the plot
                        var = subVar_trend[i]
                        mat = mat +
                            # Writes the unit of the variable
                            annotate('text', x=X[i], y=max(Y) + 0.63,
                                     label=bquote('[%.'*ans^{-1}*']'),
                                     hjust=0.5, vjust=0.5, 
                                     size=2, color='grey40') +
                            # Writes the type of the variable
                            annotate('text', x=X[i], y=max(Y) + 0.9,
                                     label=bquote(.(var)),
                                     hjust=0.5, vjust=0.5, 
                                     size=3.25, color='grey20') +
                            # Writes the unit of the averaged variable
                            annotate('text', x=Xm[i], y=max(Y) + 0.63,
                                     label=bquote('['*m^3*'.'*s^{-1}*']'),
                                     hjust=0.5, vjust=0.5, 
                                     size=2, color='grey40') +
                            # Writes the type of the averaged variable
                            annotate('text', x=Xm[i], y=max(Y) + 0.9,
                                     label=expr(bar(!!var)),
                                     hjust=0.5, vjust=0.5, 
                                     size=3.25, color='grey20')
                    }

                    # For all the station on the page
                    for (k in 1:nsubCode) {
                        # Gets the code
                        code = subCode[k]
                        # Extracts label for the period dates
                        label =
                            Periods_trend_per[Code_trend_per == code][1]
                        # Gets the start and end of the period
                        # for the station
                        periodStart = substr(label, 1, 4)
                        periodEnd = substr(label, 14, 17)
                        
                        mat = mat +
                            # Writes the starting value
                            annotate('text', x=Xc, y=k + 0.13,
                                     label=bquote(bold(.(periodStart))),
                                     hjust=0.5, vjust=0.5, 
                                     size=3, color='grey40') + 
                            # Writes the ending value
                            annotate('text', x=Xc, y=k - 0.13,
                                     label=bquote(bold(.(periodEnd))),
                                     hjust=0.5, vjust=0.5, 
                                     size=3, color='grey40')
                    }
Heraut Louis's avatar
Heraut Louis committed
                }

Heraut Louis's avatar
Heraut Louis committed
                ### Mean ###
                # For all the trend period
                for (j in 1:nPeriod_mean) {
                    # Extracts the info to plot associated to the
                    # right period
                    Periods_mean_per =
                        subPeriods_mean[subNPeriod_mean == j]
                    NPeriods_mean_per =
                        subNPeriod_mean[subNPeriod_mean == j]
                    Var_mean_per =
                        subVar_mean[subNPeriod_mean == j]
                    Type_mean_per =
                        subType_mean[subNPeriod_mean == j]
                    Code_mean_per =
                        subCode_mean[subNPeriod_mean == j]
                    DataMean_mean_per =
                        subDataMean_mean[subNPeriod_mean == j]
                    BreakMean_mean_per =
                        subBreakMean_mean[subNPeriod_mean == j]
                    Fill_mean_per =
                        subFill_mean[subNPeriod_mean == j]
                    Color_mean_per =
                        subColor_mean[subNPeriod_mean == j]

                    # Converts the variable list into levels for factor
                    levels = unlist(Var_mean_per)
                    # Converts the vector of hydrological variable to
                    # a vector of integer associated to those variable
                    Xtmp_mean =
                        as.integer(factor(as.character(Var_mean_per),
                                          levels=levels))
                    # Computes X position of the column for
                    # the period dates
                    Xc_mean = j + (j - 1)*nbpMod + X[length(X)]
                    # Computes X positions of columns for
                    # the mean of variables
                    Xm_mean =
                        Xtmp_mean + (j - 1)*nbpMod + j + X[length(X)]
                    # Computes X positions of columns for
                    # the difference of mean between periods (break)
                    Xr_mean =
                        Xtmp_mean + (j - 1)*nbpMod*2 + j + X[length(X)]

                    # Computes Y positions of each line for each station
                    Y_mean = as.integer(factor(Code_mean_per))
                    # Reverses vertical order of stations
                    Y_mean = rev(Y_mean)
                    
                    # Position of a line to delimite periods
                    x = Xc_mean - 0.4
                    xend = Xm_mean[length(Xm_mean)] + 0.25
                    y = height + 1.1
                    yend = height + 1.1
Heraut Louis's avatar
Heraut Louis committed
                    # Drawing of the line
Heraut Louis's avatar
Heraut Louis committed
                    mat = mat +
                        annotate("segment",
                                 x=x, xend=xend,
                                 y=y, yend=yend, 
                                 color="grey40", size=0.35)
Heraut Louis's avatar
Heraut Louis committed

                    # Position of the name of the current period
                    yt = y + 0.15
                    Start = mean_period[[j]][1]
                    End = mean_period[[j]][2]
                    # Name of the period
                    periodName = bquote(bold('Période')~bold(.(as.character(j+nPeriod_trend))))
                    # Naming the period
Heraut Louis's avatar
Heraut Louis committed
                    mat = mat +
                        annotate("text", x=x, y=yt,
Heraut Louis's avatar
Heraut Louis committed
                                 label=periodName,
Heraut Louis's avatar
Heraut Louis committed
                                 hjust=0, vjust=0.5, 
                                 size=3, color='grey40')
Heraut Louis's avatar
Heraut Louis committed

                    # If this is not the first period
Heraut Louis's avatar
Heraut Louis committed
                    if (j > 1) {
Heraut Louis's avatar
Heraut Louis committed
                        # Position of a line to delimite results of
                        # difference of mean bewteen periods
                        x = Xr_mean[1] - 0.4
                        xend = Xr_mean[length(Xr_mean)] + 0.25
                        # Drawing of the line
Heraut Louis's avatar
Heraut Louis committed
                        mat = mat +
Heraut Louis's avatar
Heraut Louis committed
                            annotate("segment",
                                     x=x, xend=xend,
                                     y=y, yend=yend, 
                                     color="grey40", size=0.35)
                        # Naming the breaking columns
                        breakName =  bquote(bold('Écart')~bold(.(as.character(j-1+nPeriod_trend)))*bold('-')*bold(.(as.character(j+nPeriod_trend))))
                        # Writes the name
                        mat = mat +
                            annotate("text", x=x, y=yt,
                                     label=breakName,
                                     hjust=0, vjust=0.5, 
                                     size=3, color='grey40')
Heraut Louis's avatar
Heraut Louis committed
                    }

Heraut Louis's avatar
Heraut Louis committed
                    # For all the variable
                    for (i in 1:length(Xm_mean)) {
                        mat = mat +
                            # Plots circles for averaged variables
                            gg_circle(r=0.45, xc=Xm_mean[i], yc=Y[i],
                                      fill='white', color='grey40') +
                            # Plots circles for the column of period dates
                            gg_circle(r=0.45, xc=Xc_mean, yc=Y[i],
                                      fill='white', color='grey40')
                        
                        # If this is not the first period
                        if (j > 1) {
                            mat = mat +
                                # Plots circles for breaking results
                                gg_circle(r=0.45, xc=Xr_mean[i], yc=Y[i],
                                          fill=Fill_mean_per[i],
                                          color=Color_mean_per[i])
                        }
                    }

                    # For all averaged variables on this period
                    for (i in 1:length(DataMean_mean_per)) {
                        # Extracts values of averaged variables
                        dataMean = DataMean_mean_per[i]
Heraut Louis's avatar
Heraut Louis committed
                        # Converts it to the right format with two
                        # significant figures
Heraut Louis's avatar
Heraut Louis committed
                        dataMeanC = signif(dataMean, 2)
                        # Writes averaged variables values
Heraut Louis's avatar
Heraut Louis committed
                        mat = mat +
Heraut Louis's avatar
Heraut Louis committed
                            annotate('text', x=Xm_mean[i], y=Y[i],
                                     label=dataMeanC,
Heraut Louis's avatar
Heraut Louis committed
                                     hjust=0.5, vjust=0.5, 
Heraut Louis's avatar
Heraut Louis committed
                                     size=3, color='grey40')
                        # If this is not the first period
                        if (j > 1) {
                            # Extracts values of breaking between periods
                            BreakMean = BreakMean_mean_per[i]
                            # Converts it to the right format with two
                            # significant figures
                            BreakMeanC = signif(BreakMean*100, 2)
                            # Writes breaking values
                            mat = mat +
                                annotate('text', x=Xr_mean[i], y=Y[i],
                                         label=BreakMeanC,
                                         hjust=0.5, vjust=0.5, 
                                         size=3, color='white')   
                        }
Heraut Louis's avatar
Heraut Louis committed
                    }
Heraut Louis's avatar
Heraut Louis committed
                    # Writes a name for the period dates column
Heraut Louis's avatar
Heraut Louis committed
                    mat = mat +
Heraut Louis's avatar
Heraut Louis committed
                        annotate('text', x=Xc_mean, y=max(Y) + 0.9,
                                 label=bquote(bold('Début')),
Heraut Louis's avatar
Heraut Louis committed
                                 hjust=0.5, vjust=0.5, 
Heraut Louis's avatar
Heraut Louis committed
                                 size=3, color='grey20') + 
                        annotate('text', x=Xc_mean, y=max(Y) + 0.63,
                                 label=bquote(bold('Fin')),
Heraut Louis's avatar
Heraut Louis committed
                                 hjust=0.5, vjust=0.5, 
Heraut Louis's avatar
Heraut Louis committed
                                 size=3, color='grey20')
                    
                    # For all variables
                    for (i in 1:nbpMod) {
                        # Extract the variable of the plot
                        var = subVar_mean[i]
Heraut Louis's avatar
Heraut Louis committed
                        mat = mat +
Heraut Louis's avatar
Heraut Louis committed
                            # Writes the unit of the averaged variable
                            annotate('text',
                                     x=Xm_mean[i], y=max(Y) + 0.63,
                                     label=bquote('['*m^3*'.'*s^{-1}*']'),
Heraut Louis's avatar
Heraut Louis committed
                                     hjust=0.5, vjust=0.5, 
Heraut Louis's avatar
Heraut Louis committed
                                     size=2, color='grey40') +
Heraut Louis's avatar
Heraut Louis committed
                            # Writes the type of the averaged variable
                            annotate('text',
                                     x=Xm_mean[i], y=max(Y) + 0.9,
                                     label=expr(bar(!!var)),
                                     hjust=0.5, vjust=0.5, 
Heraut Louis's avatar
Heraut Louis committed
                                     size=3.25, color='grey20')
Heraut Louis's avatar
Heraut Louis committed

                        # If this is not the first period
                        if (j > 1) {
                            mat = mat +
                                # Writes the unit of the breaking variable
                                annotate('text', x=Xr_mean[i],
                                         y=max(Y) + 0.63,
                                         label=bquote('[%]'),
                                         hjust=0.5, vjust=0.5, 
                                         size=2, color='grey40') +
                                # Writes the type of the breaking variable
                                annotate('text', x=Xr_mean[i],
                                         y=max(Y) + 0.9,
                                         label=paste("d", var, sep=''),
                                         hjust=0.5, vjust=0.5,
                                         size=3.25, color='grey20')
                        }
Heraut Louis's avatar
Heraut Louis committed
                    }

Heraut Louis's avatar
Heraut Louis committed
                    # For all the station on the page
                    for (k in 1:nsubCode) {
                        # Gets the code
                        code = subCode[k]
                        # Extracts label for the period dates
                        label = Periods_mean_per[Code_mean_per == code][1]
                        # Gets the start and end of the period
                        # for the station
                        periodStart = substr(label, 1, 4)
                        periodEnd = substr(label, 14, 17)

                        mat = mat +
                            # # Writes the starting value
                            annotate('text', x=Xc_mean, y=k + 0.13,
                                     label=bquote(bold(.(periodStart))),
                                     hjust=0.5, vjust=0.5, 
                                     size=3, color='grey40') + 
                            # Writes the ending value
                            annotate('text', x=Xc_mean, y=k - 0.13,
                                     label=bquote(bold(.(periodEnd))),
                                     hjust=0.5, vjust=0.5, 
                                     size=3, color='grey40')
                    }            
                }
                
                ### Code ###
                # For all the station
Heraut Louis's avatar
Heraut Louis committed
                for (k in 1:nsubCode) {
Heraut Louis's avatar
Heraut Louis committed
                    # Gets the code
Heraut Louis's avatar
Heraut Louis committed
                    code = subCode[k]
Heraut Louis's avatar
Heraut Louis committed
                    # Gets the name of the station
                    name = df_meta[df_meta$code == code,]$nom
                    # Fixes a limit for the max number
                    # of characters available
                    ncharMax = 38
                    # If the number of character of the name is greater
                    # than the limit
                    if (nchar(name) > ncharMax) {
                        # Cuts the name and add '...'
                        name = paste(substr(name, 1, ncharMax),
                                     '...', sep='')
                    }
Heraut Louis's avatar
Heraut Louis committed

                    mat = mat +
Heraut Louis's avatar
Heraut Louis committed
                        # Writes the code of the station
                        annotate('text', x=0.3, y=k + 0.14,
                                 label=bquote(bold(.(code))),
                                 hjust=1, vjust=0.5, 
                                 size=3.5, color="#00A3A8") +
                        # Writes the name of the station
                        annotate('text', x=0.3, y=k - 0.14,
                                 label=name,
                                 hjust=1, vjust=0.5, 
                                 size=3.5, color="#00A3A8")
Heraut Louis's avatar
Heraut Louis committed
                }

Heraut Louis's avatar
Heraut Louis committed
                ### Environment ###
Heraut Louis's avatar
Heraut Louis committed
                mat = mat +
Heraut Louis's avatar
Heraut Louis committed
                    # Fixed coordinate system
                    coord_fixed() +
                    # X axis
                    scale_x_continuous(limits=c(1 - rel(6), 
                                                width + rel(0.5)),
                                       expand=c(0, 0)) + 
                    # Y axis
                    scale_y_continuous(limits=c(1 - rel(0.5), 
                                                height + rel(2)),
                                       expand=c(0, 0))
                
                # Paper format in A3 if needed
                if (A3) {
                    width = 42
                    height = 29.7
                    dpi = 300
                    # Otherwise in A4
                } else {
                    width = 29.7
                    height = 21
                    dpi = 100
                }
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
                # Saving
                ggsave(plot=mat,