Skip to content
GitLab
Select Git revision
  • master default
Find file BlameHistoryPermalink
panel.R 85.07 KiB
1 
# Usefull library
2 
library(ggplot2)
3 
library(scales)
4 
library(qpdf)
5 
library(gridExtra)
6 
library(gridtext)
7 
library(dplyr)
8 
library(grid)
9 
library(ggh4x)
10 
library(RColorBrewer)
11 
library(rgdal)
12
13
14
15
16 
display_type = function (type, bold=FALSE) {
17
18 
    if (type == "QA") {
19 
        if (bold) {
20 
            disp = bquote(Q[A])
21 
        } else {
22 
            disp = bquote(bold(Q[A]))
23 
        }
24 
    } else if (type == "QMNA") {
25 
        if (bold) {
26 
            disp = bquote(Q[MNA])
27 
        } else {
28 
            disp = bquote(bold(Q[MNA]))
29 
        }
30 
    } else if (type == "VCN10") {
31 
        if (bold) {
32 
            disp = bquote(V[CN10])
33 
        } else {
34 
            disp = bquote(bold(V[CN10]))
35 
        }
36 
    }
37
38 
    return (disp)
39 
}
40
41
42 
# Personal theme
43 
theme_ash =
44 
    theme(
45 
        # White background
46 
        panel.background=element_rect(fill='white'),
47 
        # Font
48 
        text=element_text(family='sans'),
49 
        # Border of plot
50 
        panel.border = element_rect(color="grey85",
51 
                                    fill=NA,
52 
                                    size=0.7),
53 
        # Grid
54 
        panel.grid.major.x=element_blank(),
55 
        panel.grid.major.y=element_blank(),
56 
        # Ticks marker
57 
        axis.ticks.x=element_line(color='grey75', size=0.3),
58 
        axis.ticks.y=element_line(color='grey75', size=0.3),
59 
        # Ticks label
60 
        axis.text.x=element_text(color='grey40'),
61 
        axis.text.y=element_text(color='grey40'),
62 
        # Ticks length
63 
        axis.ticks.length=unit(1.5, 'mm'),
64 
        # Ticks minor
65 
        ggh4x.axis.ticks.length.minor=rel(0.5),
66 
        # Title
67 
        plot.title=element_text(size=9, vjust=-2,
68 
                                hjust=-1E-3, color='grey20'),
69 
        # Axis title
70 
        axis.title.x=element_blank(),
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
axis.title.y=element_blank(), # Axis line axis.line.x=element_blank(), axis.line.y=element_blank(), ) time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missRect=FALSE, unit2day=365.25, trend_period=NULL, mean_period=NULL, axis_xlim=NULL, last=FALSE, first=FALSE, color=NULL) { # If 'type' is square root apply it to data if (type == 'sqrt(Q)') { df_data_code$Qm3s = sqrt(df_data_code$Qm3s) } # Compute max of flow maxQ = max(df_data_code$Qm3s, na.rm=TRUE) # Get the magnitude of the max of flow power = get_power(maxQ) # Normalize the max flow by it's magnitude maxQtmp = maxQ/10^power # Compute the spacing between y ticks if (maxQtmp >= 5) { dbrk = 1.0 } else if (maxQtmp < 5 & maxQtmp >= 3) { dbrk = 0.5 } else if (maxQtmp < 3 & maxQtmp >= 2) { dbrk = 0.4 } else if (maxQtmp < 2 & maxQtmp >= 1) { dbrk = 0.2 } else if (maxQtmp < 1) { dbrk = 0.1 } # Get the spacing in the right magnitude dbrk = dbrk * 10^power # Fix the accuracy for label accuracy = NULL # Time span in the unit of time dDate = as.numeric(df_data_code$Date[length(df_data_code$Date)] - df_data_code$Date[1]) / unit2day # Compute the spacing between x ticks if (dDate >= 100) { datebreak = 25 dateminbreak = 5 } else if (dDate < 100 & dDate >= 50) { datebreak = 10 dateminbreak = 1 } else if (dDate < 50) { datebreak = 5 dateminbreak = 1 } # Open new plot p = ggplot() + theme_ash # If it is the lats plot of the pages or not if (last) { if (first) { p = p + theme(plot.margin=margin(5, 5, 5, 5, unit="mm")) } else { p = p + theme(plot.margin=margin(0, 5, 5, 5, unit="mm")) } # If it is the first plot of the pages or not } else { if (first) {
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
p = p + theme(plot.margin=margin(5, 5, 0, 5, unit="mm")) } else { p = p + theme(plot.margin=margin(0, 5, 0, 5, unit="mm")) } } ## Sub period background ## if (!is.null(trend_period)) { # trend_period = as.list(trend_period) # Imin = 10^99 # for (per in trend_period) { # I = interval(per[1], per[2]) # if (I < Imin) { # Imin = I # trend_period_min = as.Date(per) # } # } # p = p + # geom_rect(aes(xmin=min(df_data_code$Date), # ymin=0, # xmax=trend_period_min[1], # ymax= maxQ*1.1), # linetype=0, fill='grey97') + # geom_rect(aes(xmin=trend_period_min[2], # ymin=0, # xmax=max(df_data_code$Date), # ymax= maxQ*1.1), # linetype=0, fill='grey97') # Convert trend period to list if it is not trend_period = as.list(trend_period) # Fix a disproportionate minimum for period Imin = 10^99 # For all the sub period of analysis in 'trend_period' for (per in trend_period) { # Compute time interval of period I = interval(per[1], per[2]) # If it is the smallest interval if (I < Imin) { # Store it Imin = I # Fix min period of analysis trend_period_min = as.Date(per) } } # Search for the index of the closest existing date # to the start of the min period of analysis idMinPer = which.min(abs(df_data_code$Date - trend_period_min[1])) # Same for the end of the min period of analysis idMaxPer = which.min(abs(df_data_code$Date - trend_period_min[2])) # Get the start and end date associated minPer = df_data_code$Date[idMinPer] maxPer = df_data_code$Date[idMaxPer] # If it is not a flow or sqrt of flow time serie if (type != 'sqrt(Q)' & type != 'Q') { # If there is an 'axis_lim' if (!is.null(axis_xlim)) { # If the temporary start of period is smaller # than the fix start of x axis limit if (minPer < axis_xlim[1]) { # Set the start of the period to the start of # the x axis limit minPer = axis_xlim[1] } }
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
} # If it is not a flow or sqrt of flow time serie if (type != 'sqrt(Q)' & type != 'Q') { # If there is an 'axis_lim' if (!is.null(axis_xlim)) { # If the temporary end of period plus one year # is smaller than the fix end of x axis limit if (maxPer + years(1) < axis_xlim[2]) { # Add one year the the temporary end of period maxPer = maxPer + years(1) } else { # Set the start of the period to the start of # the x axis limit maxPer = axis_xlim[2] } # Add one year the the temporary end of period # if there is no 'axis_lim' } else { maxPer = maxPer + years(1) } } # Draw rectangle to delimiting the sub period p = p + geom_rect(aes(xmin=minPer, ymin=0, xmax=maxPer, ymax= maxQ*1.1), linetype=0, fill='grey97') } ## Mean step ## # If there is a 'mean_period' if (!is.null(mean_period)) { # Convert 'mean_period' to list mean_period = as.list(mean_period) # Number of mean period nPeriod_mean = length(mean_period) # Blank tibble to store variable in order to plot # rectangle for mean period plot_mean = tibble() # Blank tibble to store variable in order to plot # upper limit of rectangle for mean period plot_line = tibble() # For all mean period for (j in 1:nPeriod_mean) { # 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,] # Min for the sub period xmin = min(df_data_code_per$Date) # If the min over the sub period is greater # than the min of the entier period and # it is not the first sub period if (xmin > min(df_data_code$Date) & j != 1) { # Substract 6 months to be in the middle of # the previous year xmin = xmin - months(6) } # If it is not a flow or sqrt of flow time serie and # it is the first period if (type != 'sqrt(Q)' & type != 'Q' & j == 1) { # If there is an x axis limit
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
if (!is.null(axis_xlim)) { # If the min of the period is before the x axis min if (xmin < axis_xlim[1]) { # The min for the sub period is the x axis xmin = axis_xlim[1] } } } # Max for the sub period xmax = max(df_data_code_per$Date) # If the max over the sub period is smaller # than the max of the entier period and # it is not the last sub period if (xmax < max(df_data_code$Date) & j != nPeriod_mean) { # Add 6 months to be in the middle of # the following year xmax = xmax + months(6) } # If it is not a flow or sqrt of flow time serie and # it is the last period if (type != 'sqrt(Q)' & type != 'Q' & j == nPeriod_mean) { # If there is an x axis limit if (!is.null(axis_xlim)) { # If the max of the period plus 1 year # is smaller thant the max of the x axis limit if (xmax + years(1) < axis_xlim[2]) { # Add one year to the max to include # the entire last year graphically xmax = xmax + years(1) } else { # The max of this sub period is the max # of the x axis limit xmax = axis_xlim[2] } # If there is no axis limit } else { # Add one year to the max to include # the entire last year graphically xmax = xmax + years(1) } } # Mean of the flow over the sub period ymax = mean(df_data_code_per$Qm3s, na.rm=TRUE) # Create temporary tibble with variable # to create rectangle for mean step plot_meantmp = tibble(xmin=xmin, xmax=xmax, ymin=0, ymax=ymax, period=j) # Bind it to the main tibble to store it with other period plot_mean = bind_rows(plot_mean, plot_meantmp) # Create vector for the upper limit of the rectangle abs = c(xmin, xmax) ord = c(ymax, ymax) # Create temporary tibble with variable # to create upper limit for rectangle plot_linetmp = tibble(abs=abs, ord=ord, period=j) # Bind it to the main tibble to store it with other period plot_line = bind_rows(plot_line, plot_linetmp) } # Plot rectangles p = p + geom_rect(data=plot_mean, aes(xmin=xmin, ymin=ymin, xmax=xmax, ymax=ymax), linetype=0, fill='grey93') # Plot upper line for rectangle p = p + geom_line(data=plot_line,
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
aes(x=abs, y=ord, group=period), color='grey85', size=0.15) # for all the sub periods except the last one for (i in 1:(nPeriod_mean-1)) { # The y limit of rectangle is the max of # the two neighboring mean step rectangle yLim = max(c(plot_mean$ymax[i], plot_mean$ymax[i+1])) # The x limit is the x max of the ith rectangle xLim = plot_mean$xmax[i] # Make a tibble to store data plot_lim = tibble(x=c(xLim, xLim), y=c(0, yLim)) # Plot the limit of rectangles p = p + geom_line(data=plot_lim, aes(x=x, y=y), linetype='dashed', size=0.15, color='grey85') } } ### Grid ### # If there is no axis limit if (is.null(axis_xlim)) { # The min and the max is set by # the min and the max of the date data xmin = min(df_data_code$Date) xmax = max(df_data_code$Date) } else { # Min and max is set with the limit axis parameter xmin = axis_xlim[1] xmax = axis_xlim[2] } # Create a vector for all the y grid position ygrid = seq(0, maxQ*10, dbrk) # Blank vector to store position ord = c() abs = c() # For all the grid element for (i in 1:length(ygrid)) { # Store grid position ord = c(ord, rep(ygrid[i], times=2)) abs = c(abs, xmin, xmax) } # Create a tibble to store all the position plot_grid = tibble(abs=as.Date(abs), ord=ord) # Plot the y grid p = p + geom_line(data=plot_grid, aes(x=abs, y=ord, group=ord), color='grey85', size=0.15) ### Data ### # If it is a square root flow or flow if (type == 'sqrt(Q)' | type == 'Q') { # Plot the data as line p = p + geom_line(aes(x=df_data_code$Date, y=df_data_code$Qm3s), color='grey20', size=0.3, lineend="round") } else { # Plot the data as point p = p + geom_point(aes(x=df_data_code$Date, y=df_data_code$Qm3s), shape=21, color='grey50', fill='grey97', size=1) } ### Missing data ###
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490
# If the option is TRUE if (missRect) { # Remove NA data NAdate = df_data_code$Date[is.na(df_data_code$Qm3s)] # Get the difference between each point of date data without NA dNAdate = diff(NAdate) # If difference of day is not 1 then # it is TRUE for the beginning of each missing data period NAdate_Down = NAdate[append(Inf, dNAdate) != 1] # If difference of day is not 1 then # it is TRUE for the ending of each missing data period NAdate_Up = NAdate[append(dNAdate, Inf) != 1] # Plot the missing data period p = p + geom_rect(aes(xmin=NAdate_Down, ymin=0, xmax=NAdate_Up, ymax=maxQ*1.1), linetype=0, fill='Wheat', alpha=0.4) } ### Trend ### # If there is trends if (!is.null(df_trend_code)) { # Extract starting period of trends Start = df_trend_code$period_start # Get the name of the different period UStart = levels(factor(Start)) # Same for ending End = df_trend_code$period_end UEnd = levels(factor(End)) # Compute the max of different start and end # so the number of different period nPeriod_trend = max(length(UStart), length(UEnd)) # Blank tibble to store trend data and legend data plot_trend = tibble() leg_trend = tibble() # For all the different period for (i in 1:nPeriod_trend) { # Get the trend associated to the first period df_trend_code_per = df_trend_code[df_trend_code$period_start == Start[i] & df_trend_code$period_end == End[i],] # Number of trend selected Ntrend = nrow(df_trend_code_per) # If the number of trend is greater than a unique one if (Ntrend > 1) { # Extract only the first hence it is the same period df_trend_code_per = df_trend_code_per[1,] } # Search for the index of the closest existing date # to the start of the trend period of analysis iStart = which.min(abs(df_data_code$Date - Start[i])) # Same for the end iEnd = which.min(abs(df_data_code$Date - End[i])) # Get the start and end date associated xmin = df_data_code$Date[iStart] xmax = df_data_code$Date[iEnd] # If there is a x axis limit if (!is.null(axis_xlim)) { # If the min of the current period
491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560
# is smaller than the min of the x axis limit if (xmin < axis_xlim[1]) { # The min of the period is the min # of the x axis limit xmin = axis_xlim[1] } # Same for end if (xmax > axis_xlim[2]) { xmax = axis_xlim[2] } } # Create vector to store x data abs = c(xmin, xmax) # Convert the number of day to the unit of the period abs_num = as.numeric(abs) / unit2day # Compute the y of the trend ord = abs_num * df_trend_code_per$trend + df_trend_code_per$intercept # Create temporary tibble with variable to plot trend # for each period plot_trendtmp = tibble(abs=abs, ord=ord, period=i) # Bind it to the main tibble to store it with other period plot_trend = bind_rows(plot_trend, plot_trendtmp) # If there is a x axis limit if (!is.null(axis_xlim)) { # The x axis limit is selected codeDate = axis_xlim } else { # The entire date data is selected codeDate = df_data_code$Date } # The flow data is extract codeQ = df_data_code$Qm3s # Position of the x beginning and end of the legend symbol x = gpct(2, codeDate, shift=TRUE) xend = x + gpct(3, codeDate) # Position of the y beginning and end of the legend symbol dy = gpct(7, codeQ, ref=0) y = gpct(100, codeQ, ref=0) - (i-1)*dy yend = y # Position of x for the beginning of the associated text xt = xend + gpct(1, codeDate) # Position of the background rectangle of the legend xminR = x - gpct(1, codeDate) yminR = y - gpct(4, codeQ, ref=0) xmaxR = x + gpct(24, codeDate) ymaxR = y + gpct(5, codeQ, ref=0) # Get the tendance analyse trend = df_trend_code_per$trend # Compute the magnitude of the trend power = get_power(trend) # Convert it to character powerC = as.character(power) # Get the power of ten of magnitude brk = 10^power # Convert trend to character for sientific expression trendC = as.character(round(trend / brk, 2)) # Create temporary tibble with variable to plot legend leg_trendtmp = tibble(x=x, xend=xend, y=y, yend=yend, xt=xt,
561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630
trendC=trendC, powerC=powerC, xminR=xminR, yminR=yminR, xmaxR=xmaxR, ymaxR=ymaxR, period=i) # Bind it to the main tibble to store it with other period leg_trend = bind_rows(leg_trend, leg_trendtmp) } # For all periods for (i in 1:nPeriod_trend) { # Extract the trend of the current sub period leg_trend_per = leg_trend[leg_trend$period == i,] # Plot the background for legend p = p + geom_rect(data=leg_trend_per, aes(xmin=xminR, ymin=yminR, xmax=xmaxR, ymax=ymaxR), linetype=0, fill='white', alpha=0.5) } # For all periods for (i in 1:nPeriod_trend) { # Extract the trend of the current sub period leg_trend_per = leg_trend[leg_trend$period == i,] # Get the character variable for naming the trend trendC = leg_trend_per$trendC powerC = leg_trend_per$powerC # Create the name of the trend label = bquote(bold(.(trendC)~'x'~'10'^{.(powerC)})~'['*m^{3}*'.'*s^{-1}*'.'*an^{-1}*']') # Plot the trend symbole and value of the legend p = p + annotate("segment", x=leg_trend_per$x, xend=leg_trend_per$xend, y=leg_trend_per$y, yend=leg_trend_per$yend, color=color[i], linetype='solid', lwd=1) + annotate("text", label=label, size=3, x=leg_trend_per$xt, y=leg_trend_per$y, hjust=0, vjust=0.4, color=color[i]) } # For all periods for (i in 1:nPeriod_trend) { # Extract the trend of the current sub period plot_trend_per = plot_trend[plot_trend$period == i,] # Plot the line of white background of each trend p = p + geom_line(data=plot_trend_per, aes(x=abs, y=ord), color='white', linetype='solid', size=1, lineend="round") } # For all periods for (i in 1:nPeriod_trend) { # Extract the trend of the current sub period
631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700
plot_trend_per = plot_trend[plot_trend$period == i,] # Plot the line of trend p = p + geom_line(data=plot_trend_per, aes(x=abs, y=ord), color=color[i], linetype='solid', size=0.5, lineend="round") } } # Title p = p + ggtitle(bquote(bold(.(type))~~'['*m^{3}*'.'*s^{-1}*']')) # If the is no x axis limit if (is.null(axis_xlim)) { # Parameters of the x axis contain the limit of the date data p = p + scale_x_date(date_breaks=paste( as.character(datebreak), 'year', sep=' '), date_minor_breaks=paste( as.character(dateminbreak), 'year', sep=' '), guide='axis_minor', date_labels="%Y", limits=c(min(df_data_code$Date), max(df_data_code$Date)), expand=c(0, 0)) } else { # Parameters of the x axis contain the x axis limit p = p + scale_x_date(date_breaks=paste( as.character(datebreak), 'year', sep=' '), date_minor_breaks=paste( as.character(dateminbreak), 'year', sep=' '), guide='axis_minor', date_labels="%Y", limits=axis_xlim, expand=c(0, 0)) } # Parameters of the y axis p = p + scale_y_continuous(breaks=seq(0, maxQ*10, dbrk), limits=c(0, maxQ*1.1), expand=c(0, 0), labels=label_number(accuracy=accuracy)) return(p) } matrice_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice=NULL, outdirTmp='', outnameTmp='matrix', title=NULL, A3=FALSE) { nbp = length(list_df2plot) # Get all different stations code Code = levels(factor(df_meta$code)) nCode = length(Code) # nbp = length(list_df2plot) # # Get all different stations code # Code = levels(factor(df_meta$code))
701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770
# nCode = length(Code) 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) nPeriod_max = 0 for (code in Code) { df_trend_code = df_trend[df_trend$code == code,] Start = df_trend_code$period_start UStart = levels(factor(Start)) End = df_trend_code$period_end UEnd = levels(factor(End)) nPeriod = max(length(UStart), length(UEnd)) if (nPeriod > nPeriod_max) { nPeriod_max = nPeriod } } # print(nPeriod_trend) # print(nPeriod_max) 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 (j in 1:nCode) { code = Code[j] df_trend_code = df_trend[df_trend$code == code,] Start = df_trend_code$period_start UStart = levels(factor(Start)) End = df_trend_code$period_end UEnd = levels(factor(End)) nPeriod = max(length(UStart), length(UEnd)) Periods = c() for (i in 1:nPeriod_trend) { Periods = append(Periods, paste(Start[i], End[i], sep=' / ')) } Start_code[[j]] = Start End_code[[j]] = End Code_code[[j]] = code Periods_code[[j]] = Periods } TrendMean_code = array(rep(1, nPeriod_trend*nbp*nCode), dim=c(nPeriod_trend, nbp, nCode)) for (j in 1:nPeriod_trend) {
771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840
for (k in 1:nCode) { code = Code[k] for (i in 1:nbp) { df_data = list_df2plot[[i]]$data df_trend = list_df2plot[[i]]$trend p_threshold = list_df2plot[[i]]$p_threshold df_data_code = df_data[df_data$code == 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 TrendMean_code[j, i, k] = trendMean } } } minTrendMean = apply(TrendMean_code, c(1, 2), min, na.rm=TRUE) maxTrendMean = apply(TrendMean_code, c(1, 2), max, na.rm=TRUE) Periods_trend = c() NPeriod_trend = c() Type_trend = list() Code_trend = c() Pthresold_trend = c() TrendMean_trend = c() DataMean_trend = c() Fill_trend = c() Color_trend = c() for (j in 1:nPeriod_trend) { for (code in Code) { for (i in 1:nbp) { df_data = list_df2plot[[i]]$data df_trend = list_df2plot[[i]]$trend p_threshold = list_df2plot[[i]]$p_threshold type = list_df2plot[[i]]$type df_data_code = df_data[df_data$code == 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]
841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910
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){ color_res = get_color(trendMean, minTrendMean[j, i], maxTrendMean[j, i], palette_name='perso', reverse=TRUE) fill = color_res color = 'white' Pthresold = p_thresold } else { fill = 'white' color = 'grey85' Pthresold = NA } Periods_trend = append(Periods_trend, Periods) NPeriod_trend = append(NPeriod_trend, j) Type_trend = append(Type_trend, type) 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)) { Periods_mean = c() NPeriod_mean = c() Type_mean = list() 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) BreakMean_code = array(rep(1, nPeriod_mean*nbp*nCode), dim=c(nPeriod_mean, nbp, nCode)) dataMeantmp = array(rep(NA, nbp*nCode), dim=c(nbp, nCode)) # For all mean period for (j in 1:nPeriod_mean) {
911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980
for (k in 1:nCode) { code = Code[k] for (i in 1:nbp) { df_data = list_df2plot[[i]]$data type = list_df2plot[[i]]$type 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,] # Min for the sub period Datemin = min(df_data_code_per$Date) # Max for the sub period Datemax = max(df_data_code_per$Date) Periods = paste(Datemin, Datemax, sep=' / ') # Mean of the flow over the sub period dataMean = mean(df_data_code_per$Qm3s, na.rm=TRUE) if (j > 1) { Break = dataMean - dataMeantmp[i, k] } else { Break = NA } BreakMean = Break / dataMeantmp[i, k] BreakMean_code[j, i, k] = BreakMean dataMeantmp[i, k] = dataMean Periods_mean = append(Periods_mean, Periods) NPeriod_mean = append(NPeriod_mean, j) Type_mean = append(Type_mean, type) Code_mean = append(Code_mean, code) DataMean_mean = append(DataMean_mean, dataMean) BreakMean_mean = append(BreakMean_mean, BreakMean) } } } minBreakMean = apply(BreakMean_code, c(1, 2), min, na.rm=TRUE) maxBreakMean = apply(BreakMean_code, c(1, 2), max, na.rm=TRUE) Fill_mean = c() Color_mean = c() ii = 1 for (j in 1:nPeriod_mean) { for (k in 1:nCode) { code = Code[k]
981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050
for (i in 1:nbp) { BreakMean = BreakMean_mean[ii] color_res = get_color(BreakMean, minBreakMean[j, i], maxBreakMean[j, i], palette_name='perso', reverse=TRUE) fill = color_res color = 'white' Fill_mean = append(Fill_mean, fill) Color_mean = append(Color_mean, color) ii = ii + 1 } } } } if (is.null(slice)) { slice = nCode } firstLetter = levels(factor(substr(Code, 1, 1))) for (fL in firstLetter) { print(paste('matrix for region :', fL)) subCodefL = Code[substr(Code, 1, 1) == fL] nsubCodefL = length(subCodefL) nMat = as.integer(nsubCodefL/slice) + 1 print(nMat) for (imat in 1:nMat) { subCode = subCodefL[(slice*(imat-1)+1):(slice*imat)] subCode = subCode[!is.na(subCode)] # print('aaa') # print(subCode) nsubCode = length(subCode) CodefL_trend = substr(Code_trend, 1, 1) == fL subPeriods_trend = Periods_trend[CodefL_trend] subNPeriod_trend = NPeriod_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] # print(subNPeriod_trend) CodefL_mean = substr(Code_mean, 1, 1) == fL subPeriods_mean = Periods_mean[CodefL_mean] subNPeriod_mean = NPeriod_mean[CodefL_mean] subType_mean = Type_mean[CodefL_mean]
1051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120
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] title = df_meta[df_meta$code == subCode[1],]$region_hydro # print('bbb') ### Plot ### height = length(subCode) width = nbp * 2 * nPeriod_trend + nPeriod_trend + nPeriod_mean * nbp + nPeriod_mean + nbp options(repr.plot.width=width, repr.plot.height=height) mat = ggplot() + theme_ash + 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") ) xt = 1 - 6 yt = height + 2 Title = bquote(bold(.(title))) mat = mat + annotate("text", x=xt, y=yt, label=Title, hjust=0, vjust=1, size=6, color="#00A3A8") # print('ccc') ### Trend ### for (j in 1:nPeriod_trend) { 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] Xtmp = as.integer(factor(as.character(Type_trend_per))) Xc = j + (j - 1)*nbp*2 Xm = Xtmp + (j - 1)*nbp*2 + j X = Xtmp + (j - 1)*nbp*2 + nbp + j Y = as.integer(factor(Code_trend_per)) x = Xc - 0.4 xend = X[length(X)] + 0.25 y = height + 1
1121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190
yend = height + 1 mat = mat + annotate("segment", x=x, xend=xend, y=y, yend=yend, color="grey40", size=0.35) yt = y + 0.15 Start = trend_period[[j]][1] End = trend_period[[j]][2] periodName = bquote(bold('Priode')~bold(.(as.character(j)))) # bquote(bold(.(Start))~'/'~bold(.(End))) mat = mat + annotate("text", x=x, y=yt, label=periodName, hjust=0, vjust=0.5, size=3, color='grey40') for (i in 1:length(X)) { mat = mat + gg_circle(r=0.45, xc=X[i], yc=Y[i], fill=Fill_trend_per[i], color=Color_trend_per[i]) + gg_circle(r=0.45, xc=Xm[i], yc=Y[i], fill='white', color='grey40') + gg_circle(r=0.45, xc=Xc, yc=Y[i], fill='white', color='grey40') } for (i in 1:length(TrendMean_trend_per)) { trendMean = TrendMean_trend_per[i] trendC = signif(trendMean*100, 2) if (!is.na(Pthresold_trend_per[i])) { Tcolor = 'white' } else { Tcolor = 'grey85' } dataMean = signif(DataMean_trend_per[i], 2) mat = mat + annotate('text', x=X[i], y=Y[i], label=trendC, hjust=0.5, vjust=0.5, size=3, color=Tcolor) + annotate('text', x=Xm[i], y=Y[i], label=dataMean, hjust=0.5, vjust=0.5, size=3, color='grey40') } mat = mat + annotate('text', x=Xc, y=max(Y) + 0.85, label=bquote(bold('Dbut')), hjust=0.5, vjust=0.5, size=3, color='grey20') + annotate('text', x=Xc, y=max(Y) + 0.6, label=bquote(bold('Fin')), hjust=0.5, vjust=0.5, size=3, color='grey20') for (i in 1:nbp) {
1191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260
type = list_df2plot[[i]]$type mat = mat + annotate('text', x=X[i], y=max(Y) + 0.7, label=bquote(.(type)), hjust=0.5, vjust=0.5, size=3.5, color='grey20') + annotate('text', x=Xm[i], y=max(Y) + 0.7, label=bquote(''*.(type)), hjust=0.5, vjust=0.5, size=3.5, color='grey20') } for (k in 1:nsubCode) { code = subCode[k] label = Periods_trend[NPeriod_trend == j & subCode_trend == code][1] periodStart = substr(label, 1, 4) periodEnd = substr(label, 14, 17) mat = mat + annotate('text', x=Xc, y=k + 0.13, label=bquote(bold(.(periodStart))), hjust=0.5, vjust=0.5, size=3, color='grey40') + annotate('text', x=Xc, y=k - 0.13, label=bquote(bold(.(periodEnd))), hjust=0.5, vjust=0.5, size=3, color='grey40') } } # print('ddd') ### Mean ### for (j in 1:nPeriod_mean) { 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] # print('ddd1') Xtmp_mean = as.integer(factor(as.character(Type_mean_per))) Xc_mean = j + (j - 1)*nbp + X[length(X)] Xm_mean = Xtmp_mean + (j - 1)*nbp + j + X[length(X)] Xr_mean = Xtmp_mean + (j - 1)*nbp*2 + j + X[length(X)] Y_mean = as.integer(factor(Code_mean_per)) x = Xc_mean - 0.4 xend = Xm_mean[length(Xm_mean)] + 0.25 y = height + 1 yend = height + 1 mat = mat + annotate("segment",
1261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330
x=x, xend=xend, y=y, yend=yend, color="grey40", size=0.35) yt = y + 0.15 Start = mean_period[[j]][1] End = mean_period[[j]][2] periodName = bquote(bold('Priode')~bold(.(as.character(j+nPeriod_trend)))) mat = mat + annotate("text", x=x, y=yt, label=periodName, hjust=0, vjust=0.5, size=3, color='grey40') # print('ddd2') if (j > 1) { x = Xr_mean[1] - 0.4 xend = Xr_mean[length(Xr_mean)] + 0.25 mat = mat + annotate("segment", x=x, xend=xend, y=y, yend=yend, color="grey40", size=0.35) breakName = bquote(bold('Rupture')~bold(.(as.character(j-1+nPeriod_trend)))*bold('-')*bold(.(as.character(j+nPeriod_trend)))) mat = mat + annotate("text", x=x, y=yt, label=breakName, hjust=0, vjust=0.5, size=3, color='grey40') } # print('ddd3') for (i in 1:length(Xm_mean)) { mat = mat + gg_circle(r=0.45, xc=Xm_mean[i], yc=Y[i], fill='white', color='grey40') + gg_circle(r=0.45, xc=Xc_mean, yc=Y[i], fill='white', color='grey40') if (j > 1) { mat = mat + gg_circle(r=0.45, xc=Xr_mean[i], yc=Y[i], fill=Fill_mean_per[i], color=Color_mean_per[i]) } } # print('ddd4') for (i in 1:length(DataMean_mean_per)) { dataMean = signif(DataMean_mean_per[i], 2) # print(i) mat = mat + annotate('text', x=Xm_mean[i], y=Y[i], label=dataMean, hjust=0.5, vjust=0.5, size=3, color='grey40') if (j > 1) {
1331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400
BreakMean = BreakMean_mean_per[i] BreakC = signif(BreakMean*100, 2) # print(BreakMean) mat = mat + annotate('text', x=Xr_mean[i], y=Y[i], label=BreakC, hjust=0.5, vjust=0.5, size=3, color='white') } } # print('ddd5') mat = mat + annotate('text', x=Xc_mean, y=max(Y) + 0.85, label=bquote(bold('Dbut')), hjust=0.5, vjust=0.5, size=3, color='grey20') + annotate('text', x=Xc_mean, y=max(Y) + 0.6, label=bquote(bold('Fin')), hjust=0.5, vjust=0.5, size=3, color='grey20') # print('ddd6') for (i in 1:nbp) { type = list_df2plot[[i]]$type mat = mat + annotate('text', x=Xm_mean[i], y=max(Y) + 0.7, label=bquote(''*.(type)), hjust=0.5, vjust=0.5, size=3.5, color='grey20') if (j > 1) { mat = mat + annotate('text', x=Xr_mean[i], y=max(Y) + 0.7, label=bquote('d'*.(type)), hjust=0.5, vjust=0.5, size=3.5, color='grey20') } } # print('ddd7') for (k in 1:nsubCode) { code = subCode[k] label = Periods_mean[NPeriod_mean == j & subCode_mean == code][1] periodStart = substr(label, 1, 4) periodEnd = substr(label, 14, 17) mat = mat + annotate('text', x=Xc_mean, y=k + 0.13, label=bquote(bold(.(periodStart))), hjust=0.5, vjust=0.5, size=3, color='grey40') + annotate('text', x=Xc_mean, y=k - 0.13, label=bquote(bold(.(periodEnd))), hjust=0.5, vjust=0.5, size=3, color='grey40') } } # print('eee')
1401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470
### Code ### for (k in 1:nsubCode) { code = subCode[k] name = df_meta[df_meta$code == code,]$nom ncharMax = 40 if (nchar(name) > ncharMax) { name = paste(substr(name, 1, ncharMax), '...', sep='') } mat = mat + annotate('text', x=0.3, y=k + 0.14, label=bquote(bold(.(code))), hjust=1, vjust=0.5, size=3.5, color="#00A3A8") + annotate('text', x=0.3, y=k - 0.14, label=name, hjust=1, vjust=0.5, size=3.5, color="#00A3A8") } # print('fff') ### Environment ### mat = mat + coord_fixed() + scale_x_continuous(limits=c(1 - rel(6), width + rel(0.5)), expand=c(0, 0)) + scale_y_continuous(limits=c(1 - rel(0.5), height + rel(2)), expand=c(0, 0)) # Saving matrix plot if (A3) { width = 42 height = 29.7 dpi = 300 } else { width = 29.7 height = 21 dpi = 100 } # print('ggg') ggsave(plot=mat, path=outdirTmp, filename=paste(outnameTmp, '_', fL, imat, '.pdf', sep=''), width=width, height=height, units='cm', dpi=dpi) # print('hhh') } } } map_panel = function (list_df2plot, df_meta, computer_data_path, fr_shpdir, fr_shpname, bs_shpdir, bs_shpname, rv_shpdir, rv_shpname, idPer=1, outdirTmp='', codeLight=NULL, margin=NULL, showSea=TRUE) { fr_shppath = file.path(computer_data_path, fr_shpdir, fr_shpname) rv_shppath = file.path(computer_data_path, rv_shpdir, rv_shpname)
1471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540
bs_shppath = file.path(computer_data_path, bs_shpdir, bs_shpname) # France fr_spdf = readOGR(dsn=fr_shppath, verbose=FALSE) proj4string(fr_spdf) = CRS("+proj=longlat +ellps=WGS84") # Trasnformation en Lambert93 france = spTransform(fr_spdf, CRS("+init=epsg:2154")) df_france = tibble(fortify(france)) # Bassin hydrographique bassin = readOGR(dsn=bs_shppath, verbose=FALSE) df_bassin = tibble(fortify(bassin)) # Rseau hydrographique # river = readOGR(dsn=rv_shppath, verbose=FALSE) ### trop long ### # river = river[which(river$Classe == 1),] # df_river = tibble(fortify(river)) nbp = length(list_df2plot) # Get all different stations code Code = levels(factor(df_meta$code)) nCode = length(Code) df_trend = list_df2plot[[1]]$trend nPeriod_max = 0 for (code in Code) { df_trend_code = df_trend[df_trend$code == code,] Start = df_trend_code$period_start UStart = levels(factor(Start)) End = df_trend_code$period_end UEnd = levels(factor(End)) nPeriod = max(length(UStart), length(UEnd)) if (nPeriod > nPeriod_max) { nPeriod_max = nPeriod } } 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 (j in 1:nCode) { code = Code[j] df_trend_code = df_trend[df_trend$code == code,] Start = df_trend_code$period_start UStart = levels(factor(Start)) End = df_trend_code$period_end UEnd = levels(factor(End)) nPeriod = max(length(UStart), length(UEnd)) Periods = c() for (i in 1:nPeriod_max) { Periods = append(Periods, paste(substr(Start[i], 1, 4), substr(End[i], 1, 4), sep=' / '))
1541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610
} 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 (j in 1:nPeriod_max) { for (k in 1:nCode) { code = Code[k] for (i in 1:nbp) { df_data = list_df2plot[[i]]$data df_trend = list_df2plot[[i]]$trend p_threshold = list_df2plot[[i]]$p_threshold df_data_code = df_data[df_data$code == 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 TrendMean_code[j, i, k] = trendMean } } } 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 (i in 1:nbp) { if (i > 1 & !is.null(codeLight)) { break } outname = paste('map_', i, sep='') print(paste('map :', outname)) type = list_df2plot[[i]]$type map = ggplot() + theme_void() +
1611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680
# theme(plot.background=element_rect(fill=NA, # color="#EC4899")) + coord_fixed() + geom_polygon(data=df_france, aes(x=long, y=lat, group=group), color=NA, fill="grey97") + # geom_path(data=df_river, # aes(x=long, y=lat, group=group), # color="grey85", size=0.3) + geom_polygon(data=df_bassin, aes(x=long, y=lat, group=group), color="grey70", fill=NA, size=0.1) + geom_polygon(data=df_france, aes(x=long, y=lat, group=group), color="grey40", fill=NA, size=0.2) 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) 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=3) 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=3) } } 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) }
1681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750
lon = c() lat = c() fill = c() shape = c() trend = c() for (code in Code) { df_data = list_df2plot[[i]]$data df_trend = list_df2plot[[i]]$trend p_threshold = list_df2plot[[i]]$p_threshold df_data_code = df_data[df_data$code == 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 = color_res 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)
1751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820
shape = c(shape, shapetmp) trend = c(trend, trendMean) } 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() + annotate('text', x=-0.3, y=0.15, label=bquote(bold(.(type))), hjust=0, vjust=0, size=10, color="#00A3A8") + geom_line(aes(x=c(-0.3, 3.3), y=c(0.05, 0.05)), size=0.6, color="#00A3A8") + 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,
1821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890
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') # print(minTrendMean[idPer, i]) # print(maxTrendMean[idPer, i]) yTrend = (trend - minTrendMean[idPer, i]) / (maxTrendMean[idPer, i] - minTrendMean[idPer, i]) * valNorm xTrend = rnorm(length(yTrend), mean=1.75, sd=0.1)
1891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960
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 sec", 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) } hydrogramme = function (df_data_code, margin=NULL) { monthData = as.numeric(format(df_data_code$Date, "%m")) monthMean = c() for (i in 1:12) { data = df_data_code$Qm3s[monthData == i] monthMean[i] = mean(data, na.rm=TRUE) } monthNum = 1:12 monthName = c("jan", "fv", "mar", "avr", "mai", "jun", "jul", "ao", "sep", "oct", "nov", "dc") monthName = factor(monthName, levels=monthName) hyd = ggplot() + theme_ash +
1961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030
theme( # plot.background=element_rect(fill=NA, color="#EC4899"), panel.border=element_blank(), axis.text.x=element_text(margin=unit(c(0, 0, 0, 0), "mm"), angle=90, vjust=0.5, hjust=1), axis.ticks.x=element_blank(), axis.line.y=element_line(color='grey85', size=0.3)) if (is.null(margin)) { hyd = hyd + theme(plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm")) } else { hyd = hyd + theme(plot.margin=margin) } hyd = hyd + geom_bar(aes(x=monthNum, y=monthMean), stat='identity', fill="grey70", width=0.75, size=0.2) + scale_x_continuous(breaks=monthNum, labels=monthName, limits=c(0, max(monthNum)+0.5), expand=c(0, 0)) + scale_y_continuous(limits=c(0, max(monthMean)), expand=c(0, 0)) return (hyd) } info_panel = function(list_df2plot, df_meta, computer_data_path, fr_shpdir, fr_shpname, bs_shpdir, bs_shpname, rv_shpdir, rv_shpname, codeLight, df_data_code=NULL) { if (!is.null(df_data_code)) { hyd = hydrogramme(df_data_code, margin=margin(t=3, r=0, b=0, l=5, unit="mm")) } else { hyd = void } # databin = list_df2plot[[1]]$data[list_df2plot[[1]]$data$code == codeLight,] # yearLast = format(databin$Date[nrow(databin)], "%Y") # yearFirst = format(databin$Date[1], "%Y") # Nyear = as.numeric(yearLast) - as.numeric(yearFirst) + 1 map = map_panel(list_df2plot, df_meta, computer_data_path=computer_data_path, fr_shpdir=fr_shpdir, fr_shpname=fr_shpname, bs_shpdir=bs_shpdir, bs_shpname=bs_shpname, rv_shpdir=rv_shpdir, rv_shpname=rv_shpname, codeLight=codeLight, margin=margin(t=5, r=2, b=0, l=0, unit="mm"), showSea=FALSE) df_meta_code = df_meta[df_meta$code == codeLight,] nom = df_meta_code$nom nom = gsub("-", "-&nbsp;", nom) duration = as.numeric(format(as.Date(df_meta_code$fin), "%Y")) - as.numeric(format(as.Date(df_meta_code$debut), "%Y")) debut = format(as.Date(df_meta_code$debut), "%d/%m/%Y")
2031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100
fin = format(as.Date(df_meta_code$fin), "%d/%m/%Y") text1 = paste( "<b>", codeLight, '</b> - ', nom, sep='') text2 = paste( "<b>", "Gestionnaire : ", df_meta_code$gestionnaire, "<br>", "Rgion hydro : ", df_meta_code$region_hydro, "</b>", sep='') text3 = paste( "<b>", "Superficie : ", df_meta_code$surface_km2_BH, " [km<sup>2</sup>] <br>", "Altitude : ", df_meta_code$altitude_m_BH, " [m]<br>", "X = ", df_meta_code$L93X_m_BH, " [m ; Lambert 93]<br>", "Y = ", df_meta_code$L93Y_m_BH, " [m ; Lambert 93]", "</b>", sep='') text4 = paste( "<b>", "Date de dbut : ", debut, "<br>", "Date de fin : ", fin, "<br>", "Nombre d'annes : ", duration, " [ans]", "</b>", sep='') gtext1 = richtext_grob(text1, x=0, y=1, margin=unit(c(t=5, r=5, b=10, l=5), "mm"), hjust=0, vjust=1, gp=gpar(col="#00A3A8", fontsize=14)) gtext2 = richtext_grob(text2, x=0, y=1, margin=unit(c(t=0, r=0, b=0, l=5), "mm"), hjust=0, vjust=1, gp=gpar(col="grey20", fontsize=8)) gtext3 = richtext_grob(text3, x=0, y=1, margin=unit(c(t=0, r=0, b=0, l=5), "mm"), hjust=0, vjust=1, gp=gpar(col="grey20", fontsize=9)) gtext4 = richtext_grob(text4, x=0, y=1, margin=unit(c(t=0, r=0, b=0, l=0), "mm"), hjust=0, vjust=1, gp=gpar(col="grey20", fontsize=9)) void = void + theme(plot.background=element_rect(fill=NA, color="#EC4899"), plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm")) P = list(gtext1, gtext2, gtext3, gtext4, hyd, map) # P = list(void, void, void, void, void, void, void) LM = matrix(c(1, 1, 1, 6, 7, 7, 7, 6, 2, 2, 5, 6, # 2, 2, 5, 6, 3, 4, 5, 6, 3, 4, 5, 6), nrow=5, byrow=TRUE)
2101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170
heights = rep(1, times=nrow(LM)) heights[2] = 0.1 heights[3] = 0.8 plot = grid.arrange(grobs=P, layout_matrix=LM, heights=heights) return(plot) } histogram = function (data_bin, df_meta, figdir='', filedir_opt='') { # Get all different stations code Code = levels(factor(df_meta$code)) nCode = length(Code) # If there is not a dedicated figure directory it creats one outdir = file.path(figdir, filedir_opt, sep='') if (!(file.exists(outdir))) { dir.create(outdir) } datebreak = 10 dateminbreak = 1 res_hist = hist(data_bin, breaks='years', plot=FALSE) counts = res_hist$counts counts_pct = counts/nCode * 100 breaks = as.Date(res_hist$breaks) mids = as.Date(res_hist$mids) p = ggplot() + theme_ash + theme(panel.grid.major.y=element_line(color='grey85', size=0.15), axis.title.y=element_blank()) + geom_bar(aes(x=mids, y=counts_pct), stat='identity', fill="#00A3A8") + scale_x_date(date_breaks=paste(as.character(datebreak), 'year', sep=' '), date_minor_breaks=paste(as.character(dateminbreak), 'year', sep=' '), guide='axis_minor', date_labels="%Y", limits=c(min(data_bin)-years(0), max(data_bin)+years(0)), expand=c(0, 0)) + scale_y_continuous(limits=c(0, max(counts_pct)*1.1), expand=c(0, 0)) ggsave(plot=p, path=outdir, filename=paste('hist_break_date', '.pdf', sep=''), width=10, height=10, units='cm', dpi=100) } cumulative = function (data_bin, df_meta, dyear=10, figdir='', filedir_opt='') { # Get all different stations code Code = levels(factor(df_meta$code)) nCode = length(Code)
2171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204220522062207220822092210221122122213221422152216221722182219222022212222222322242225222622272228222922302231223222332234223522362237223822392240
# If there is not a dedicated figure directory it creats one outdir = file.path(figdir, filedir_opt, sep='') if (!(file.exists(outdir))) { dir.create(outdir) } datebreak = 10 dateminbreak = 1 res_hist = hist(data_bin, breaks='years', plot=FALSE) counts = res_hist$counts cumul = cumsum(counts) cumul_pct = cumul/nCode * 100 breaks = as.Date(res_hist$breaks) mids = as.Date(res_hist$mids) mids = c(mids[1] - years(dyear), mids[1] - years(1), mids, mids[length(mids)] + years(dyear)) cumul_pct = c(0, 0, cumul_pct, cumul_pct[length(cumul_pct)]) mids = mids + months(6) breaks = breaks + 1 breaks = breaks[-length(breaks)] DB = c() for (i in 1:length(breaks)) { DB = c(DB, rep(breaks[i], times=counts[i])) } q50 = as.Date(quantile(DB, probs=0.5)) + years(1) print(paste('mediane :', q50)) p = ggplot() + theme_ash + theme(panel.grid.major.y=element_line(color='grey85', size=0.15), axis.title.y=element_blank()) + geom_line(aes(x=mids, y=cumul_pct), color="#00A3A8") + geom_line(aes(x=c(q50, q50), y=c(0, 100)), color="wheat", lty='dashed') + scale_x_date(date_breaks=paste(as.character(datebreak), 'year', sep=' '), date_minor_breaks=paste(as.character(dateminbreak), 'year', sep=' '), guide='axis_minor', date_labels="%Y", limits=c(min(mids)-years(0), max(mids)+years(0)), expand=c(0, 0)) + scale_y_continuous(limits=c(-1, 101), expand=c(0, 0)) ggsave(plot=p, path=outdir, filename=paste('cumul_break_date', '.pdf', sep=''), width=10, height=10, units='cm', dpi=100) } # get_color = function (value, min, max, ncolor=256, palette_name='perso', reverse=FALSE, nbTick=10) {
2241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310
# if (palette_name == 'perso') { # colorList = c('#0f3b57', # '#1d7881', # '#80c4a9', # '#e2dac6', #mid # '#fadfad', # '#d08363', # '#7e392f') # } else { # colorList = brewer.pal(11, palette_name) # } # nSample = length(colorList) # palette = colorRampPalette(colorList)(ncolor) # Sample_hot = 1:(as.integer(nSample/2)+1) # Sample_cold = (as.integer(nSample/2)+1):nSample # palette_hot = colorRampPalette(colorList[Sample_hot])(ncolor) # palette_cold = colorRampPalette(colorList[Sample_cold])(ncolor) # if (reverse) { # palette = rev(palette) # palette_hot = rev(palette_hot) # palette_cold = rev(palette_cold) # } # if (value < 0) { # idNorm = (value - min) / (0 - min) # id = round(idNorm*(ncolor - 1) + 1, 0) # color = palette_cold[id] # } else { # idNorm = (value - 0) / (max - 0) # id = round(idNorm*(ncolor - 1) + 1, 0) # color = palette_hot[id] # } # if (min < 0 & max < 0) { # paletteShow = palette_cold # idTick = c() # for (i in 1:nbTick) { # id = round((ncolor-1)/(nbTick-1)*(i-1)) + 1 # idTick = c(idTick, id) # } # labTick = seq(min, max, length.out=nbTick) # colTick = paletteShow[idTick] # } else if (min > 0 & max > 0) { # paletteShow = palette_hot # idTick = c() # for (i in 1:nbTick) { # id = round((ncolor-1)/(nbTick-1)*(i-1)) + 1 # idTick = c(idTick, id) # } # labTick = seq(min, max, length.out=nbTick) # colTick = paletteShow[idTick] # } else { # paletteShow = palette # nbSemiTick = round(nbTick/2) + 1 # idSemiTick = c() # for (i in 1:nbSemiTick) { # id = round((ncolor-1)/(nbSemiTick-1)*(i-1)) + 1 # idSemiTick = c(idSemiTick, id) # } # labTick_hot = seq(0, max, length.out=nbSemiTick) # labTick_cold = seq(min, 0, length.out=nbSemiTick)
2311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380
# colTick_hot = palette_hot[idSemiTick] # colTick_cold = palette_cold[idSemiTick] # idTick = as.integer(seq(1, ncolor, length.out=nbTick+1)) # labTick = c(labTick_cold, labTick_hot[-1]) # colTick = c(colTick_cold, colTick_hot[-1]) # } # return(list(color=color, palette=paletteShow, # idTick=idTick, labTick=labTick, colTick=colTick)) # } get_color = function (value, min, max, ncolor=256, palette_name='perso', reverse=FALSE) { if (palette_name == 'perso') { colorList = c('#0f3b57', '#1d7881', '#80c4a9', '#e2dac6', #mid '#fadfad', '#d08363', '#7e392f') } else { colorList = brewer.pal(11, palette_name) } nSample = length(colorList) palette = colorRampPalette(colorList)(ncolor) Sample_hot = 1:(as.integer(nSample/2)+1) Sample_cold = (as.integer(nSample/2)+1):nSample palette_hot = colorRampPalette(colorList[Sample_hot])(ncolor) palette_cold = colorRampPalette(colorList[Sample_cold])(ncolor) if (reverse) { palette = rev(palette) palette_hot = rev(palette_hot) palette_cold = rev(palette_cold) } if (is.na(value)) { return (NA) } if (value < 0) { idNorm = (value - min) / (0 - min) id = round(idNorm*(ncolor - 1) + 1, 0) color = palette_cold[id] } else { idNorm = (value - 0) / (max - 0) id = round(idNorm*(ncolor - 1) + 1, 0) color = palette_hot[id] } return(color) } get_palette = function (min, max, ncolor=256, palette_name='perso', reverse=FALSE, nbTick=10) { if (palette_name == 'perso') { colorList = c('#0f3b57', '#1d7881', '#80c4a9', '#e2dac6', #mid '#fadfad', '#d08363', '#7e392f')
2381238223832384238523862387238823892390239123922393239423952396239723982399240024012402240324042405240624072408240924102411241224132414241524162417241824192420242124222423242424252426242724282429243024312432243324342435243624372438243924402441244224432444244524462447244824492450
} else { colorList = brewer.pal(11, palette_name) } nSample = length(colorList) palette = colorRampPalette(colorList)(ncolor) Sample_hot = 1:(as.integer(nSample/2)+1) Sample_cold = (as.integer(nSample/2)+1):nSample palette_hot = colorRampPalette(colorList[Sample_hot])(ncolor) palette_cold = colorRampPalette(colorList[Sample_cold])(ncolor) if (reverse) { palette = rev(palette) palette_hot = rev(palette_hot) palette_cold = rev(palette_cold) } if (min < 0 & max < 0) { paletteShow = palette_cold } else if (min > 0 & max > 0) { paletteShow = palette_hot } else { paletteShow = palette } posTick = seq(0, 1, length.out=nbTick) labTick = c() colTick = c() for (i in 1:nbTick) { lab = (i-1)/(nbTick-1) * (max - min) + min labTick = c(labTick, lab) col = get_color(lab, min=min, max=max, ncolor=ncolor, palette_name=palette_name, reverse=reverse) colTick = c(colTick, col) } return(list(palette=paletteShow, posTick=posTick, labTick=labTick, colTick=colTick)) } void = ggplot() + geom_blank(aes(1,1)) + theme( plot.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), axis.line = element_blank() ) palette_tester = function (n=256) { X = 1:n Y = rep(0, times=n) palette = colorRampPalette(c(
2451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520
# '#1a4157', # '#00af9d', # '#fbdd7e', # '#fdb147', # '#fd4659' # '#543005', # '#8c510a', # '#bf812d', # '#dfc27d', # '#f6e8c3', # '#f5f5f5', # '#c7eae5', # '#80cdc1', # '#35978f', # '#01665e', # '#003c30' '#0f3b57', '#1d7881', '#80c4a9', '#e2dac6', #mid '#fadfad', '#d08363', '#7e392f' # '#193830', # '#2A6863', # '#449C93', # '#7ACEB9', # '#BCE6DB', # '#EFE0B0', # '#D4B86A', # '#B3762A', # '#7F4A23', # '#452C1A' ))(n) p = ggplot() + theme( plot.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), axis.line = element_blank() ) + geom_line(aes(x=X, y=Y), color=palette[X], size=60) + scale_y_continuous(expand=c(0, 0)) ggsave(plot=p, filename=paste('palette_test', '.pdf', sep=''), width=10, height=10, units='cm', dpi=100) } # palette_tester() get_power = function (value) { if (value > 1) { power = nchar(as.character(as.integer(value))) - 1
252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559
} else { dec = gsub('0.', '', as.character(value), fixed=TRUE) ndec = nchar(dec) nnum = nchar(as.character(as.numeric(dec))) power = -(ndec - nnum + 1) } return(power) } gg_circle = function(r, xc, yc, color="black", fill=NA, ...) { x = xc + r*cos(seq(0, pi, length.out=100)) ymax = yc + r*sin(seq(0, pi, length.out=100)) ymin = yc + r*sin(seq(0, -pi, length.out=100)) annotate("ribbon", x=x, ymin=ymin, ymax=ymax, color=color, fill=fill, ...) } gpct = function (pct, L, ref=NULL, shift=FALSE) { if (is.null(ref)) { minL = min(L, na.rm=TRUE) } else { minL = ref } maxL = max(L, na.rm=TRUE) spanL = maxL - minL xL = pct/100 * as.numeric(spanL) if (shift) { xL = xL + minL } return (xL) }