An error occurred while loading the file. Please try again.
-
Dorchies David authored
Refs #111
dd7ef429
# \\\# 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## 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 variablemap_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='', codeLight=NULL, margin=NULL, showSea=TRUE, verbose=TRUE) { # Extract shapefiles df_france = df_shapefile$france df_bassin = df_shapefile$bassin df_river = df_shapefile$river # Number of variable/plot nbp = length(list_df2plot) # Get all different stations code Code = levels(factor(df_meta$code)) nCode = length(Code) # Gets a trend example df_trend = list_df2plot[[1]]$trend nPeriod_max = 0 for (code in 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 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)) # If the number of period for the trend is greater # than the current max period, stocks it if (nPeriod > nPeriod_max) { nPeriod_max = nPeriod } }7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
# Blank list to store time info by station code
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 all the code
for (j in 1:nCode) {
# Gets the code
code = Code[j]
# 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))
# Vector to store trend period
Periods = c()
# For all the period
for (i in 1:nPeriod_max) {
# Stocks period
Periods = append(Periods,
paste(substr(Start[i], 1, 4),
substr(End[i], 1, 4),
sep=' / '))
}
# Stores time info by station
Start_code[[j]] = Start
End_code[[j]] = End
Code_code[[j]] = code
Periods_code[[j]] = Periods
}
# Blank array to store mean of the trend for each
# station, perdiod and variable
TrendMean_code = array(rep(1, nPeriod_max*nbp*nCode),
dim=c(nPeriod_max, nbp, nCode))
# For all the period
for (j in 1:nPeriod_max) {
# For all the code
for (k in 1:nCode) {
# Gets the code
code = Code[k]
# For all variable
for (i in 1:nbp) {
# 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
p_threshold = list_df2plot[[i]]$p_threshold
# 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,]
# Gets the associated time info
Start = Start_code[Code_code == code][[1]][j]
End = End_code[Code_code == code][[1]][j]
Periods = Periods_code[Code_code == code][[1]][j]
# Extracts the corresponding data for the period
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
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,]
}
# Computes the mean of the data on the period
dataMean = mean(df_data_code_per$Qm3s, na.rm=TRUE)
# Normalises the trend value by the mean of the data
trendMean = df_trend_code_per$trend / dataMean
# If the p value is under the threshold
if (df_trend_code_per$p <= p_threshold){
# Stores the mean trend
TrendMean_code[j, i, k] = trendMean
# Otherwise
} else {
# Do not stocks it
TrendMean_code[j, i, k] = NA
}
}
}
}
# Compute the min and the max of the mean trend for all the station
minTrendMean = apply(TrendMean_code, c(1, 2), min, na.rm=TRUE)
maxTrendMean = apply(TrendMean_code, c(1, 2), max, na.rm=TRUE)
# Number of ticks for the colorbar
nbTick = 10
# 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 & !is.null(codeLight)) {
# Stop the for loop over the variable
break
}
# Extracts the variable of the plot
type = list_df2plot[[i]]$type
# Createsa name for the map
outname = paste('map_', type, sep='')
# If there is the verbose option
if (verbose) {
# Prints the name of the map
print(paste('Map for variable :', type))
}
# 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
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
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(plot.background=element_rect(fill=NA,
# color="#EC4899")) +
# 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)
}
map = map +
# Plot the hydrological basin
geom_polygon(data=df_bassin,
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 the sea needs to be shown
if (showSea) {
# Leaves space around the France
xlim = c(280000, 790000)
ylim = c(6110000, 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(7, 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(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
}
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
map = map +
# 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
annotate("text",
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()
trend = c()
p_threshold_Ok = c()
# For all code
for (code in Code) {
# 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
p_threshold = list_df2plot[[i]]$p_threshold
# 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,]
# Gets the associated time info
Start = Start_code[Code_code == code][[1]][idPer]
End = End_code[Code_code == code][[1]][idPer]
# 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
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
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,]
}
# Computes the mean of the data on the period
dataMean = mean(df_data_code_per$Qm3s, na.rm=TRUE)
# Normalises the trend value by the mean of the data
trendMean = df_trend_code_per$trend / dataMean
# Computes the color associated to the mean trend
color_res = get_color(trendMean,
minTrendMean[idPer, i],
maxTrendMean[idPer, i],
palette_name='perso',
reverse=TRUE,
ncolor=256)
# Computes the colorbar info
palette_res = get_palette(minTrendMean[idPer, i],
maxTrendMean[idPer, i],
palette_name='perso',
reverse=TRUE,
ncolor=256,
nbTick=nbTick)
# If it is significative
if (df_trend_code_per$p <= p_threshold){
# The computed color is stored
filltmp = color_res
# If the mean tend is positive
if (trendMean >= 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
} else {
# 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)
trend = c(trend, trendMean)
# If the trend analysis is significative a TRUE is stored
p_threshold_Ok = c(p_threshold_Ok,
df_trend_code_per$p <= p_threshold)
}
# Creates a tibble to stores all the data to plot
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490
plot_map = tibble(lon=lon, lat=lat, fill=fill,
shape=shape, code=Code)
# If there is no specified station code to highlight (mini map)
if (is.null(codeLight)) {
map = map +
# Plots the trend point
geom_point(data=plot_map,
aes(x=lon, y=lat),
shape=shape, size=5, stroke=1,
color='grey50', fill=fill)
# 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='grey70', fill='grey70') +
# For the station to highlight
geom_point(data=plot_map_code,
aes(x=lon, y=lat),
shape=21, size=1.5, stroke=0.5,
color='grey40', fill='grey40')
}
# 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
# Formatting of label in pourcent
labTick = as.character(round(labTick*100, 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)
# New plot with void theme
title = ggplot() + theme_void() +
# Plots separation line
geom_line(aes(x=c(-0.3, 3.3), y=c(0.05, 0.05)),
size=0.6, color="#00A3A8") +
# Writes title
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) +
# X axis
scale_x_continuous(limits=c(-1, 1 + 3),
expand=c(0, 0)) +
# Y axis
491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560
scale_y_continuous(limits=c(0, 10),
expand=c(0, 0)) +
# Margin
theme(plot.margin=margin(t=5, r=5, 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)
pal = pal +
# Name of the colorbar
annotate('text',
x=-0.3, y= valNorm + 23,
label="Tendance",
hjust=0, vjust=0.5,
size=6, color='grey40') +
# Unit legend of the colorbar
annotate('text',
x=-0.2, y= valNorm + 13,
label=bquote(bold("% par an")),
hjust=0, vjust=0.5,
size=4, color='grey40')
# For all the ticks
for (j in 1:nbTick) {
pal = pal +
# Adds the value
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 +
# Up triangle in the marker legend
geom_point(aes(x=0, y=-20),
shape=24, size=4, stroke=1,
color='grey50', fill='grey97') +
# Up triangle text legend
annotate('text',
x=0.3, y=-20,
label=bquote(bold("Hausse significative à 10%")),
hjust=0, vjust=0.5,
size=3, color='grey40')
pal = pal +
# Circle in the marker legend
geom_point(aes(x=0, y=-29),
shape=21, size=4, stroke=1,
color='grey50', fill='grey97') +
# Circle text legend
annotate('text',
x=0.3, y=-29,
label=bquote(bold("Non significatif à 10%")),
hjust=0, vjust=0.7,
size=3, color='grey40')
pal = pal +
# Down triangle in the marker legend
geom_point(aes(x=0, y=-40),
shape=25, size=4, stroke=1,
color='grey50', fill='grey97') +
# Down triangle text legend
annotate('text',
x=0.3, y=-40,
label=bquote(bold("Baisse significative à 10%")),
561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630
hjust=0, vjust=0.5,
size=3, color='grey40')
# Normalises all the trend values for each station
# according to the colorbar
yTrend = (trend - minTrendMean[idPer, i]) /
(maxTrendMean[idPer, i] - minTrendMean[idPer, i]) * valNorm
# Takes only the significative ones
yTrend = yTrend[p_threshold_Ok]
## Random distribution ##
# xTrend = rnorm(length(yTrend), mean=1.75, sd=0.1)
## Histogram distribution ##
# Computes the histogram of the trend
res_hist = hist(yTrend, 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
xTrend = c()
yTrend = c()
# Start X position of the distribution
start_hist = 1.25
# X separation bewteen point
hist_sep = 0.15
# 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
xTrend = c(xTrend,
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
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
# }
# }
# Makes a tibble to plot the trend distribution
plot_trend = tibble(xTrend=xTrend, yTrend=yTrend)
pal = pal +