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/map.R
#
# Deals with the creation of a map for presenting the trend analysis of hydrological variables
## 1. MAP PANEL ______________________________________________________
# Generates a map plot of the tendancy of a hydrological variable
map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1,
trend_period, mean_period, colorForce=FALSE,
outdirTmp='', codeLight=NULL,
df_france = df_shapefile$france
df_bassin = df_shapefile$bassin
nbp = length(list_df2plot)
# Get all different stations code
Code = levels(factor(df_meta$code))
nCode = length(Code)
if (!is.null(trend_period)) {
# Convert 'trend_period' to list
trend_period = as.list(trend_period)
# Number of trend period
nPeriod_trend = length(trend_period)
# Extracts the min and the max of the mean trend for all the station
res = short_trendExtremes(list_df2plot, Code, nPeriod_trend, nbp, nCode, colorForce)
# 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)
res = short_meanExtremes(list_df2plot, Code, nPeriod_mean, nbp, nCode)
minBreakValue = res$min
maxBreakValue = res$max
breakValue_code = res$value
} else {
for (j in 1:nPeriod_mean) {
# For all variable
for (i in 1:nbp) {
# If there is a specified station code to highlight (mini map)
# and there has already been one loop
if ((i > 1 | j > 1) & !is.null(codeLight)) {
# Stop the for loop over the variable
break
}
# Extracts the variable of the plot
var = list_df2plot[[i]]$var
# Extracts the type of variable of the plot
type = list_df2plot[[i]]$type
# Explanations about the variable
glose = list_df2plot[[i]]$glose
# Createsa name for the map
if (j > 1) {
outname = paste('map_d', var, sep='')
} else {
outname = paste('map_', var, sep='')
}
if (j > 1) {
mapName = 'difference'
} else {
mapName = 'tendence'
}
print(paste('Map of ', mapName, ' for : ', var,
" (",
round(n_loop/N_loop*100, 0),
" %)",
# If there is no specified station code to highlight
# (mini map)
if (is.null(codeLight)) {
# Sets the size of the countour
sizefr = 0.45
sizebs = 0.4
sizerv = 0.3
} else {
sizefr = 0.35
sizebs = 0.3
sizerv = 0.2
}
# Stores the coordonate system
cf = coord_fixed()
# Makes it the default one to remove useless warning
cf$default = TRUE
# Open a new plot with the personalise theme
map = ggplot() + theme_void() +
# theme(
# # Ticks marker
# axis.ticks.x=element_line(color='grey75', size=0.3),
# axis.ticks.y=element_line(color='grey75', size=0.3),
# # Ticks label
# axis.text.x=element_text(color='grey40'),
# axis.text.y=element_text(color='grey40')) +
# Fixed coordinate system (remove useless warning)
cf +
# Plot the background of France
geom_polygon(data=df_france,
aes(x=long, y=lat, group=group),
color=NA, fill="grey97")
# If the river shapefile exists
if (!is.null(df_river)) {
# Plot the river
map = map +
geom_path(data=df_river,
aes(x=long, y=lat, group=group),
color="grey85", size=sizerv)
}
# Plot the hydrological basin
geom_polygon(data=df_bassin,
aes(x=long, y=lat, group=group),
color="grey70", fill=NA, size=sizebs) +
# Plot the hydrological sub-basin
geom_polygon(data=df_subbassin,
aes(x=long, y=lat, group=group),
color="grey70", fill=NA, size=sizebs) +
# Plot the countour of France
geom_polygon(data=df_france,
aes(x=long, y=lat, group=group),
color="grey40", fill=NA, size=sizefr)
if (is.null(codeLight)) {
xBasin = c(410000, 520000, 630000,
620000, 510000, 450000,
390000)
yBasin = c(6280000, 6290000, 6320000,
6385000, 6450000, 6530000,
nameBasin = c('Adour', 'Garonne', 'Tarn-Aveyron',
'Lot', 'Dordogne', 'Charente',
'Fleuves-\nCôtiers')
nBasin = length(xBasin)
plot_basin = tibble(x=xBasin, y=yBasin, label=nameBasin)
map = map +
geom_shadowtext(data=plot_basin,
aes(x=x, y=y, label=label),
fontface="bold",
color="grey85",
bg.colour="grey97",
hjust=0.5, vjust=0.5, size=5)
}
# If the sea needs to be shown
if (showSea) {
# Leaves space around the France
xlim = c(295000, 790000)
ylim = c(6125000, 6600000)
# Otherwise
} else {
# Leaves minimal space around France
xlim = c(305000, 790000)
ylim = c(6135000, 6600000)
}
# If there is no specified station code to highlight (mini map)
if (is.null(codeLight)) {
# Sets a legend scale start
xmin = gpct(4, xlim, shift=TRUE)
# Sets graduations
xint = c(0, 10*1E3, 50*1E3, 100*1E3)
# Sets the y postion
ymin = gpct(5, ylim, shift=TRUE)
# Sets the height of graduations
ymax = ymin + gpct(1, ylim)
# Size of the value
size = 3
# Size of the 'km' unit
sizekm = 2.5
# If there is a specified station code
} else {
# Same but with less graduation and smaller size
xmin = gpct(2, xlim, shift=TRUE)
xint = c(0, 100*1E3)
ymin = gpct(1, ylim, shift=TRUE)
ymax = ymin + gpct(3, ylim)
size = 2
sizekm = 1.5
}
# Adds the base line of the scale
geom_line(aes(x=c(xmin, max(xint)+xmin),
y=c(ymin, ymin)),
color="grey40", size=0.2) +
# Adds the 'km' unit
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
x=max(xint)+xmin+gpct(1, xlim), y=ymin,
vjust=0, hjust=0, label="km",
color="grey40", size=sizekm)
# For all graduations
for (x in xint) {
map = map +
# Draws the tick
annotate("segment",
x=x+xmin, xend=x+xmin, y=ymin, yend=ymax,
color="grey40", size=0.2) +
# Adds the value
annotate("text",
x=x+xmin, y=ymax+gpct(0.5, ylim),
vjust=0, hjust=0.5, label=x/1E3,
color="grey40", size=size)
}
map = map +
# Allows to crop shapefile without graphical problem
coord_sf(xlim=xlim, ylim=ylim,
expand=FALSE)
# If there is no margins specified
if (is.null(margin)) {
# Sets all margins to 0
map = map +
theme(plot.margin=margin(t=0, r=0, b=0, l=0,
unit="mm"))
# Otherwise
} else {
# Sets margins to the given ones
map = map +
theme(plot.margin=margin)
}
# Blank vector to store data about station
lon = c()
lat = c()
fill = c()
shape = c()
Value = c()
# For all code
for (k in 1:nCode) {
# Gets the code
code = Code[k]
if (j > 1) {
value = breakValue_code[j, i, k]
minValue = minBreakValue[j, i]
maxValue = maxBreakValue[j, i]
} else if (is.null(trend_period)) {
value = NA
minValue = NULL
maxValue = NULL
# 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
# Gets the risk of the test
alpha = list_df2plot[[i]]$alpha
# 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,]
# Extract start and end of trend periods
Start = df_trend_code$period_start[idPer_trend]
End = df_trend_code$period_end[idPer_trend]
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
# Extracts the corresponding data for the period
df_data_code_per =
df_data_code[df_data_code$Date >= Start
& df_data_code$Date <= End,]
# Same for trend
df_trend_code_per =
df_trend_code[df_trend_code$period_start == Start
& df_trend_code$period_end == End,]
# Computes the number of trend analysis selected
Ntrend = nrow(df_trend_code_per)
# If there is more than one trend on the same period
if (Ntrend > 1) {
# Takes only the first because they are similar
df_trend_code_per = df_trend_code_per[1,]
}
# If it is a flow variable
if (type == 'sévérité') {
# Computes the mean of the data on the period
dataMean = mean(df_data_code_per$Value,
na.rm=TRUE)
# Normalises the trend value by the mean
# of the data
value = df_trend_code_per$trend / dataMean
# If it is a date variable
} else if (type == 'saisonnalité') {
value = df_trend_code_per$trend
}
minValue = minTrendValue[idPer_trend, i]
maxValue = maxTrendValue[idPer_trend, i]
}
# Computes the color associated to the mean trend
color_res = get_color(value,
minValue,
maxValue,
palette_name='perso',
reverse=TRUE,
ncolor=256)
# Computes the colorbar info
palette_res = get_palette(minValue,
maxValue,
palette_name='perso',
reverse=TRUE,
ncolor=256,
nbTick=nbTick)
# The marker is a circle
shapetmp = 21
} else {
# If it is significative
# The computed color is stored
filltmp = color_res
# If the mean tend is positive
if (value >= 0) {
# Uses a triangle up for the shape
# of the marker
shapetmp = 24
# If negative
} else {
# Uses a triangle down for the shape
# of the marker
shapetmp = 25
}
} else if (pVal > alpha & colorForce) {
# The computed color is stored
filltmp = color_res
# The marker is a circle
shapetmp = 21
# If it is not significative
# The fill color is grey
filltmp = 'grey97'
# The marker is a circle
shapetmp = 21
# Extracts the localisation of the current station
lontmp =
df_meta[df_meta$code == code,]$L93X_m_BH
lattmp =
df_meta[df_meta$code == code,]$L93Y_m_BH
# Stores all the parameters
lon = c(lon, lontmp)
lat = c(lat, lattmp)
fill = c(fill, filltmp)
shape = c(shape, shapetmp)
Value = c(Value, value)
# If the trend analysis is significative a TRUE is stored
# Creates a tibble to stores all the data to plot
plot_map = tibble(lon=lon, lat=lat, fill=fill,
# If there is no specified station code to highlight
# (mini map)
if (is.null(codeLight)) {
plot_map_NOk = plot_map[!plot_map$OkVal,]
plot_map_Ok = plot_map[plot_map$OkVal,]
if (nrow(plot_map_NOk) > 0) {
map = map +
# Plots the point that are not
# significant first
geom_point(data=plot_map_NOk,
aes(x=lon, y=lat),
shape=shape[!OkVal],
size=5, stroke=1,
color='grey50', fill=fill[!OkVal])
}
if (nrow(plot_map_Ok) > 0) {
map = map +
# Plots the point that are significant last
geom_point(data=plot_map_Ok,
aes(x=lon, y=lat),
shape=shape[OkVal], size=5, stroke=1,
color='grey50', fill=fill[OkVal])
}
} else {
# Extract data of all stations not to highlight
plot_map_codeNo = plot_map[plot_map$code != codeLight,]
# Extract data of the station to highlight
plot_map_code = plot_map[plot_map$code == codeLight,]
# Plots only the localisation
map = map +
# For all stations not to highlight
geom_point(data=plot_map_codeNo,
aes(x=lon, y=lat),
shape=21, size=0.5, stroke=0.5,
color='grey50', fill='grey50') +
# For the station to highlight
geom_point(data=plot_map_code,
aes(x=lon, y=lat),
# Extracts the position of the tick of the colorbar
posTick = palette_res$posTick
# Extracts the label of the tick of the colorbar
labTick = palette_res$labTick
# Extracts the color corresponding to the tick of the colorbar
colTick = palette_res$colTick
# Spreading of the colorbar
valNorm = nbTick * 10
# Normalisation of the position of ticks
ytick = posTick / max(posTick) * valNorm
# If it is a flow variable
if (type == 'sévérité') {
# Formatting of label in pourcent
labTick = as.character(signif(labTick*100, 2))
# If it is a date variable
} else if (type == 'saisonnalité') {
# Formatting of label
labTick = as.character(signif(labTick, 2))
}
# X position of ticks all similar
xtick = rep(0, times=nbTick)
# Creates a tibble to store all parameters of colorbar
plot_palette = tibble(xtick=xtick, ytick=ytick,
colTick=colTick, labTick=labTick)
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
nbLine = as.integer(nchar(glose)/40) + 1
nbNewline = 0
nbLim = 43
gloseName = glose
nbChar = nchar(gloseName)
while (nbChar > nbLim) {
nbNewline = nbNewline + 1
posSpace = which(strsplit(gloseName, "")[[1]] == " ")
idNewline = which.min(abs(posSpace - nbLim * nbNewline))
posNewline = posSpace[idNewline]
gloseName = paste(substring(gloseName,
c(1, posNewline + 1),
c(posNewline,
nchar(gloseName))),
collapse="\n")
Newline = substr(gloseName,
posNewline + 2,
nchar(gloseName))
nbChar = nchar(Newline)
}
Yline = 0.6 + 0.47*nbNewline
Ytitle = Yline + 0.15
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
# New plot with void theme
title = ggplot() + theme_void() +
# Plots separation lines
geom_line(aes(x=c(-0.3, 3.9), y=c(0.05, 0.05)),
size=0.6, color="#00A3A8") +
geom_line(aes(x=c(-0.3, 3.9), y=c(Yline, Yline)),
size=0.6, color="#00A3A8") +
# Writes title
geom_shadowtext(data=tibble(x=-0.3, y=Ytitle,
label=var),
aes(x=x, y=y, label=label),
fontface="bold",
color="#00A3A8",
bg.colour="white",
hjust=0, vjust=0, size=10) +
# Writes title
geom_shadowtext(data=tibble(x=-0.3, y=0.2,
label=gloseName),
aes(x=x, y=y, label=label),
fontface="bold",
color="#00A3A8",
bg.colour="white",
hjust=0, vjust=0, size=3) +
# X axis
scale_x_continuous(limits=c(-0.3, 1 + 3),
expand=c(0, 0)) +
# Y axis
scale_y_continuous(limits=c(0, 10),
expand=c(0, 0)) +
# Margin
theme(plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm"))
# New plot with void theme
pal = ggplot() + theme_void() +
# Plots the point of the colorbar
geom_point(data=plot_palette,
aes(x=xtick, y=ytick),
shape=21, size=5, stroke=1,
color='white', fill=colTick)
periodName_trend = paste(
format(as.Date(trend_period[[idPer_trend]][1]),
'%Y'),
format(as.Date(trend_period[[idPer_trend]][2]),
'%Y'),
sep='-')
periodName1_mean = paste(
format(as.Date(mean_period[[1]][1]),
'%Y'),
format(as.Date(mean_period[[1]][2]),
'%Y'),
sep='-')
periodName2_mean = paste(
format(as.Date(mean_period[[2]][1]),
'%Y'),
format(as.Date(mean_period[[2]][2]),
'%Y'),
sep='-')
if (j > 1) {
ValueName1 = "Écarts observés entre"
ValueName2 = paste(periodName1_mean,
" et ",
periodName2_mean,
sep='')
# If it is a flow variable
if (type == 'sévérité') {
unit = bquote(bold("(%)"))
# If it is a date variable
} else if (type == 'saisonnalité') {
unit = bquote(bold("(jour)"))
}
} else {
ValueName1 = "Tendances observées"
ValueName2 = paste("sur la période ",
periodName_trend, sep='')
# If it is a flow variable
if (type == 'sévérité') {
unit = bquote(bold("(% par an)"))
# If it is a date variable
} else if (type == 'saisonnalité') {
unit = bquote(bold("(jour par an)"))
}
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
size=6, color='grey40') +
# Second line
annotate('text',
x=-0.3, y= valNorm + 26,
label=ValueName2,
hjust=0, vjust=0.5,
size=6, color='grey40') +
# Unit legend of the colorbar
annotate('text',
x=-0.3, y= valNorm + 14,
label=unit,
hjust=0, vjust=0.5,
size=4, color='grey40')
# For all the ticks
for (id in 1:nbTick) {
pal = pal +
# Adds the value
annotate('text', x=xtick[id]+0.3,
y=ytick[id],
label=bquote(bold(.(labTick[id]))),
hjust=0, vjust=0.7,
size=3, color='grey40')
}
if (j == 1) {
upLabel = bquote(bold("Hausse significative à 10%"))
noneLabel = bquote(bold("Non significatif à 10%"))
downLabel = bquote(bold("Baisse significative à 10%"))
yUp = -20
yNone = -29
yDown = -40
pal = pal +
# Up triangle in the marker legend
geom_point(aes(x=0, y=yUp),
shape=24, size=4, stroke=1,
color='grey50', fill='grey97') +
# Up triangle text legend
annotate('text',
x=0.3, y=yUp,
label=upLabel,
hjust=0, vjust=0.5,
size=3, color='grey40')
pal = pal +
# Circle in the marker legend
geom_point(aes(x=0, y=yNone),
shape=21, size=4, stroke=1,
color='grey50', fill='grey97') +
# Circle text legend
annotate('text',
x=0.3, y=yNone,
label=noneLabel,
hjust=0, vjust=0.7,
size=3, color='grey40')
pal = pal +
# Down triangle in the marker legend
geom_point(aes(x=0, y=yDown),
shape=25, size=4, stroke=1,
color='grey50', fill='grey97') +
# Down triangle text legend
annotate('text',
x=0.3, y=yDown,
label=downLabel,
hjust=0, vjust=0.5,
size=3, color='grey40')
}
# Normalises all the trend values for each station
# according to the colorbar
if (j > 1) {
yValue = (Value - minBreakValue[j, i]) / (maxBreakValue[j, i] - minBreakValue[j, i]) * valNorm
} else {
yValue = (Value - minTrendValue[idPer_trend, i]) / (maxTrendValue[idPer_trend, i] - minTrendValue[idPer_trend, i]) * valNorm
}
# Histogram distribution
# Computes the histogram of values
midsOk = res_hist$mids
# Histogram distribution
# Computes the histogram of values
res_hist = hist(yValueNOk, breaks=ytick, plot=FALSE)
# Extracts the number of counts per cells
countsNOk = res_hist$counts
counts = countsOk + countsNOk
# Blank vectors to store position of points of
# the distribution to plot
xValue = c()
yValue = c()
# Start X position of the distribution
start_hist = 1
# X separation bewteen point
hist_sep = 0.15
# Gets the maximun number of point of the distribution
maxCount = max(counts, na.rm=TRUE)
# Limit of the histogram
lim_hist = 2
# If the number of point will exceed the limit
if (maxCount * hist_sep > lim_hist) {
# Computes the right amount of space between points
hist_sep = lim_hist / maxCount
}
# For all cells of the histogram
# If the count in the current cell is not zero
if (counts[ii] != 0) {
# Stores the X positions of points of the
# distribution for the current cell
xValue = c(
xValue,
seq(start_hist,
start_hist+(counts[ii]-1)*hist_sep,
by=hist_sep))
}
# Stores the Y position which is the middle of the
# current cell the number of times it has been counted
yValue = c(yValue, rep(midsOk[ii],
times=counts[ii]))
color = c(color, rep('grey50',
times=countsOk[ii]))
color = c(color, rep('grey85',
times=countsNOk[ii]))
if (j > 1) {
shape = 21
} else {
if (midsOk[ii] > 0) {
shapetmp = 25
} else {
shapetmp = 24
}
shape = c(shape, rep(shapetmp,
times=countsOk[ii]))
shape = c(shape, rep(21,
times=countsNOk[ii]))
}
}
# Makes a tibble to plot the distribution
plot_value = tibble(xValue=xValue, yValue=yValue)
# Plots the point of the distribution
geom_point(data=plot_value,
aes(x=xValue, y=yValue),
alpha=1)
if (type == 'sévérité') {
labelArrow = 'Plus sévère'
} else if (type == 'saisonnalité') {
labelArrow = 'Plus tôt'
}
# Position of the arrow
xArrow = 3.3
# Arrow to show a worsening of the situation
geom_segment(aes(x=xArrow, y=valNorm*0.75,
xend=xArrow, yend=valNorm*0.25),
color='grey50', size=0.3,
arrow=arrow(length=unit(2, "mm"))) +
# Text associated to the arrow
x=xArrow+0.1, y=valNorm*0.5,
label=labelArrow,
angle=90,
hjust=0.5, vjust=1,
size=3, color='grey50')
pal = pal +
# X axis of the colorbar
scale_x_continuous(limits=c(-0.3, 4),
expand=c(0, 0)) +
# Y axis of the colorbar
scale_y_continuous(limits=c(-47, valNorm + 48),
expand=c(0, 0)) +
# Margin of the colorbar
theme(plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm"))
}
subsection = var
n_page = df_page$n[nrow(df_page)] + 1
df_page = bind_rows(
df_page,
tibble(section=section,
subsection=subsection,
n=n_page))
}
# If there is a foot note
if (foot_note) {
if (j > 1) {
footName = 'carte des écarts observés'
} else {
footName = 'carte des tendances observées'
}
if (is.null(df_page)) {
n_page = n_loop
}
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
AEAGlogo_file, INRAElogo_file,
FRlogo_file, foot_height)
# Stores the map, the title and the colorbar in a list
P = list(map, title, pal, foot)
LM = matrix(c(1, 1, 1, 2,
1, 1, 1, 3,
4, 4, 4, 4),
nrow=3, byrow=TRUE)
} else {
foot_height = 0
# Stores the map, the title and the colorbar in a list
P = list(map, title, pal)
LM = matrix(c(1, 1, 1, 2,
1, 1, 1, 3),
nrow=2, byrow=TRUE)
}
id_foot = 4
LMcol = ncol(LM)
LMrow = nrow(LM)
LM = rbind(rep(99, times=LMcol), LM, rep(99, times=LMcol))
LMrow = nrow(LM)
LM = cbind(rep(99, times=LMrow), LM, rep(99, times=LMrow))
LMcol = ncol(LM)
row_height = (height - 2*margin_size - foot_height) / (LMrow - 3)
Hcut = LM[, 2]
heightLM = rep(row_height, times=LMrow)
heightLM[Hcut == id_foot] = foot_height
Wcut = LM[(nrow(LM)-1),]
widthLM = rep(col_width, times=LMcol)
# Arranges the graphical object
plot = grid.arrange(grobs=P, layout_matrix=LM,
heights=heightLM, widths=widthLM)
# If there is no specified station code to highlight
# (mini map)
if (is.null(codeLight)) {
# Saving matrix plot
ggsave(plot=plot,
path=outdirTmp,
filename=paste(outname, '.pdf', sep=''),
width=width, height=height, units='cm', dpi=100)
}
# If there is no specified station code to highlight
# (mini map)
if (is.null(codeLight)) {
return (df_page)