map.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/>.
# ///
#
#
# plotting/map.R
#
# Deals with the creation of a map for presenting the trend analysis of hydrological variables


## 1. MAP PANEL
# Generates a map plot of the tendancy of a hydrological variable
map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='', codeLight=NULL, margin=NULL, showSea=TRUE, verbose=TRUE) {

    # Extract shapefiles
    df_france = df_shapefile$france
    df_bassin = df_shapefile$bassin
    df_subbassin = df_shapefile$subbassin
    df_river = df_shapefile$river

    # Number of variable/plot
    nbp = length(list_df2plot)
  
    # Get all different stations code
    Code = levels(factor(df_meta$code))
    nCode = length(Code)

    # Gets a trend example
    df_trend = list_df2plot[[1]]$trend
    
    nPeriod_max = 0
    for (code in Code) {
        # 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))

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

    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) {
        # Gets the code
        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
                
            }
        }
    }
    
    
    # Blank array to store mean of the trend for each
    # station, perdiod and variable
    TrendMean_code = array(rep(1, nPeriod_max*nbp*nCode),
                           dim=c(nPeriod_max, nbp, nCode))
    # For all the period
    for (j in 1:nPeriod_max) {
        # For all the code
        for (k in 1:nCode) {
            # Gets the code
            code = Code[k]
            # For all variable
            for (i in 1:nbp) {
                # Extracts the data corresponding to the
                # current variable
                df_data = list_df2plot[[i]]$data
                # Extracts the trend corresponding to the
                # current variable
                df_trend = list_df2plot[[i]]$trend
                p_threshold = list_df2plot[[i]]$p_threshold
                # Extracts the data corresponding to the code
                df_data_code = df_data[df_data$code == code,]
                # Extracts the trend corresponding to the code
                df_trend_code = df_trend[df_trend$code == code,]

                # Gets the associated time info
                Start = tab_Start[k, i, j]
                End = tab_End[k, i, j]
                Periods = tab_Periods[k, i, j]

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

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

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

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

    # Number of ticks for the colorbar
    nbTick = 10
    # For all variable
    for (i in 1:nbp) {
        # If there is a specified station code to highlight (mini map)
        # and there has already been one loop
        if (i > 1 & !is.null(codeLight)) {
            # Stop the for loop over the variable
            break
        }
        # Extracts the variable of the plot
        var = list_df2plot[[i]]$var
        # Createsa name for the map
        outname = paste('map_', var, sep='')
        # If there is the verbose option
        if (verbose) {
            # Prints the name of the map
            print(paste('Map for variable : ', var,
                        "   (", round(i/nbp*100, 1), " %)", 
                        sep=''))
        } 

        # If there is no specified station code to highlight (mini map)
        if (is.null(codeLight)) {
            # Sets the size of the countour
            sizefr = 0.45
            sizebs = 0.4
            sizerv = 0.3
        } else {
            sizefr = 0.35
            sizebs = 0.3
            sizerv = 0.2
        }

        # Stores the coordonate system 
        cf = coord_fixed()
        # Makes it the default one to remove useless warning
        cf$default = TRUE

        # Open a new plot with the personalise theme
        map = ggplot() + theme_void() +
            # theme(plot.background=element_rect(fill=NA,
                                               # color="#EC4899")) +
            # Fixed coordinate system (remove useless warning)
            cf +
            # Plot the background of France
            geom_polygon(data=df_france,
                         aes(x=long, y=lat, group=group),
                         color=NA, fill="grey97")
        # If the river shapefile exists
        if (!is.null(df_river)) {
            # Plot the river
            map = map +
                geom_path(data=df_river,
                          aes(x=long, y=lat, group=group),
                          color="grey85", size=sizerv)   
        }
        
        map = map +
            # Plot the hydrological basin
            geom_polygon(data=df_bassin,
                         aes(x=long, y=lat, group=group),
                         color="grey70", fill=NA, size=sizebs) +
            # Plot the hydrological sub-basin
            geom_polygon(data=df_subbassin,
                         aes(x=long, y=lat, group=group),
                         color="grey70", fill=NA, size=sizebs) +
            # Plot the countour of France
            geom_polygon(data=df_france,
                         aes(x=long, y=lat, group=group),
                         color="grey40", fill=NA, size=sizefr)

        # If the sea needs to be shown
        if (showSea) {
            # Leaves space around the France
            xlim = c(280000, 790000)
            ylim = c(6110000, 6600000)
        # Otherwise
        } else {
            # Leaves minimal space around France
            xlim = c(305000, 790000)
            ylim = c(6135000, 6600000)
        }

        # If there is no specified station code to highlight (mini map)
        if (is.null(codeLight)) {
            # Sets a legend scale start
            xmin = gpct(7, xlim, shift=TRUE)
            # Sets graduations
            xint = c(0, 10*1E3, 50*1E3, 100*1E3)
            # Sets the y postion
            ymin = gpct(5, ylim, shift=TRUE)
            # Sets the height of graduations
            ymax = ymin + gpct(1, ylim)
            # Size of the value
            size = 3
            # Size of the 'km' unit
            sizekm = 2.5
        # If there is a specified station code
        } else {
            # Same but with less graduation and smaller size
            xmin = gpct(2, xlim, shift=TRUE)
            xint = c(0, 100*1E3)
            ymin = gpct(1, ylim, shift=TRUE)
            ymax = ymin + gpct(3, ylim)
            size = 2
            sizekm = 1.5
        }
        
        map = map +
            # Adds the base line of the scale
            geom_line(aes(x=c(xmin, max(xint)+xmin),
                          y=c(ymin, ymin)),
                      color="grey40", size=0.2) +
            # Adds the 'km' unit
            annotate("text",
                     x=max(xint)+xmin+gpct(1, xlim), y=ymin,
                     vjust=0, hjust=0, label="km",
                     color="grey40", size=sizekm)
        # For all graduations
        for (x in xint) {
            map = map +
                # Draws the tick
                annotate("segment",
                         x=x+xmin, xend=x+xmin, y=ymin, yend=ymax,
                         color="grey40", size=0.2) +
                # Adds the value
                annotate("text",
                         x=x+xmin, y=ymax+gpct(0.5, ylim),
                         vjust=0, hjust=0.5, label=x/1E3,
                         color="grey40", size=size)
        }
        
        map = map +
            # Allows to crop shapefile without graphical problem
            coord_sf(xlim=xlim, ylim=ylim,
                     expand=FALSE)
        
        # If there is no margins specified
        if (is.null(margin)) {
            # Sets all margins to 0
            map = map + 
                theme(plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm"))
            # Otherwise
        } else {
            # Sets margins to the given ones
            map = map + 
                theme(plot.margin=margin)
        }

        # Blank vector to store data about station
        lon = c()
        lat = c()
        fill = c()
        shape = c()
        trend = c()
        p_threshold_Ok = c()
        # For all code
        for (k in 1:nCode) {
            # Gets the code
            code = Code[k]
            # Extracts the data corresponding to the current variable
            df_data = list_df2plot[[i]]$data
            # Extracts the trend corresponding to the
            # current variable
            df_trend = list_df2plot[[i]]$trend
            p_threshold = list_df2plot[[i]]$p_threshold
            # Extracts the data corresponding to the code
            df_data_code = df_data[df_data$code == code,]
            # Extracts the trend corresponding to the code
            df_trend_code = df_trend[df_trend$code == code,]

            # Gets the associated time info
            # Start = Start_code[Code_code == code][[1]][idPer]
            # End = End_code[Code_code == code][[1]][idPer]

            # Gets the associated time info
            Start = tab_Start[k, i, idPer]
            End = tab_End[k, i, idPer]
            Periods = tab_Periods[k, i, idPer]

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

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

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

            # Computes the color associated to the mean trend
            color_res = get_color(trendMean, 
                                  minTrendMean[idPer, i],
                                  maxTrendMean[idPer, i],
                                  palette_name='perso',
                                  reverse=TRUE,
                                  ncolor=256)
            # Computes the colorbar info 
            palette_res = get_palette(minTrendMean[idPer, i],
                                      maxTrendMean[idPer, i],
                                      palette_name='perso',
                                      reverse=TRUE,
                                      ncolor=256,
                                      nbTick=nbTick)

            # If it is significative
            if (df_trend_code_per$p <= p_threshold){
                # The computed color is stored
                filltmp = color_res
                # If the mean tend is positive
                if (trendMean >= 0) {
                    # Uses a triangle up for the shape of the marker
                    shapetmp = 24
                # If negative
                } else {
                    # Uses a triangle down for the shape of the marker
                    shapetmp = 25
                }
            # If it is not significative
            } else {
                # The fill color is grey
                filltmp = 'grey97'
                # The marker is a circle
                shapetmp = 21
            }

            # Extracts the localisation of the current station
            lontmp =
                df_meta[df_meta$code == code,]$L93X_m_BH            
            lattmp =
                df_meta[df_meta$code == code,]$L93Y_m_BH 

            # Stores all the parameters
            lon = c(lon, lontmp)
            lat = c(lat, lattmp)
            fill = c(fill, filltmp)
            shape = c(shape, shapetmp)
            trend = c(trend, trendMean)
            # If the trend analysis is significative a TRUE is stored
            p_threshold_Ok = c(p_threshold_Ok,
                               df_trend_code_per$p <= p_threshold)
        }
        # Creates a tibble to stores all the data to plot
        plot_map = tibble(lon=lon, lat=lat, fill=fill,
                          shape=shape, code=Code)

        # If there is no specified station code to highlight (mini map)
        if (is.null(codeLight)) {
            map = map +
                # Plots the trend point
                geom_point(data=plot_map,
                           aes(x=lon, y=lat),
                           shape=shape, size=5, stroke=1,
                           color='grey50', fill=fill)
        # If there is a specified station code
        } else {
            # Extract data of all stations not to highlight
            plot_map_codeNo = plot_map[plot_map$code != codeLight,]
            # Extract data of the station to highlight
            plot_map_code = plot_map[plot_map$code == codeLight,]

            # Plots only the localisation
            map = map +
                # For all stations not to highlight
                geom_point(data=plot_map_codeNo,
                           aes(x=lon, y=lat),
                           shape=21, size=0.5, stroke=0.5,
                           color='grey50', fill='grey50') +
                # For the station to highlight
                geom_point(data=plot_map_code,
                           aes(x=lon, y=lat),
                           shape=21, size=1.5, stroke=0.5,
                           color='grey40', fill='grey40')
        }
        
        # Extracts the position of the tick of the colorbar
        posTick = palette_res$posTick
        # Extracts the label of the tick of the colorbar
        labTick = palette_res$labTick
        # Extracts the color corresponding to the tick of the colorbar
        colTick = palette_res$colTick

        # Spreading of the colorbar
        valNorm = nbTick * 10
        # Normalisation of the position of ticks
        ytick = posTick / max(posTick) * valNorm
        # Formatting of label in pourcent
        labTick = as.character(round(labTick*100, 2))

        # X position of ticks all similar
        xtick = rep(0, times=nbTick)

        # Creates a tibble to store all parameters of colorbar
        plot_palette = tibble(xtick=xtick, ytick=ytick,
                              colTick=colTick, labTick=labTick)

        # New plot with void theme
        title = ggplot() + theme_void() +
            # Plots separation line
            geom_line(aes(x=c(-0.3, 3.3), y=c(0.05, 0.05)),
                      size=0.6, color="#00A3A8") +
            # Writes title
            geom_shadowtext(data=tibble(x=-0.3, y=0.2,
                                        label=var),
                            aes(x=x, y=y, label=label),
                            fontface="bold",
                            color="#00A3A8",
                            bg.colour="white",
                            hjust=0, vjust=0, size=10) +
            # X axis
            scale_x_continuous(limits=c(-1, 1 + 3),
                               expand=c(0, 0)) +
            # Y axis
            scale_y_continuous(limits=c(0, 10),
                               expand=c(0, 0)) +
            # Margin
            theme(plot.margin=margin(t=5, r=5, b=0, l=0, unit="mm"))
        
        # New plot with void theme
        pal = ggplot() + theme_void() +
            # Plots the point of the colorbar
            geom_point(data=plot_palette,
                       aes(x=xtick, y=ytick),
                       shape=21, size=5, stroke=1,
                       color='white', fill=colTick)

        pal = pal +
            # Name of the colorbar
            annotate('text',
                     x=-0.3, y= valNorm + 23,
                     label="Tendance",
                     hjust=0, vjust=0.5,
                     size=6, color='grey40') +
            # Unit legend of the colorbar
            annotate('text',
                     x=-0.2, y= valNorm + 13,
                     label=bquote(bold("% par an")),
                     hjust=0, vjust=0.5,
                     size=4, color='grey40')
        # For all the ticks
        for (j in 1:nbTick) {
            pal = pal +
                # Adds the value
                annotate('text', x=xtick[j]+0.3,
                         y=ytick[j],
                         label=bquote(bold(.(labTick[j]))),
                         hjust=0, vjust=0.7, 
                         size=3, color='grey40')
        }

        pal = pal +
            # Up triangle in the marker legend
            geom_point(aes(x=0, y=-20),
                       shape=24, size=4, stroke=1,
                       color='grey50', fill='grey97') +
            # Up triangle text legend
            annotate('text',
                     x=0.3, y=-20,
                     label=bquote(bold("Hausse significative à 10%")),
                     hjust=0, vjust=0.5,
                     size=3, color='grey40')

        pal = pal +
            # Circle in the marker legend
            geom_point(aes(x=0, y=-29),
                       shape=21, size=4, stroke=1,
                       color='grey50', fill='grey97') +
            # Circle text legend
            annotate('text',
                     x=0.3, y=-29,
                     label=bquote(bold("Non significatif à 10%")),
                     hjust=0, vjust=0.7,
                     size=3, color='grey40')

        pal = pal +
            # Down triangle in the marker legend
            geom_point(aes(x=0, y=-40),
                       shape=25, size=4, stroke=1,
                       color='grey50', fill='grey97') +
            # Down triangle text legend
            annotate('text',
                     x=0.3, y=-40,
                     label=bquote(bold("Baisse significative à 10%")),
                     hjust=0, vjust=0.5,
                     size=3, color='grey40')

        # Normalises all the trend values for each station
        # according to the colorbar
        yTrend = (trend - minTrendMean[idPer, i]) /
            (maxTrendMean[idPer, i] - minTrendMean[idPer, i]) * valNorm
        # Takes only the significative ones
        yTrend = yTrend[p_threshold_Ok]

        ## Random distribution ##
        # xTrend = rnorm(length(yTrend), mean=1.75, sd=0.1)

        ## Histogram distribution ##
        # Computes the histogram of the trend
        res_hist = hist(yTrend, breaks=ytick, plot=FALSE)
        # Extracts the number of counts per cells
        counts = res_hist$counts
        # Extracts limits of cells 
        breaks = res_hist$breaks
        # Extracts middle of cells 
        mids = res_hist$mids

        # Blank vectors to store position of points of
        # the distribution to plot
        xTrend = c()
        yTrend = c()
        # Start X position of the distribution
        start_hist = 1.25
        # X separation bewteen point
        hist_sep = 0.15
        # For all cells of the histogram
        for (ii in 1:length(mids)) {
            # If the count in the current cell is not zero
            if (counts[ii] != 0) {
                # Stores the X positions of points of the distribution
                # for the current cell
                xTrend = c(xTrend,
                           seq(start_hist,
                               start_hist+(counts[ii]-1)*hist_sep,
                               by=hist_sep))
            }
            # Stores the Y position which is the middle of the
            # current cell the number of times it has been counted
            yTrend = c(yTrend, rep(mids[ii], times=counts[ii]))
        }

        ## No touch distribution ##
        # start_hist = 1.25
        # min_xsep = 0.15
        # min_ysep = 4

        # xTrend = rep(start_hist, times=length(yTrend))
        
        # for (ii in 1:length(yTrend)) {

        #     yTrendtmp = yTrend
        #     yTrendtmp[ii] = 1E99
                
        #     isinf_ysep = abs(yTrendtmp - yTrend[ii]) < min_ysep

        #     if (any(isinf_ysep) & !all(xTrend[which(isinf_ysep)] > start_hist)) {
        #         xTrend[ii] = max(xTrend[which(isinf_ysep)]) + min_xsep
        #     }
        # }
        
        # Makes a tibble to plot the trend distribution
        plot_trend = tibble(xTrend=xTrend, yTrend=yTrend)
        
        pal = pal +
            # Plots the point of the trend distribution
            geom_point(data=plot_trend,
                       aes(x=xTrend, y=yTrend),
                       # shape=21, size=1,
                       # color="grey20", fill="grey20")
                       alpha=0.4)


        pal = pal +
            # Arrow to show a worsening of the situation
            geom_segment(aes(x=2.7, y=valNorm*0.75,
                             xend=2.7, yend=valNorm*0.25),
                         color='grey50', size=0.3,
                         arrow=arrow(length=unit(2, "mm"))) +
            # Text associated to the arrow
            annotate('text',
                     x=2.8, y=valNorm*0.5,
                     label= "Plus sévère",
                     angle=90,
                     hjust=0.5, vjust=1,
                     size=3, color='grey50')
        
        
        pal = pal +
            # X axis of the colorbar
            scale_x_continuous(limits=c(-1, 1 + 3),
                               expand=c(0, 0)) +
            # Y axis of the colorbar
            scale_y_continuous(limits=c(-60, valNorm + 35),
                               expand=c(0, 0)) +
            # Margin of the colorbar
            theme(plot.margin=margin(t=0, r=5, b=5, l=0, unit="mm"))

        # Stores the map, the title and the colorbarin a list
        Map = list(map, title, pal)
        # Arranges the graphical object
        plot = grid.arrange(grobs=Map, layout_matrix=
                                           matrix(c(1, 1, 1, 2,
                                                    1, 1, 1, 3),
                                                  nrow=2, byrow=TRUE))

        # If there is no specified station code to highlight (mini map)
        if (is.null(codeLight)) {
            # Saving matrix plot
            ggsave(plot=plot,
                   path=outdirTmp,
                   filename=paste(outname, '.pdf', sep=''),
                   width=29.7, height=21, units='cm', dpi=100)
        }
    }
    # Returns the map object
    return (map)
}