# \\\
# Copyright 2021-2022 Louis H�raut*1
#
# *1 INRAE, France
# louis.heraut@inrae.fr
#
# This file is part of ash R toolbox.
#
# ash R toolbox is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# ash R toolbox is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with ash R toolbox. If not, see <https://www.gnu.org/licenses/>.
# ///
#
#
# plotting/map.R
#
#
## 1. MAP PANEL
map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='', codeLight=NULL, margin=NULL, showSea=TRUE, verbose=TRUE) {
df_france = df_shapefile$france
df_bassin = df_shapefile$bassin
df_river = df_shapefile$river
# Number of type/variable
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=' / '))
}
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
if (df_trend_code_per$p <= p_threshold){
TrendMean_code[j, i, k] = trendMean
} else {
TrendMean_code[j, i, k] = NA
}
}
}
}
minTrendMean = apply(TrendMean_code, c(1, 2), min, na.rm=TRUE)
maxTrendMean = apply(TrendMean_code, c(1, 2), max, na.rm=TRUE)
ncolor = 256
nbTick = 10
for (i in 1:nbp) {
if (i > 1 & !is.null(codeLight)) {
break
}
type = list_df2plot[[i]]$type
outname = paste('map_', type, sep='')
if (verbose) {
print(paste('Map for variable :', type))
}
if (is.null(codeLight)) {
sizefr = 0.45
sizebs = 0.4
sizerv = 0.3
} else {
sizefr = 0.35
sizebs = 0.3
sizerv = 0.2
}
map = ggplot() + theme_void() +
# theme(plot.background=element_rect(fill=NA,
# color="#EC4899")) +
coord_fixed() +
geom_polygon(data=df_france,
aes(x=long, y=lat, group=group),
color=NA, fill="grey97")
if (!is.null(df_river)) {
map = map +
geom_path(data=df_river,
aes(x=long, y=lat, group=group),
color="grey85", size=sizerv)
}
map = map +
geom_polygon(data=df_bassin,
aes(x=long, y=lat, group=group),
color="grey70", fill=NA, size=sizebs) +
geom_polygon(data=df_france,
aes(x=long, y=lat, group=group),
color="grey40", fill=NA, size=sizefr)
if (showSea) {
xlim = c(280000, 790000)
ylim = c(6110000, 6600000)
} else {
xlim = c(305000, 790000)
ylim = c(6135000, 6600000)
}
if (is.null(codeLight)) {
xmin = gpct(7, xlim, shift=TRUE)
xint = c(0, 10*1E3, 50*1E3, 100*1E3)
ymin = gpct(5, ylim, shift=TRUE)
ymax = ymin + gpct(1, ylim)
size = 3
sizekm = 2.5
} else {
xmin = gpct(5, xlim, shift=TRUE)
xint = c(0, 100*1E3)
ymin = gpct(1, ylim, shift=TRUE)
ymax = ymin + gpct(3, ylim)
size = 2
sizekm = 1.5
}
map = map +
geom_line(aes(x=c(xmin, max(xint)+xmin),
y=c(ymin, ymin)),
color="grey40", size=0.2) +
annotate("text",
x=max(xint)+xmin+gpct(1, xlim), y=ymin,
vjust=0, hjust=0, label="km",
color="grey40", size=sizekm)
for (x in xint) {
map = map +
annotate("segment",
x=x+xmin, xend=x+xmin, y=ymin, yend=ymax,
color="grey40", size=0.2) +
annotate("text",
x=x+xmin, y=ymax+gpct(0.5, ylim),
vjust=0, hjust=0.5, label=x/1E3,
color="grey40", size=size)
}
map = map +
coord_sf(xlim=xlim, ylim=ylim,
expand=FALSE)
if (is.null(margin)) {
map = map +
theme(plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm"))
} else {
map = map +
theme(plot.margin=margin)
}
lon = c()
lat = c()
fill = c()
shape = c()
trend = c()
p_threshold_Ok = c()
for (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 = 'grey97'
shapetmp = 21
}
lontmp =
df_meta[df_meta$code == code,]$L93X_m_BH
lattmp =
df_meta[df_meta$code == code,]$L93Y_m_BH
lon = c(lon, lontmp)
lat = c(lat, lattmp)
fill = c(fill, filltmp)
shape = c(shape, shapetmp)
trend = c(trend, trendMean)
p_threshold_Ok = c(p_threshold_Ok,
df_trend_code_per$p <= p_threshold)
}
plot_map = tibble(lon=lon, lat=lat, fill=fill,
shape=shape, code=Code)
if (is.null(codeLight)) {
map = map +
geom_point(data=plot_map,
aes(x=lon, y=lat),
shape=shape, size=5, stroke=1,
color='grey50', fill=fill)
} else {
plot_map_codeNo = plot_map[plot_map$code != codeLight,]
plot_map_code = plot_map[plot_map$code == codeLight,]
map = map +
geom_point(data=plot_map_codeNo,
aes(x=lon, y=lat),
shape=21, size=0.5, stroke=0.5,
color='grey70', fill='grey70') +
geom_point(data=plot_map_code,
aes(x=lon, y=lat),
shape=21, size=1.5, stroke=0.5,
color='grey40', fill='grey40')
}
posTick = palette_res$posTick
labTick = palette_res$labTick
colTick = palette_res$colTick
nbTickmod = length(posTick)
valNorm = nbTickmod * 10
ytick = posTick / max(posTick) * valNorm
labTick = as.character(round(labTick*100, 2))
xtick = rep(0, times=nbTickmod)
plot_palette = tibble(xtick=xtick, ytick=ytick,
colTick=colTick, labTick=labTick)
title = ggplot() + theme_void() +
geom_line(aes(x=c(-0.3, 3.3), y=c(0.05, 0.05)),
size=0.6, color="#00A3A8") +
geom_shadowtext(data=tibble(x=-0.3, y=0.2,
label=type),
aes(x=x, y=y, label=label),
fontface="bold",
color="#00A3A8",
bg.colour="white",
hjust=0, vjust=0, size=10) +
scale_x_continuous(limits=c(-1, 1 + 3),
expand=c(0, 0)) +
scale_y_continuous(limits=c(0, 10),
expand=c(0, 0)) +
theme(plot.margin=margin(t=5, r=5, b=0, l=0, unit="mm"))
pal = ggplot() + theme_void() +
geom_point(data=plot_palette,
aes(x=xtick, y=ytick),
shape=21, size=5, stroke=1,
color='white', fill=colTick)
pal = pal +
annotate('text',
x=-0.3, y= valNorm + 23,
label="Tendance",
hjust=0, vjust=0.5,
size=6, color='grey40') +
annotate('text',
x=-0.2, y= valNorm + 13,
label=bquote(bold("% par an")),
hjust=0, vjust=0.5,
size=4, color='grey40')
for (j in 1:nbTickmod) {
pal = pal +
annotate('text', x=xtick[j]+0.3,
y=ytick[j],
label=bquote(bold(.(labTick[j]))),
hjust=0, vjust=0.7,
size=3, color='grey40')
}
pal = pal +
geom_point(aes(x=0, y=-20),
shape=24, size=4, stroke=1,
color='grey50', fill='grey97') +
annotate('text',
x=0.3, y=-20,
label=bquote(bold("Hausse significative � 10%")),
hjust=0, vjust=0.5,
size=3, color='grey40')
pal = pal +
geom_point(aes(x=0, y=-29),
shape=21, size=4, stroke=1,
color='grey50', fill='grey97') +
annotate('text',
x=0.3, y=-29,
label=bquote(bold("Non significatif � 10%")),
hjust=0, vjust=0.7,
size=3, color='grey40')
pal = pal +
geom_point(aes(x=0, y=-40),
shape=25, size=4, stroke=1,
color='grey50', fill='grey97') +
annotate('text',
x=0.3, y=-40,
label=bquote(bold("Baisse significative � 10%")),
hjust=0, vjust=0.5,
size=3, color='grey40')
yTrend = (trend - minTrendMean[idPer, i]) /
(maxTrendMean[idPer, i] - minTrendMean[idPer, i]) * valNorm
yTrend = yTrend[p_threshold_Ok]
## Random distribution ##
# xTrend = rnorm(length(yTrend), mean=1.75, sd=0.1)
## Histogram distribution ##
res_hist = hist(yTrend, breaks=ytick, plot=FALSE)
counts = res_hist$counts
breaks = res_hist$breaks
mids = res_hist$mids
xTrend = c()
yTrend = c()
start_hist = 1.25
hist_sep = 0.15
for (ii in 1:length(mids)) {
if (counts[ii] != 0) {
xTrend = c(xTrend,
seq(start_hist,
start_hist+(counts[ii]-1)*hist_sep,
by=hist_sep))
}
yTrend = c(yTrend, rep(mids[ii], times=counts[ii]))
}
## No touch distribution ##
# start_hist = 1.25
# min_xsep = 0.15
# min_ysep = 4
# xTrend = rep(start_hist, times=length(yTrend))
# for (ii in 1:length(yTrend)) {
# yTrendtmp = yTrend
# yTrendtmp[ii] = 1E99
# isinf_ysep = abs(yTrendtmp - yTrend[ii]) < min_ysep
# if (any(isinf_ysep) & !all(xTrend[which(isinf_ysep)] > start_hist)) {
# xTrend[ii] = max(xTrend[which(isinf_ysep)]) + min_xsep
# }
# }
plot_trend = tibble(xTrend=xTrend, yTrend=yTrend)
pal = pal +
geom_point(data=plot_trend,
aes(x=xTrend, y=yTrend),
# shape=21, size=1,
# color="grey20", fill="grey20")
alpha=0.4)
pal = pal +
geom_segment(aes(x=2.7, y=valNorm*0.75,
xend=2.7, yend=valNorm*0.25),
color='grey50', size=0.3,
arrow=arrow(length=unit(2, "mm"))) +
annotate('text',
x=2.8, y=valNorm*0.5,
label= "Plus s�v�re",
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)
}