diff --git a/plotting/layout.R b/plotting/layout.R index e0c659253eb640c0340fc6ec3e253229620780ef..ec53ea93d97ba5f8ac6167fb4a6e59bc2c002885 100644 --- a/plotting/layout.R +++ b/plotting/layout.R @@ -143,7 +143,7 @@ datasheet_layout = function (df_data, df_meta, layout_matrix, missRect=FALSE, time_header=NULL, info_header=TRUE, foot_note=FALSE, info_ratio=1, time_ratio=2, - var_ratio=3, foot_height=0.5, + var_ratio=3, foot_height=1, df_shapefile=NULL, resources_path=NULL, AEAGlogo_file=NULL, @@ -532,7 +532,7 @@ summary_panel = function (df_page, foot_note, foot_height, resources_path, AEAGl n_page = df_page$n[df_page$section == sec_name & df_page$subsection == subsec_name][1] - line = paste(idS, ".", idSS, ". ", + line = paste("<b>", idS, ".", idSS, ".</b> ", subsec_name, "<br>", sep='') page = paste("p.", n_page, "<br>", sep='') @@ -696,7 +696,7 @@ foot_panel = function (name, n_page, resources_path, AEAGlogo_file, INRAElogo_fi gp=gpar(col="#00A3A8", fontsize=8)) gtext_date = richtext_grob(text_date, - x=1, y=0.4, + x=1, y=0.55, margin=unit(c(t=0, r=0, b=0, l=0), "mm"), hjust=1, vjust=0.5, gp=gpar(col="#00A3A8", fontsize=6)) @@ -707,13 +707,15 @@ foot_panel = function (name, n_page, resources_path, AEAGlogo_file, INRAElogo_fi AEAGlogo_img = readPNG(AEAGlogo_path) AEAGlogo_grob = rasterGrob(AEAGlogo_img, + y=0.49, + vjust=0.5, width=unit(0.7*foot_height, "cm")) INRAElogo_img = readPNG(INRAElogo_path) INRAElogo_grob = rasterGrob(INRAElogo_img, - y=0.57, + y=0.565, vjust=0.5, - width=unit(1.1*foot_height, "cm")) + width=unit(1.08*foot_height, "cm")) FRlogo_img = readPNG(FRlogo_path) FRlogo_grob = rasterGrob(FRlogo_img, @@ -722,7 +724,7 @@ foot_panel = function (name, n_page, resources_path, AEAGlogo_file, INRAElogo_fi P = list(void, FRlogo_grob, INRAElogo_grob, AEAGlogo_grob, - gtext_page, gtext_date) + gtext_page, gtext_date) # Creates the matrix layout LM = matrix(c(1, 2, 3, 4, 5, @@ -732,7 +734,7 @@ foot_panel = function (name, n_page, resources_path, AEAGlogo_file, INRAElogo_fi # And sets the relative width of each plot widths = rep(1, times=ncol(LM)) - widths[2] = 0.18 + widths[2] = 0.2 widths[3] = 0.25 widths[4] = 0.2 diff --git a/plotting/map.R b/plotting/map.R index d86006ad4c3b09575f021b5fca9b9764aaa8e2a1..5f6ede5dfda8aa7536e6d8852a7cba87431aae36 100644 --- a/plotting/map.R +++ b/plotting/map.R @@ -200,10 +200,6 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1, minTrendValue = apply(TrendValue_code, c(1, 2), min, na.rm=TRUE) maxTrendValue = apply(TrendValue_code, c(1, 2), max, na.rm=TRUE) - - - - # If there is a 'mean_period' if (!is.null(mean_period)) { # Blank vectors to store info about breaking analysis @@ -362,8 +358,15 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1, # Open a new plot with the personalise theme map = ggplot() + theme_void() + - # theme(plot.background=element_rect(fill=NA, - # color="#EC4899")) + + + # 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 @@ -393,6 +396,29 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1, 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, + 6360000) + 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 @@ -478,7 +504,7 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1, fill = c() shape = c() Value = c() - alpha_Ok = c() + OkVal = c() # For all code for (k in 1:nCode) { # Gets the code @@ -560,27 +586,35 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1, ncolor=256, nbTick=nbTick) - # If it is significative - if (pvalue <= alpha){ + if (j > 1) { # 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 + # 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 } else { - # Uses a triangle down for the shape - # of the marker - shapetmp = 25 + # The fill color is grey + filltmp = 'grey97' + # The marker is a circle + shapetmp = 21 } - # 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 @@ -596,23 +630,39 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1, shape = c(shape, shapetmp) Value = c(Value, value) # If the trend analysis is significative a TRUE is stored - alpha_Ok = c(alpha_Ok, - pvalue <= alpha) + OkVal = c(OkVal, pvalue <= alpha) } # Creates a tibble to stores all the data to plot plot_map = tibble(lon=lon, lat=lat, fill=fill, - shape=shape, code=Code) + shape=shape, code=Code, OkVal=OkVal) # 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 + + 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,] @@ -629,7 +679,7 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1, # For the station to highlight geom_point(data=plot_map_code, aes(x=lon, y=lat), - shape=21, size=1.5, stroke=0.25, + shape=21, size=2, stroke=0.5, color='grey97', fill='#00A3A8') } @@ -736,34 +786,27 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1, size=3, color='grey40') } - yUp = -20 - yNone = -29 - - if (j > 1) { - upLabel = bquote(bold("Hausse")) - noneLabel = NULL - downLabel = bquote(bold("Baisse")) - yDown = -29 - } else { + 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 + + # 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') - if (!is.null(noneLabel)) { pal = pal + # Circle in the marker legend geom_point(aes(x=0, y=yNone), @@ -775,20 +818,21 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1, 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') + } - 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) { @@ -798,7 +842,7 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1, } # Takes only the significative ones - yValue = yValue[alpha_Ok] + yValue = yValue[OkVal] # Histogram distribution # Computes the histogram of values diff --git a/processing/analyse.R b/processing/analyse.R index 6ee5525f3be62205b29a3c63b268581f0fb8eba6..267b891d953ac72eb36d511ca33d7028a95623ca 100644 --- a/processing/analyse.R +++ b/processing/analyse.R @@ -292,7 +292,7 @@ get_VCN10trend = function (df_data, df_meta, period, alpha, sampleSpan, yearLac_ return (res_VCN10trend) } -### 1.4. DEB date +### 1.4. tDEB date which_underfirst = function (L, UpLim, select_longest=TRUE) { ID = which(L <= UpLim) @@ -328,7 +328,7 @@ which_underfirst = function (L, UpLim, select_longest=TRUE) { return (id) } -get_DEBtrend = function (df_data, df_meta, period, alpha, sampleSpan, yearLac_day, thresold_type='VCN10', select_longest=TRUE) { +get_tDEBtrend = function (df_data, df_meta, period, alpha, sampleSpan, yearLac_day, thresold_type='VCN10', select_longest=TRUE) { # Get all different stations code Code = levels(factor(df_meta$code)) @@ -375,7 +375,7 @@ get_DEBtrend = function (df_data, df_meta, period, alpha, sampleSpan, yearLac_da # Set the max interval period as the minimal possible Imax = 0 # Blank tibble for data to return - df_DEBtrendB = tibble() + df_tDEBtrendB = tibble() # For all periods for (per in period) { @@ -408,8 +408,8 @@ get_DEBtrend = function (df_data, df_meta, period, alpha, sampleSpan, yearLac_da # print(df_QT) - df_DEBEx = tibble() - df_DEBlist = list(data=tibble(), info=tibble()) + df_tDEBEx = tibble() + df_tDEBlist = list(data=tibble(), info=tibble()) # For all the code for (k in 1:nCode) { @@ -422,13 +422,13 @@ get_DEBtrend = function (df_data, df_meta, period, alpha, sampleSpan, yearLac_da df_data_roll_code = df_data_roll[df_data_roll$code == code,] # Prepare the data to fit the entry of extract.Var - df_DEBlist_code = prepare(df_data_roll_code, + df_tDEBlist_code = prepare(df_data_roll_code, colnamegroup=c('code')) QT_code = df_QT$Thresold[df_QT$code == code] # Compute the yearly min over the averaged data - df_DEBEx_code = extract.Var(data.station=df_DEBlist_code, + df_tDEBEx_code = extract.Var(data.station=df_tDEBlist_code, funct=which_underfirst, period=per, timestep='year', @@ -436,26 +436,26 @@ get_DEBtrend = function (df_data, df_meta, period, alpha, sampleSpan, yearLac_da UpLim=QT_code, select_longest=select_longest) - df_DEBEx_code$group1 = k - df_DEBlist_code$data$group = k - df_DEBlist_code$info$group = k + df_tDEBEx_code$group1 = k + df_tDEBlist_code$data$group = k + df_tDEBlist_code$info$group = k - # Converts index of the DEB to the julian date associated - df_DEBEx_code = prepare_date(df_DEBEx_code, - df_DEBlist_code) + # Converts index of the tDEB to the julian date associated + df_tDEBEx_code = prepare_date(df_tDEBEx_code, + df_tDEBlist_code) # Store the results - df_DEBEx = bind_rows(df_DEBEx, df_DEBEx_code) + df_tDEBEx = bind_rows(df_tDEBEx, df_tDEBEx_code) - df_DEBlist$data = bind_rows(df_DEBlist$data, - df_DEBlist_code$data) + df_tDEBlist$data = bind_rows(df_tDEBlist$data, + df_tDEBlist_code$data) - df_DEBlist$info = bind_rows(df_DEBlist$info, - df_DEBlist_code$info) + df_tDEBlist$info = bind_rows(df_tDEBlist$info, + df_tDEBlist_code$info) } # Compute the trend analysis - df_DEBtrend = Estimate.stats(data.extract=df_DEBEx, + df_tDEBtrend = Estimate.stats(data.extract=df_tDEBEx, level=alpha) # Get the associated time interval @@ -464,25 +464,25 @@ get_DEBtrend = function (df_data, df_meta, period, alpha, sampleSpan, yearLac_da if (I > Imax) { # Store it and the associated data and info Imax = I - df_DEBlistB = df_DEBlist - df_DEBExB = df_DEBEx + df_tDEBlistB = df_tDEBlist + df_tDEBExB = df_tDEBEx } # Specify the period of analyse - df_DEBtrend = get_period(per, df_DEBtrend, df_DEBEx, - df_DEBlist) + df_tDEBtrend = get_period(per, df_tDEBtrend, df_tDEBEx, + df_tDEBlist) # Store the trend - df_DEBtrendB = bind_rows(df_DEBtrendB, df_DEBtrend) + df_tDEBtrendB = bind_rows(df_tDEBtrendB, df_tDEBtrend) } # Clean results of trend analyse - res_DEBtrend = clean(df_DEBtrendB, df_DEBExB, df_DEBlistB) - return (res_DEBtrend) + res_tDEBtrend = clean(df_tDEBtrendB, df_tDEBExB, df_tDEBlistB) + return (res_tDEBtrend) } -### 1.5. CEN date +### 1.5. tCEN date # Realises the trend analysis of the date of the minimum 10 day # average flow over the year (VCN10) hydrological variable -get_CENtrend = function (df_data, df_meta, period, alpha, sampleSpan, yearLac_day) { +get_tCENtrend = function (df_data, df_meta, period, alpha, sampleSpan, yearLac_day) { # Get all different stations code Code = levels(factor(df_meta$code)) @@ -516,24 +516,24 @@ get_CENtrend = function (df_data, df_meta, period, alpha, sampleSpan, yearLac_da # Set the max interval period as the minimal possible Imax = 0 # Blank tibble for data to return - df_CENtrendB = tibble() + df_tCENtrendB = tibble() # For all periods for (per in period) { # Prepare the data to fit the entry of extract.Var - df_CENlist = prepare(df_data_roll, colnamegroup=c('code')) + df_tCENlist = prepare(df_data_roll, colnamegroup=c('code')) # Compute the yearly min over the averaged data - df_CENEx = extract.Var(data.station=df_CENlist, + df_tCENEx = extract.Var(data.station=df_tCENlist, funct=which.min, period=per, timestep='year', pos.datetime=1) - # Converts index of the CEN to the julian date associated - df_CENEx = prepare_date(df_CENEx, df_CENlist) + # Converts index of the tCEN to the julian date associated + df_tCENEx = prepare_date(df_tCENEx, df_tCENlist) # Compute the trend analysis - df_CENtrend = Estimate.stats(data.extract=df_CENEx, + df_tCENtrend = Estimate.stats(data.extract=df_tCENEx, level=alpha) # Get the associated time interval @@ -542,19 +542,19 @@ get_CENtrend = function (df_data, df_meta, period, alpha, sampleSpan, yearLac_da if (I > Imax) { # Store it and the associated data and info Imax = I - df_CENlistB = df_CENlist - df_CENExB = df_CENEx + df_tCENlistB = df_tCENlist + df_tCENExB = df_tCENEx } # Specify the period of analyse - df_CENtrend = get_period(per, df_CENtrend, df_CENEx, - df_CENlist) + df_tCENtrend = get_period(per, df_tCENtrend, df_tCENEx, + df_tCENlist) # Store the trend - df_CENtrendB = bind_rows(df_CENtrendB, df_CENtrend) + df_tCENtrendB = bind_rows(df_tCENtrendB, df_tCENtrend) } # Clean results of trend analyse - res_CENtrend = clean(df_CENtrendB, df_CENExB, df_CENlistB) - return (res_CENtrend) + res_tCENtrend = clean(df_tCENtrendB, df_tCENExB, df_tCENlistB) + return (res_tCENtrend) } diff --git a/processing/extract.R b/processing/extract.R index fc4f61ffedaacbf9df1ac8b5235a765b15ce75a0..aac50b48f9d6d5deaf167919361989a6781e5b15 100644 --- a/processing/extract.R +++ b/processing/extract.R @@ -477,7 +477,7 @@ extract_data = function (computer_data_path, filedir, filename, ## 4. SHAPEFILE MANAGEMENT # Generates a list of shapefiles to draw a hydrological map of # the France -ini_shapefile = function (computer_data_path, fr_shpdir, fr_shpname, bs_shpdir, bs_shpname, sbs_shpdir, sbs_shpname, rv_shpdir, rv_shpname, riv=TRUE) { +ini_shapefile = function (computer_data_path, fr_shpdir, fr_shpname, bs_shpdir, bs_shpname, sbs_shpdir, sbs_shpname, rv_shpdir, rv_shpname, is_river=TRUE) { # Path for shapefile fr_shppath = file.path(computer_data_path, fr_shpdir, fr_shpname) @@ -501,7 +501,7 @@ ini_shapefile = function (computer_data_path, fr_shpdir, fr_shpname, bs_shpdir, df_subbassin = tibble(fortify(subbassin)) # If the river shapefile needs to be load - if (riv) { + if (is_river) { # Hydrographic network river = readOGR(dsn=rv_shppath, verbose=FALSE) ### trop long ### river = river[which(river$Classe == 1),] diff --git a/script.R b/script.R index a8ce827b0026981f29f7c456b3aff311f3d4b7de..b049e48badc65443583ac29ef1fb7b91c57f939b 100644 --- a/script.R +++ b/script.R @@ -114,6 +114,9 @@ sampleSpan = c('05-01', '11-30') ## MAP +# Is the hydrological network needs to be plot +is_river = TRUE + # Path to the shapefile for france contour from 'computer_data_path' fr_shpdir = 'map/france' fr_shpname = 'gadm36_FRA_0.shp' @@ -261,17 +264,17 @@ res_VCN10trend = get_VCN10trend(df_data, df_meta, yearLac_day=yearLac_day) # Start date for low water trend -res_DEBtrend = get_DEBtrend(df_data, df_meta, +res_tDEBtrend = get_tDEBtrend(df_data, df_meta, period=trend_period, alpha=alpha, sampleSpan=sampleSpan, thresold_type='VCN10', select_longest=TRUE, yearLac_day=yearLac_day) -# res_DEBtrend = read_listofdf(resdir, 'res_DEBtrend') +# res_tDEBtrend = read_listofdf(resdir, 'res_tDEBtrend') # Center date for low water trend -res_CENtrend = get_CENtrend(df_data, df_meta, +res_tCENtrend = get_tCENtrend(df_data, df_meta, period=trend_period, alpha=alpha, sampleSpan=sampleSpan, @@ -295,7 +298,7 @@ df_shapefile = ini_shapefile(computer_data_path, fr_shpdir, fr_shpname, bs_shpdir, bs_shpname, sbs_shpdir, sbs_shpname, - rv_shpdir, rv_shpname, riv=TRUE) + rv_shpdir, rv_shpname, is_river=is_river) ### 4.1. Simple time panel to criticize station data # Plot time panel of debit by stations @@ -322,29 +325,33 @@ datasheet_layout(toplot=c( res_QAtrend$data, res_QMNAtrend$data, res_VCN10trend$data, - res_DEBtrend$data, - res_CENtrend$data), + res_tDEBtrend$data, + res_tCENtrend$data + ), df_trend=list( res_QAtrend$trend, res_QMNAtrend$trend, res_VCN10trend$trend, - res_DEBtrend$trend, - res_CENtrend$trend), + res_tDEBtrend$trend, + res_tCENtrend$trend + ), var=list( 'QA', 'QMNA', 'VCN10', - 'DEB', - 'CEN'), + 'tDEB', + 'tCEN' + ), type=list( 'sévérité', 'sévérité', 'sévérité', 'saisonnalité', - 'saisonnalité'), + 'saisonnalité' + ), layout_matrix=matrix(c(1, 2, 3, 4, 5), ncol=1), @@ -357,7 +364,7 @@ datasheet_layout(toplot=c( info_ratio=2, time_ratio=2, var_ratio=3, - foot_height=1, + foot_height=1.25, df_shapefile=df_shapefile, figdir=figdir, filename_opt='',