map.R

# \\\
# Copyright 2021-2022 Louis Hraut*1
#
# *1   INRAE, France
#      louis.heraut@inrae.fr
#
# This file is part of ash R toolbox.
#
# ash R toolbox is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# ash R toolbox is distributed in the hope that it will be useful, but 
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with ash R toolbox.  If not, see <https://www.gnu.org/licenses/>.
# ///
#
#
# plotting/map.R
#
#


## 1. MAP PANEL
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_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
        }
    }

    # Blank list to store time info by station code
    Start_code = vector(mode='list', length=nCode)
    End_code = vector(mode='list', length=nCode)
    Code_code = vector(mode='list', length=nCode)
    Periods_code = vector(mode='list', length=nCode)
    # For all the code
    for (j in 1:nCode) {
        # Gets the code
        code = Code[j]
        # 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))
        
        # Vector to store trend period
        Periods = c()
        # For all the period
        for (i in 1:nPeriod_max) {
            # Stocks period
            Periods = append(Periods, 
                             paste(substr(Start[i], 1, 4),
                                   substr(End[i], 1, 4),
                                   sep=' / '))
        }

        Start_code[[j]] = Start
        End_code[[j]] = End
        Code_code[[j]] = code
        Periods_code[[j]] = Periods
        
    }

    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,]

                Start = Start_code[Code_code == code][[1]][j]
                End = End_code[Code_code == code][[1]][j]
                Periods = Periods_code[Code_code == code][[1]][j]

                df_data_code_per =
                    df_data_code[df_data_code$Date >= Start 
                                 & df_data_code$Date <= End,]

                df_trend_code_per = 
                    df_trend_code[df_trend_code$period_start == Start 
                                  & df_trend_code$period_end == End,]

                Ntrend = nrow(df_trend_code_per)
                if (Ntrend > 1) {
                    df_trend_code_per = df_trend_code_per[1,]
                }
                
                dataMean = mean(df_data_code_per$Qm3s, na.rm=TRUE)
                trendMean = df_trend_code_per$trend / dataMean

                if (df_trend_code_per$p <= p_threshold){
                    TrendMean_code[j, i, k] = trendMean
                } else {
                    TrendMean_code[j, i, k] = NA
                }
            }
        }
    }

    minTrendMean = apply(TrendMean_code, c(1, 2), min, na.rm=TRUE)
    maxTrendMean = apply(TrendMean_code, c(1, 2), max, na.rm=TRUE)    

    ncolor = 256
    nbTick = 10
    # For all variable
    for (i in 1:nbp) {
        
        if (i > 1 & !is.null(codeLight)) {
            break
        }
        # Extract the variable of the plot
        type = list_df2plot[[i]]$type
        outname = paste('map_', type, sep='')
        if (verbose) {
            print(paste('Map for variable :', type))
        } 

        if (is.null(codeLight)) {
            sizefr = 0.45
            sizebs = 0.4
            sizerv = 0.3
        } else {
            sizefr = 0.35
            sizebs = 0.3
            sizerv = 0.2
        }
        
        map = ggplot() + theme_void() +
            
            # theme(plot.background=element_rect(fill=NA,
                                               # color="#EC4899")) +
            # Fixed coordinate system
            coord_fixed() +
            
            geom_polygon(data=df_france,
                         aes(x=long, y=lat, group=group),
                         color=NA, fill="grey97")

        if (!is.null(df_river)) {
            map = map +
                geom_path(data=df_river,
                          aes(x=long, y=lat, group=group),
                          color="grey85", size=sizerv)   
        }
        
        map = map +
            geom_polygon(data=df_bassin,
                         aes(x=long, y=lat, group=group),
                         color="grey70", fill=NA, size=sizebs) +
            
            geom_polygon(data=df_france,
                         aes(x=long, y=lat, group=group),
                         color="grey40", fill=NA, size=sizefr)
        
        if (showSea) {
            xlim = c(280000, 790000)
            ylim = c(6110000, 6600000)
        } else {
            xlim = c(305000, 790000)
            ylim = c(6135000, 6600000)
        }

        if (is.null(codeLight)) {
            xmin = gpct(7, xlim, shift=TRUE)
            xint = c(0, 10*1E3, 50*1E3, 100*1E3)
            ymin = gpct(5, ylim, shift=TRUE)
            ymax = ymin + gpct(1, ylim)
            size = 3
            sizekm = 2.5
        } else {
            xmin = gpct(5, 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 +
            
            geom_line(aes(x=c(xmin, max(xint)+xmin),
                          y=c(ymin, ymin)),
                      color="grey40", size=0.2) +
            annotate("text",
                     x=max(xint)+xmin+gpct(1, xlim), y=ymin,
                     vjust=0, hjust=0, label="km",
                     color="grey40", size=sizekm)

        for (x in xint) {
            map = map +
                annotate("segment",
                         x=x+xmin, xend=x+xmin, y=ymin, yend=ymax,
                         color="grey40", size=0.2) +
                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 +          
            coord_sf(xlim=xlim, ylim=ylim,
                     expand=FALSE)
            
        if (is.null(margin)) {
            map = map +
                theme(plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm"))
        } else {
            map = map +
                theme(plot.margin=margin)
        }
        
        lon = c()
        lat = c()
        fill = c()
        shape = c()
        trend = c()
        p_threshold_Ok = c()
        # For all code
        for (code in Code) {
            # 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,]
            
            Start = Start_code[Code_code == code][[1]][idPer]
            End = End_code[Code_code == code][[1]][idPer]
            
            df_data_code_per =
                df_data_code[df_data_code$Date >= Start 
                             & df_data_code$Date <= End,]
            
            df_trend_code_per = 
                df_trend_code[df_trend_code$period_start == Start 
                              & df_trend_code$period_end == End,]
            
            Ntrend = nrow(df_trend_code_per)
            if (Ntrend > 1) {
                df_trend_code_per = df_trend_code_per[1,]
            }
            
            dataMean = mean(df_data_code_per$Qm3s, na.rm=TRUE)
            trendMean = df_trend_code_per$trend / dataMean

            color_res = get_color(trendMean, 
                                  minTrendMean[idPer, i],
                                  maxTrendMean[idPer, i],
                                  palette_name='perso',
                                  reverse=TRUE,
                                  ncolor=ncolor)

            palette_res = get_palette(minTrendMean[idPer, i],
                                      maxTrendMean[idPer, i],
                                      palette_name='perso',
                                      reverse=TRUE,
                                      ncolor=ncolor,
                                      nbTick=nbTick)
            
            if (df_trend_code_per$p <= p_threshold){
                filltmp = color_res

                if (trendMean >= 0) {
                    shapetmp = 24
                } else {
                    shapetmp = 25
                }
                
            } else {                
                filltmp = 'grey97'
                shapetmp = 21
            }

            lontmp =
                df_meta[df_meta$code == code,]$L93X_m_BH            
            lattmp =
                df_meta[df_meta$code == code,]$L93Y_m_BH 

            lon = c(lon, lontmp)
            lat = c(lat, lattmp)
            fill = c(fill, filltmp)
            shape = c(shape, shapetmp)
            trend = c(trend, trendMean)
            p_threshold_Ok = c(p_threshold_Ok,
                               df_trend_code_per$p <= p_threshold)

        }
        
        plot_map = tibble(lon=lon, lat=lat, fill=fill,
                          shape=shape, code=Code)

        if (is.null(codeLight)) {
            map = map +
                geom_point(data=plot_map,
                           aes(x=lon, y=lat),
                           shape=shape, size=5, stroke=1,
                           color='grey50', fill=fill)
        } else {
            plot_map_codeNo = plot_map[plot_map$code != codeLight,]
            plot_map_code = plot_map[plot_map$code == codeLight,]
            
            map = map +
                
                geom_point(data=plot_map_codeNo,
                           aes(x=lon, y=lat),
                           shape=21, size=0.5, stroke=0.5,
                           color='grey70', fill='grey70') +
                
                geom_point(data=plot_map_code,
                           aes(x=lon, y=lat),
                           shape=21, size=1.5, stroke=0.5,
                           color='grey40', fill='grey40')
        }
                
        posTick = palette_res$posTick
        labTick = palette_res$labTick
        colTick = palette_res$colTick
        
        nbTickmod = length(posTick)

        valNorm = nbTickmod * 10
        ytick = posTick / max(posTick) * valNorm
        
        labTick = as.character(round(labTick*100, 2))
       
        xtick = rep(0, times=nbTickmod)

        plot_palette = tibble(xtick=xtick, ytick=ytick,
                              colTick=colTick, labTick=labTick)
        
        title = ggplot() + theme_void() +
            
            geom_line(aes(x=c(-0.3, 3.3), y=c(0.05, 0.05)),
                      size=0.6, color="#00A3A8") +
            
            geom_shadowtext(data=tibble(x=-0.3, y=0.2,
                                        label=type),
                            aes(x=x, y=y, label=label),
                            fontface="bold",
                            color="#00A3A8",
                            bg.colour="white",
                            hjust=0, vjust=0, size=10) +
            
            scale_x_continuous(limits=c(-1, 1 + 3),
                               expand=c(0, 0)) +
            
            scale_y_continuous(limits=c(0, 10),
                               expand=c(0, 0)) +
            
            theme(plot.margin=margin(t=5, r=5, b=0, l=0, unit="mm"))
        
        
        pal = ggplot() + theme_void() +
            
            geom_point(data=plot_palette,
                       aes(x=xtick, y=ytick),
                       shape=21, size=5, stroke=1,
                       color='white', fill=colTick)

        pal = pal +
            
            annotate('text',
                     x=-0.3, y= valNorm + 23,
                     label="Tendance",
                     hjust=0, vjust=0.5,
                     size=6, color='grey40') +
            
            annotate('text',
                     x=-0.2, y= valNorm + 13,
                     label=bquote(bold("% par an")),
                     hjust=0, vjust=0.5,
                     size=4, color='grey40')

        for (j in 1:nbTickmod) {
            pal = pal +
                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 +
            
            geom_point(aes(x=0, y=-20),
                       shape=24, size=4, stroke=1,
                       color='grey50', fill='grey97') +
            
            annotate('text',
                     x=0.3, y=-20,
                     label=bquote(bold("Hausse significative  10%")),
                     hjust=0, vjust=0.5,
                     size=3, color='grey40')

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

        pal = pal +
            
            geom_point(aes(x=0, y=-40),
                       shape=25, size=4, stroke=1,
                       color='grey50', fill='grey97') +
            
            annotate('text',
                     x=0.3, y=-40,
                     label=bquote(bold("Baisse significative  10%")),
                     hjust=0, vjust=0.5,
                     size=3, color='grey40')
        
        yTrend = (trend - minTrendMean[idPer, i]) /
            (maxTrendMean[idPer, i] - minTrendMean[idPer, i]) * valNorm

        yTrend = yTrend[p_threshold_Ok]

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

        ## Histogram distribution ##     
        res_hist = hist(yTrend, breaks=ytick, plot=FALSE)
        counts = res_hist$counts
        breaks = res_hist$breaks
        mids = res_hist$mids

        xTrend = c()
        yTrend = c()
        start_hist = 1.25
        hist_sep = 0.15
        
        for (ii in 1:length(mids)) {

            if (counts[ii] != 0) {
                xTrend = c(xTrend,
                           seq(start_hist,
                               start_hist+(counts[ii]-1)*hist_sep,
                               by=hist_sep))
            }
            
            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
        #     }
        # }
        
        
        plot_trend = tibble(xTrend=xTrend, yTrend=yTrend)
        
        pal = pal +
            geom_point(data=plot_trend,
                       aes(x=xTrend, y=yTrend),
                       # shape=21, size=1,
                       # color="grey20", fill="grey20")
                       alpha=0.4)


        pal = pal +
            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"))) +
            
            annotate('text',
                     x=2.8, y=valNorm*0.5,
                     label= "Plus svre",
                     angle=90,
                     hjust=0.5, vjust=1,
                     size=3, color='grey50')
        
        
        pal = pal +
            
            scale_x_continuous(limits=c(-1, 1 + 3),
                               expand=c(0, 0)) +
            
            scale_y_continuous(limits=c(-60, valNorm + 35),
                               expand=c(0, 0)) +
            
            theme(plot.margin=margin(t=0, r=5, b=5, l=0, unit="mm"))

        
        Map = list(map, title, pal)
        
        plot = grid.arrange(grobs=Map, layout_matrix=
                                           matrix(c(1, 1, 1, 2,
                                                    1, 1, 1, 3),
                                                  nrow=2, byrow=TRUE))

        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)
        }
    }
    return (map)
}