Newer
Older
#
# *1 INRAE, France
# louis.heraut@inrae.fr
#
# This file is part of ash R toolbox.
#
# ash R toolbox is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# ash R toolbox is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with ash R toolbox. If not, see <https://www.gnu.org/licenses/>.
# ///
#
#
# plotting/matrix.R
#
# Allows the creation of a summarizing matrix of trend and break analyses
# Generates a summarizing matrix of the trend analyses of all station for different hydrological variables and periods. Also shows difference of means between specific periods.
matrix_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
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)
# Fix the maximal number of period to the minimal possible
# Extract start and end of trend periods
Start = df_trend_code$period_start
# Get the name of the different period
UStart = levels(factor(Start))
# Compute the max of different start and end
# so the number of different period
# 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) {
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
TrendValue_code = array(rep(1, nPeriod_trend*nbp*nCode),
# Extracts the data corresponding to the
# current variable
# Extracts the trend corresponding to the
# current variable
# Extracts the type of the variable
type = list_df2plot[[i]]$type
# Extracts the data corresponding to the code
df_data_code = df_data[df_data$code == code,]
# Extracts the trend corresponding to the code
Start = tab_Start[k, i, j]
End = tab_End[k, i, j]
Periods = tab_Periods[k, i, 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,]
# If there is more than one trend on the same period
# If it is a flow variable
if (type == 'sévérité') {
# Normalises the trend value by the mean of the data
trendValue = df_trend_code_per$trend / dataMean
# If it is a date variable
} else if (type == 'saisonnalité') {
# Just stocks the trend value
trendValue = df_trend_code_per$trend
}
# Computes the min and the max of the mean trend for
# all the station
minTrendValue = apply(TrendValue_code, c(1, 2), min, na.rm=TRUE)
maxTrendValue = apply(TrendValue_code, c(1, 2), max, na.rm=TRUE)
DataMean_trend = c()
Fill_trend = c()
Color_trend = c()
# Extracts the data corresponding to the current variable
# Extracts the trend corresponding to the
# current variable
df_trend = list_df2plot[[i]]$trend
p_threshold = list_df2plot[[i]]$p_threshold
# Extract the type of the variable to plot
type = list_df2plot[[i]]$type
# Extracts the data corresponding to the code
df_data_code = df_data[df_data$code == code,]
# Extracts the trend corresponding to the code
Start = tab_Start[k, i, j]
End = tab_End[k, i, j]
Periods = tab_Periods[k, i, 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,]
# If there is more than one trend on the same period
# If it is a flow variable
if (type == 'sévérité') {
# Normalises the trend value by the mean of the data
trendValue = df_trend_code_per$trend / dataMean
# If it is a date variable
} else if (type == 'saisonnalité') {
# Just stocks the trend value
trendValue = df_trend_code_per$trend
}
color_res = get_color(trendValue,
minTrendValue[j, i],
maxTrendValue[j, i],
# Specifies the color fill and contour of
# table cells
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)
Code_trend = append(Code_trend, code)
Pthresold_trend = append(Pthresold_trend, Pthresold)
TrendValue_trend = append(TrendValue_trend, trendValue)
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)) {
# Blank vectors to store info about breaking analysis
# Convert 'mean_period' to list
mean_period = as.list(mean_period)
# Number of mean period
nPeriod_mean = length(mean_period)
# Blank array to store difference of mean between two periods
breakValue_code = array(rep(1, nPeriod_mean*nbp*nCode),
# Blank array to store mean for a temporary period in order
# to compute the difference of mean with a second period
dataMeantmp = array(rep(NA, nbp*nCode),
dim=c(nbp, nCode))
# Extracts the data corresponding to
# the current variable
# Extract the type of the variable to plot
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,]
Datemin = min(df_data_code_per$Date)
Datemax = max(df_data_code_per$Date)
Periods = paste(Datemin, Datemax,
sep=' / ')
# Mean of the flow over the sub period
# If it is a flow variable
if (type == 'sévérité') {
# Normalises the break by the mean of the
# initial period
breakValue = Break / dataMeantmp[i, k]
# If it is a date variable
} else if (type == 'saisonnalité') {
# Just stocks the break value
breakValue = Break
}
# Stores temporarily the mean of the current period
Periods_mean = append(Periods_mean, Periods)
NPeriod_mean = append(NPeriod_mean, j)
Code_mean = append(Code_mean, code)
DataMean_mean = append(DataMean_mean, dataMean)
breakValue_mean = append(breakValue_mean,
breakValue)
# Computes the min and the max of the averaged trend for
color_res = get_color(breakValue,
minBreakValue[j, i],
maxBreakValue[j, i],
Fill_mean = append(Fill_mean, fill)
Color_mean = append(Color_mean, color)
# If the slice option is not specified, the info for all
# stations will be draw on the same page
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
# Gets all the different type of plots
Type = levels(factor(Type_trend))
nbType = length(Type)
# For all the type of plots
for (itype in 1:nbType) {
# Gets the type
type = Type[itype]
# Extracts each possibilities of first letter of station code
firstLetter = levels(factor(substr(Code, 1, 1)))
# Number of different first letters
nfL = length(firstLetter)
# For all the available first letter
for (ifL in 1:nfL) {
# Gets the first letter
fL = firstLetter[ifL]
# Get only station code with the same first letter
subCodefL = Code[substr(Code, 1, 1) == fL]
# Counts the number of station in it
nsubCodefL = length(subCodefL)
# Computes the number of pages needed to plot all stations
nMat = as.integer(nsubCodefL/slice) + 1
# For all the pages
for (iMat in 1:nMat) {
# Print the matrix name
print(paste('Matrix ', iMat, '/', nMat,
' of ', type,
' for region : ', fL,
" (",
round((ifL + nfL*(itype-1)) / (nfL*2) * 100,
# Extracts the station for the current page
subCode = subCodefL[(slice*(iMat-1)+1):(slice*iMat)]
# Removes NA stations
subCode = subCode[!is.na(subCode)]
# Reverses verticale order of stations
subCode = rev(subCode)
# Gets the number of station for the page
nsubCode = length(subCode)
# Creates logical vector to select only info about
# stations that will be plot on the page
CodefL_trend =
Code_trend %in% subCode & Type_trend == type
# Extracts those info
subPeriods_trend = Periods_trend[CodefL_trend]
subNPeriod_trend = NPeriod_trend[CodefL_trend]
subVar_trend = Var_trend[CodefL_trend]
subType_trend = Type_trend[CodefL_trend]
subCode_trend = Code_trend[CodefL_trend]
subPthresold_trend = Pthresold_trend[CodefL_trend]
subTrendValue_trend = TrendValue_trend[CodefL_trend]
subDataMean_trend = DataMean_trend[CodefL_trend]
subFill_trend = Fill_trend[CodefL_trend]
subColor_trend = Color_trend[CodefL_trend]
# Same for breaking analysis
CodefL_mean =
Code_mean %in% subCode & Type_mean == type
# Extracts right info
subPeriods_mean = Periods_mean[CodefL_mean]
subNPeriod_mean = NPeriod_mean[CodefL_mean]
subVar_mean = Var_mean[CodefL_mean]
subType_mean = Type_mean[CodefL_mean]
subCode_mean = Code_mean[CodefL_mean]
subDataMean_mean = DataMean_mean[CodefL_mean]
subbreakValue_mean = breakValue_mean[CodefL_mean]
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
subFill_mean = Fill_mean[CodefL_mean]
subColor_mean = Color_mean[CodefL_mean]
# Gets the number of variable to plot in
# function of the current type
nbpMod =
length(levels(factor(subVar_trend)))
### Plot ###
# Fixes the height and width of the table according to
# the number of station and the number of column to draw
height = nsubCode
width = nbpMod * 2 * nPeriod_trend + nPeriod_trend + nPeriod_mean * nbpMod + nPeriod_mean + nbpMod
# Fixes the size of the plot area to keep proportion right
options(repr.plot.width=width, repr.plot.height=height)
# Open a new plot with a personalise theme
mat = ggplot() + theme_ash +
# Modification of theme in order to remove axis
theme(
panel.border=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
plot.margin=margin(t=5, r=5, b=5, l=5, unit="mm")
)
# Extracts the name of the currently hydrological
# region plotted
title = df_meta[df_meta$code == subCode[1],]$region_hydro
sep='')
# Postion and name of the title
xt = 1 - 6
yt = height + 2
annotate("text", x=xt, y=yt,
label=Title,
hjust=0, vjust=1,
size=6, color="#00A3A8")
### Trend ###
# For all the trend period
for (j in 1:nPeriod_trend) {
# Extracts the info to plot associated to the
# right period
Periods_trend_per =
subPeriods_trend[subNPeriod_trend == j]
NPeriods_trend_per =
subNPeriod_trend[subNPeriod_trend == j]
Var_trend_per =
subVar_trend[subNPeriod_trend == j]
Type_trend_per =
subType_trend[subNPeriod_trend == j]
Code_trend_per =
subCode_trend[subNPeriod_trend == j]
Pthresold_trend_per =
subPthresold_trend[subNPeriod_trend == j]
TrendValue_trend_per =
subTrendValue_trend[subNPeriod_trend == j]
DataMean_trend_per =
subDataMean_trend[subNPeriod_trend == j]
Fill_trend_per =
subFill_trend[subNPeriod_trend == j]
Color_trend_per =
subColor_trend[subNPeriod_trend == j]
# Converts the variable list into levels for factor
# Converts the vector of hydrological variable to
# a vector of integer associated to those variable
Xtmp = as.integer(factor(as.character(Var_trend_per),
levels=levels))
# Computes X position of the column for
# the period dates
Xc = j + (j - 1)*nbpMod*2
# Computes X positions of columns for
# the mean of variables
Xm = Xtmp + (j - 1)*nbpMod*2 + j
# Computes X positions of columns for
# the averaged trend
X = Xtmp + (j - 1)*nbpMod*2 + nbpMod + j
# Computes Y positions of each line for each station
Y = as.integer(factor(Code_trend_per))
# Reverses vertical order of stations
Y = rev(Y)
# Position of a line to delimite periods
x = Xc - 0.4
xend = X[length(X)] + 0.4
y = height + 1.1
yend = height + 1.1
# Drawing of the line
mat = mat +
annotate("segment",
x=x, xend=xend,
y=y, yend=yend,
color="grey40", size=0.35)
# Position of the name of the current period
yt = y + 0.15
Start = trend_period[[j]][1]
End = trend_period[[j]][2]
# Name of the period
periodName =
bquote(bold('Période')~bold(.(as.character(j))))
# Naming the period
annotate("text", x=x, y=yt,
label=periodName,
hjust=0, vjust=0.5,
# For all the variable
for (i in 1:length(X)) {
mat = mat +
# Plots circles for averaged trends
gg_circle(r=0.45, xc=X[i], yc=Y[i],
fill=Fill_trend_per[i],
color=Color_trend_per[i]) +
# Plots circles for averaged of variables
gg_circle(r=0.45, xc=Xm[i], yc=Y[i],
fill='white', color='grey40') +
# Plots circles for the column of period dates
gg_circle(r=0.45, xc=Xc, yc=Y[i],
fill='white', color='grey40')
}
trendValue = TrendValue_trend_per[i]
type = Type_trend_per[i]
# If it is a flow variable
if (type == 'sévérité') {
Nsign_mean = 2
# Converts it to the right format with
# two significant figures
trendValueC = signif(trendValue*100, 2)
# If it is a date variable
} else if (type == 'saisonnalité') {
# Fixes the significants number for mean to 3
Nsign_mean = 3
# Converts the trend value with two
# significant figures
trendValueC = signif(trendValue, 2)
}
# If it is significative
if (!is.na(Pthresold_trend_per[i])) {
# The text color is white
Tcolor = 'white'
# Otherwise
} else {
# The text is grey
Tcolor = 'grey85'
}
# Same for averaged variables over
# the current period
dataMean = DataMean_trend_per[i]
mat = mat +
# Writes the mean trend
annotate('text', x=X[i], y=Y[i],
hjust=0.5, vjust=0.5,
size=3, color=Tcolor) +
# Writes the mean of the associated variable
annotate('text', x=Xm[i], y=Y[i],
label=dataMeanC,
hjust=0.5, vjust=0.5,
size=3, color='grey40')
}
annotate('text', x=Xc, y=max(Y) + 0.9,
label=bquote(bold('Début')),
size=3, color='grey20') +
annotate('text', x=Xc, y=max(Y) + 0.63,
label=bquote(bold('Fin')),
size=3, color='grey20')
# For all variable
for (i in 1:nbpMod) {
# Extract the variable of the plot
var = subVar_trend[i]
type = subType_trend[i]
# If it is a flow variable
if (type == 'sévérité') {
# Fixes the unit of the mean and the trend
# for the flow
unit_mean = bquote('['*m^3*'.'*s^{-1}*']')
unit_trend = bquote('[%.'*an^{-1}*']')
# If it is a date variable
} else if (type == 'saisonnalité') {
# Fixes the unit of the mean and the trend
# for the date
unit_mean = bquote('[jour]')
unit_trend = bquote('[jour.'*an^{-1}*']')
}
mat = mat +
# Writes the unit of the variable
annotate('text', x=X[i], y=max(Y) + 0.63,
hjust=0.5, vjust=0.5,
size=2, color='grey40') +
# Writes the type of the variable
annotate('text', x=X[i], y=max(Y) + 0.9,
label=bquote(.(var)),
hjust=0.5, vjust=0.5,
size=3.25, color='grey20') +
# Writes the unit of the averaged variable
annotate('text', x=Xm[i], y=max(Y) + 0.63,
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
hjust=0.5, vjust=0.5,
size=2, color='grey40') +
# Writes the type of the averaged variable
annotate('text', x=Xm[i], y=max(Y) + 0.9,
label=expr(bar(!!var)),
hjust=0.5, vjust=0.5,
size=3.25, color='grey20')
}
# For all the station on the page
for (k in 1:nsubCode) {
# Gets the code
code = subCode[k]
# Extracts label for the period dates
label =
Periods_trend_per[Code_trend_per == code][1]
# Gets the start and end of the period
# for the station
periodStart = substr(label, 1, 4)
periodEnd = substr(label, 14, 17)
mat = mat +
# Writes the starting value
annotate('text', x=Xc, y=k + 0.13,
label=bquote(bold(.(periodStart))),
hjust=0.5, vjust=0.5,
size=3, color='grey40') +
# Writes the ending value
annotate('text', x=Xc, y=k - 0.13,
label=bquote(bold(.(periodEnd))),
hjust=0.5, vjust=0.5,
size=3, color='grey40')
}
### Mean ###
# For all the trend period
for (j in 1:nPeriod_mean) {
# Extracts the info to plot associated to the
# right period
Periods_mean_per =
subPeriods_mean[subNPeriod_mean == j]
NPeriods_mean_per =
subNPeriod_mean[subNPeriod_mean == j]
Var_mean_per =
subVar_mean[subNPeriod_mean == j]
Type_mean_per =
subType_mean[subNPeriod_mean == j]
Code_mean_per =
subCode_mean[subNPeriod_mean == j]
DataMean_mean_per =
subDataMean_mean[subNPeriod_mean == j]
breakValue_mean_per =
subbreakValue_mean[subNPeriod_mean == j]
Fill_mean_per =
subFill_mean[subNPeriod_mean == j]
Color_mean_per =
subColor_mean[subNPeriod_mean == j]
# Converts the variable list into levels for factor
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
# Converts the vector of hydrological variable to
# a vector of integer associated to those variable
Xtmp_mean =
as.integer(factor(as.character(Var_mean_per),
levels=levels))
# Computes X position of the column for
# the period dates
Xc_mean = j + (j - 1)*nbpMod + X[length(X)]
# Computes X positions of columns for
# the mean of variables
Xm_mean =
Xtmp_mean + (j - 1)*nbpMod + j + X[length(X)]
# Computes X positions of columns for
# the difference of mean between periods (break)
Xr_mean =
Xtmp_mean + (j - 1)*nbpMod*2 + j + X[length(X)]
# Computes Y positions of each line for each station
Y_mean = as.integer(factor(Code_mean_per))
# Reverses vertical order of stations
Y_mean = rev(Y_mean)
# Position of a line to delimite periods
x = Xc_mean - 0.4
xend = Xm_mean[length(Xm_mean)] + 0.25
y = height + 1.1
yend = height + 1.1
mat = mat +
annotate("segment",
x=x, xend=xend,
y=y, yend=yend,
color="grey40", size=0.35)
# Position of the name of the current period
yt = y + 0.15
Start = mean_period[[j]][1]
End = mean_period[[j]][2]
# Name of the period
periodName = bquote(bold('Période')~bold(.(as.character(j+nPeriod_trend))))
# Naming the period
# Position of a line to delimite results of
# difference of mean bewteen periods
x = Xr_mean[1] - 0.4
xend = Xr_mean[length(Xr_mean)] + 0.25
# Drawing of the line
annotate("segment",
x=x, xend=xend,
y=y, yend=yend,
color="grey40", size=0.35)
# Naming the breaking columns
breakName = bquote(bold('Écart')~bold(.(as.character(j-1+nPeriod_trend)))*bold('-')*bold(.(as.character(j+nPeriod_trend))))
# Writes the name
mat = mat +
annotate("text", x=x, y=yt,
label=breakName,
hjust=0, vjust=0.5,
size=3, color='grey40')
# For all the variable
for (i in 1:length(Xm_mean)) {
mat = mat +
# Plots circles for averaged variables
gg_circle(r=0.45, xc=Xm_mean[i], yc=Y[i],
fill='white', color='grey40') +
# Plots circles for the column of period dates
gg_circle(r=0.45, xc=Xc_mean, yc=Y[i],
fill='white', color='grey40')
# If this is not the first period
if (j > 1) {
mat = mat +
# Plots circles for breaking results
gg_circle(r=0.45, xc=Xr_mean[i], yc=Y[i],
fill=Fill_mean_per[i],
color=Color_mean_per[i])
}
}
# For all averaged variables on this period
for (i in 1:length(DataMean_mean_per)) {
type = Type_mean_per[i]
# If it is a flow variable
if (type == 'sévérité') {
# The number of significant figures for
# flow mean is 2
Nsign_mean = 2
# If it is a date variable
} else if (type == 'saisonnalité') {
# The number of significant figures for
# date mean is 3
Nsign_mean = 3
}
# Extracts values of averaged variables
dataMean = DataMean_mean_per[i]
# Converts it to the right format with two
# significant figures
size=3, color='grey40')
# If this is not the first period
if (j > 1) {
# Extracts values of breaking between periods
breakValue = breakValue_mean_per[i]
# If it is a flow variable
if (type == 'sévérité') {
# Converts it to the right format with two
# significant figures
breakValueC = signif(breakValue*100, 2)
# If it is a date variable
} else if (type == 'saisonnalité') {
# Converts the break value with two
# significant figures
breakValueC = signif(breakValue, 2)
}
# Writes breaking values
mat = mat +
annotate('text', x=Xr_mean[i], y=Y[i],
annotate('text', x=Xc_mean, y=max(Y) + 0.9,
label=bquote(bold('Début')),
size=3, color='grey20') +
annotate('text', x=Xc_mean, y=max(Y) + 0.63,
label=bquote(bold('Fin')),
size=3, color='grey20')
# For all variables
for (i in 1:nbpMod) {
# Extract the variable of the plot
var = subVar_mean[i]
type = subType_mean[i]
# If it is a flow variable
if (type == 'sévérité') {
# Fixes the unit of the mean and the break
# for the flow
unit_mean = bquote('['*m^3*'.'*s^{-1}*']')
unit_break = bquote('[%]')
# If it is a date variable
# Fixes the unit of the mean and the break
# for the date
} else if (type == 'saisonnalité') {
unit_mean = bquote('[jour]')
unit_break = bquote('[jour]')
}
# Writes the unit of the averaged variable
annotate('text',
x=Xm_mean[i], y=max(Y) + 0.63,
# Writes the type of the averaged variable
annotate('text',
x=Xm_mean[i], y=max(Y) + 0.9,
label=expr(bar(!!var)),
hjust=0.5, vjust=0.5,