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,
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)
# 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)
# 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),
" %)",
sep=''))
}
# If there is no specified station code to highlight
# (mini map)
if (is.null(codeLight)) {
# Sets the size of the countour
sizefr = 0.45
sizebs = 0.4
sizerv = 0.3
} else {
sizefr = 0.35
sizebs = 0.3
sizerv = 0.2
}
# Stores the coordonate system
cf = coord_fixed()
# Makes it the default one to remove useless warning
cf$default = TRUE
# Open a new plot with the personalise theme
map = ggplot() + theme_void() +
# theme(
# # 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
248
249
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
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]
pvalue = 0
} else {
# 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]
318
319
320
321
322
323
324
325
326
327
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
361
362
363
364
365
366
367
368
369
# 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]
pvalue = df_trend_code_per$p
}
# 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
if (pvalue <= alpha){
# 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
}
# 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])
}
# If there is a specified station code
} else {
# Extract data of all stations not to highlight
plot_map_codeNo = plot_map[plot_map$code != codeLight,]
# Extract data of the station to highlight
plot_map_code = plot_map[plot_map$code == codeLight,]
# Plots only the localisation
map = map +
# For all stations not to highlight
geom_point(data=plot_map_codeNo,
aes(x=lon, y=lat),
shape=21, size=0.5, stroke=0.5,
color='grey50', fill='grey50') +
# For the station to highlight
geom_point(data=plot_map_code,
aes(x=lon, y=lat),
}
# 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é') {
# 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)
nbLine = as.integer(nchar(glose)/40) + 1
nbNewline = 0
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
# New plot with void theme
title = ggplot() + theme_void() +
geom_line(aes(x=c(-0.3, 3.9), y=c(0.05, 0.05)),
geom_line(aes(x=c(-0.3, 3.9), y=c(Yline, Yline)),
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)
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
if (!is.null(trend_period)) {
periodName_trend = paste(
format(as.Date(trend_period[[idPer_trend]][1]),
'%Y'),
format(as.Date(trend_period[[idPer_trend]][2]),
'%Y'),
sep='-')
} else {
periodName_trend = NA
}
if (!is.null(mean_period)) {
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='-')
} else {
periodName1_mean = NA
periodName2_mean = NA
}
ValueName1 = "Écarts observés entre"
ValueName2 = paste(periodName1_mean,
" et ",
periodName2_mean,
sep='')
# If it is a flow variable
if (type == 'sévérité') {
# If it is a date variable
} else if (type == 'saisonnalité') {
ValueName1 = "Tendances observées"
ValueName2 = paste("sur la période ",
periodName_trend, sep='')
# If it is a flow variable
if (type == 'sévérité') {
# If it is a date variable
} else if (type == 'saisonnalité') {
}
}
pal = pal +
# Name of the colorbar
annotate('text',
x=-0.3, y= valNorm + 37,
label=ValueName1,
hjust=0, vjust=0.5,
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',
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')
upLabel = bquote(bold("Hausse significative à 10%"))
noneLabel = bquote(bold("Non significatif à 10%"))
downLabel = bquote(bold("Baisse significative à 10%"))
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
}
# Takes only the significative ones
res_hist = hist(yValue, breaks=ytick, plot=FALSE)
# Extracts the number of counts per cells
counts = res_hist$counts
# Extracts limits of cells
breaks = res_hist$breaks
# Extracts middle of cells
mids = res_hist$mids
# Blank vectors to store position of points of
# the distribution to plot
xValue = c()
yValue = c()
# Start X position of the distribution
# 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
for (ii in 1:length(mids)) {
# If the count in the current cell is not zero
if (counts[ii] != 0) {
# Stores the X positions of points of the
# distribution for the current cell
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(mids[ii], times=counts[ii]))
}
plot_value = tibble(xValue=xValue, yValue=yValue)
pal = pal +
geom_point(data=plot_value,
aes(x=xValue, y=yValue),
shape=21, color='white',
fill='grey50', stroke=0.4,
alpha=1)
if (type == 'sévérité') {
labelArrow = 'Plus sévère'
} else if (type == 'saisonnalité') {
labelArrow = 'Plus tôt'
}
pal = pal +
# 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
annotate('text',
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_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
}
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
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)