From 38eef77e377230c25e587dd42b55b2b2aa9328ab Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 4 Jun 2021 19:35:28 +0200 Subject: [PATCH 01/11] UseCase2 and add Observed Terciles plot wind --- inst/doc/UseCase1_WindEvent_March2018.R | 38 +++- ...e2_PrecipitationDownscaling_RainFARM_RF4.R | 172 ++++++++++++++++++ ...aunch_UseCase2_PrecipitationDownscaling.sh | 10 + 3 files changed, 216 insertions(+), 4 deletions(-) create mode 100644 inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R create mode 100644 inst/doc/launch_UseCase2_PrecipitationDownscaling.sh diff --git a/inst/doc/UseCase1_WindEvent_March2018.R b/inst/doc/UseCase1_WindEvent_March2018.R index 2baa597b..657d1d06 100644 --- a/inst/doc/UseCase1_WindEvent_March2018.R +++ b/inst/doc/UseCase1_WindEvent_March2018.R @@ -63,11 +63,12 @@ for (mm in 1:3) { grid = 'r360x181') str(wind_ref$Dates) dim(wind_ref$data) + print(wind_ref$Dates$start) wind_ref_terciles <- rbind(wind_ref_terciles, quantile(MeanDims(wind_ref$data, c('lat', 'lon')), c(0.3, 0.6))) wind_ref_extremes <- rbind(wind_ref_extremes, quantile(MeanDims(wind_ref$data, c('lat', 'lon')), c(0.1, 0.9))) - #source("/esarchive/scratch/nperez/git/cstools/R/CST_BiasCorrection.R") + source("/esarchive/scratch/nperez/git/cstools/R/CST_BiasCorrection.R") library(multiApply) wind_fsct <- CST_BiasCorrection(exp = wind_hcst, obs = wind_ref, @@ -125,7 +126,8 @@ for (mm in 1:3) { width = 1000, height = 1000, units = 'px', res = 144) PlotMostLikelyQuantileMap(probs = Mean_PB, lon = wind_fsct$lon, lat = wind_fsct$lat, intylat = 2, intxlon = 2, - coast_width = 1.5, legend_scale = 0.8, cat_dim = 'bin', + coast_width = 1.5, legend_scale = 0.8, + cat_dim = 'bin', dot_size = 2, dots = filtered_obs_terciles[,,1,1,1,1], toptitle = c(paste0("Initialized on ", month.name[as.numeric(months_in_advance[mm])]))) @@ -136,7 +138,7 @@ visual <- data.frame(dec = as.vector(MeanDims(wind_fsct_BC[[3]]$data, c('lat', ' jan = as.vector(MeanDims(wind_fsct_BC[[2]]$data, c('lat', 'lon'))), feb = as.vector(MeanDims(wind_fsct_BC[[1]]$data, c('lat', 'lon')))) - wind_obs <- CST_Load(var = 'windagl100', obs = list(obs_path), + wind_obs_areave <- CST_Load(var = 'windagl100', obs = list(obs_path), sdates = '20180301', nmember = 1, leadtimemin = 1, leadtimemax = 1, storefreq = "monthly", sampleperiod = 1, @@ -152,8 +154,36 @@ agg_png("/esarchive/scratch/nperez/CSTools_manuscript/Wind/PlotForecast_IP.png", PlotForecastPDF(visual, tercile.limits = wind_ref_terciles, extreme.limits = wind_ref_extremes, var.name = "Wind Speed 100 m (m/s)", title = "Bias Corrected forecasts at IP", - fcst.names = c("December", "January", "February"), obs = as.vector(wind_obs$data)) + fcst.names = c("December", "January", "February"), obs = as.vector(wind_obs_areave$data)) dev.off() +# Plotting observed terciles: + +wind_obs_obstercile <- Apply(list(wind_obs$data), target_dims = NULL, + fun = function(x, tercile) { + if (x <= tercile[1]) { + res <- 1 + } else if (x > tercile[2]) { + res <- 3 + } else { + res <- 2 + } + return(res) + }, tercile = wind_ref_terciles[1,])$output1 + +agg_png("/esarchive/scratch/nperez/CSTools_manuscript/Wind/MostLikely_Observed_obstercile.png", + width = 1000, height = 1000, units = 'px', res = 144) + +s2dv::PlotEquiMap(wind_obs_obstercile, + lon = wind_obs$lon, lat = wind_obs$lat, + brks = c(0,1,2,3), + cols = c("#6BAED6FF", "#FFEDA0FF", "#FC4E2AFF"), + intylat = 2, intxlon = 2, + coast_width = 1.5, filled.continents = FALSE, + toptitle = c("Observed terciles March 2018")) +dev.off() + + + print("DONE") diff --git a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R new file mode 100644 index 00000000..19f59fa5 --- /dev/null +++ b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R @@ -0,0 +1,172 @@ +#!/usr/bin/env Rscript +rm(list=ls()); gc(); + +### Creation date: 3rd June 2021 +# Author: N. Pérez-Zanón +# Refer CSTools package and manuscript when using it. +# ---------------------------------------- +# This code is divided in 4 steps plus visualization of the results +# STEP 1: Compute Slope values for RainFARM downscaling method +# STEP 2: Apply Quantile Mapping Correction to Seasonal Precipitation Forecast +# STEP 3: Compute Weights for RainFARM downscaling method +# STEP 4: Apply RainFARM downscaling method +# ---------------------------------------- +# Note: STEP 3 requires a machine with internet connexion. +# In this file, the lines are commented since they have been run and the +# result saved on disk, then the result is loaded. +# ---------------------------------------- +# +# Load required libraries and setup output directory: +library(CSTools) +library(ClimProjDiags) +library(zeallot) +library(ragg) +dir_output <- '/esarchive/scratch/nperez/CSTools_manuscript/v20210603/' #slash end + +# -------------------------------------------- +# STEP 1: +# -------------------------------------------- +era5 <- list(name = 'era5', + path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc') +years <- unlist(lapply(1993:2018, function(x){ + paste0(x, sprintf("%02d",1:12), '01')})) +era5 <- CST_Load(var = 'prlr', + exp = list(era5), + sdates = years, nmember = 1, + storefreq = "daily", sampleperiod = 1, + latmin = 37.5, latmax = 53.25, lonmin = 2.5, lonmax = 18.25, + output = 'lonlat', nprocs = 1) +era5$data <- era5$data * 24 * 3600 * 1000 # with or without this line -> same result +era5 <- CST_SplitDim(era5, split_dim = 'sdate', indices = rep(1:12, 26)) +slope <- CST_RFSlope(era5, time_dim = c('sdate', 'ftime'), kmin = 5) +save(slope, file = paste0(dir_output, 'Slope.RDS'), version = 2) + + +# -------------------------------------------- +# STEP 2: +# -------------------------------------------- +StartDates <- paste0(1993:2018, '1101') +exp <- list(name = 'ecmwfS5', + path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc") +obs <- list(name = 'era5', + path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc') +#obs <- list(name = 'era5', path = '/esarchive/scratch/nperez/ERA5/daily_total_prlr_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc') + +c(exp, obs) %<-% CST_Load(var = 'prlr', exp = list(exp), obs = list(obs), + sdates = StartDates, nmember = 25, + storefreq = "daily", sampleperiod = 1, + latmin = 42, latmax = 49, lonmin = 4, lonmax = 11, + output = 'lonlat', nprocs = 1) +exp <- CST_SplitDim(exp, split_dim = c('ftime')) +obs <- CST_SplitDim(obs, split_dim = c('ftime')) +exp$data <- exp$data * 24 * 3600 * 1000 +obs$data <- obs$data * 24 * 3600 * 1000 +exp$data[which(exp$data < 0)] <- 0 +exp.qm <- CST_QuantileMapping(exp, obs, method = "QUANT", + wet.day = FALSE, + sample_dims = c('member', 'sdate', 'ftime'), + ncores = 4) +save(exp.qm, file = paste0(dir_output, 'ExpQM.RDS', version = 2) + +# -------------------------------------------- +# STEP 3: +# -------------------------------------------- +#library(raster); library(s2dv); library(CSTools) +#worldclim <- getData("worldclim", var = "prec", res = 0.5, lon = 5, lat = 45) +#wc_month <- lapply(1:12, FUN = function(x) { +# res <- crop(worldclim[[x]], +# extent(3.5, 11.5, 41.5, 49.5)) +# res <- as.array(res) +# names(dim(res)) <- c('lat', 'lon', 'month') +# return(res) +# }) +#xy <- xyFromCell(crop(worldclim[[1]], extent(3.5, 11.5, 41.5, 49.5)), +# 1:length(crop(worldclim[[1]], extent(3.5, 11.5, 41.5, 49.5)))) +#lons <- unique(xy[,1]) +#lats <- unique(xy[,2]) +#wc_month <- unlist(wc_month) +#dim(wc_month) <- c(lat = length(lats), lon = length(lons), monthly = 12) +#wc_month <- Reorder(wc_month, c('lon', 'lat', 'monthly')) +#wc_month <- s2dv_cube(data = wc_month, lon = lons, lat = lats, +# Datasets = 'WorldClim') +#weight <- CST_RFWeights(wc_month, lon = 4:11, lat = 49:42, nf = 4) +#save(weight, file = paste0(dir_output, 'weights.RDS')) +load(paste0(dir_output, 'weights.RDS')) + +# -------------------------------------------- +# STEP 4: +# -------------------------------------------- + +weights <- Subset(weight$data, along = 'monthly', indices = c(11,12,1:1:6)) +slope <- Subset(slope, along = 'monthly', indices = c(11,12,1:1:6), drop = 'non-selected') +fs <- CST_RainFARM(exp.qm, nf = 4, + weights = weights, slope = slope, + kmin = 1, nens = 10, verbose = TRUE, + time_dim = c("member", "sdate", "ftime"), nprocs = 8, + drop_realization = TRUE) +newfs <- CST_MergeDims(fs, merge_dims = c("ftime", "monthly"), + na.rm = TRUE) + +newfs$Dates[[1]] <- exp$Dates[[1]] +CST_SaveExp(newfs, destination = output_dir) + +# -------------------------------------------- +# Visualization +# -------------------------------------------- +agg_png(paste0(output_dir, "EXP_11dec.png"), + width = 500, height = 600, units = 'px',res = 144) +PlotEquiMap(exp$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, + filled.continents = FALSE, bar_limits = c(0,40), + intylat = 2, intxlon = 2, title_scale = 0.7, + toptitle = 'ECMWF-S5C3S', units = 'precipitation (mm)') +dev.off() +agg_png(paste0(output_dir, "EXPQM_11dec.png"), + width = 500, height = 600, units = 'px',res = 144) +PlotEquiMap(exp.qm$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, + filled.continents = FALSE, bar_limits = c(0,40), + intylat = 2, intxlon = 2, title_scale = 0.7, + toptitle = 'Bias Corrected', units = 'precipitation (mm)') +dev.off() +agg_png(paste0(output_dir, "Down_11dec.png"), + width = 500, height = 600, units = 'px',res = 144) +PlotEquiMap(fs$data[1,1,1,11,,,2],lon = fs$lon, lat = fs$lat, + filled.continents = FALSE, bar_limits = c(0,40), + intylat = 2, intxlon = 2, title_scale = 0.7, + toptitle = 'Downsacaled', units = 'precipitation (mm)') +dev.off() +agg_png(paste0(output_dir, "WeightsNov.png"), + width = 500, height = 600, units = 'px',res = 144) +PlotEquiMap(weight$data[,,11], lon = weight$lon, lat = weight$lat, + filled.continents = FALSE, + intylat = 2, intxlon = 2, + toptitle = 'December Weights WorldClim2') +dev.off() +agg_png(paste0(output_dir, "Slope.png"), + width = 500, height = 600, units = 'px',res = 144) +plot(1:12, slope[3:12,1:2] type = 'b', main = 'slope', pch = 3) +dev.off() + +# Plot ForecastPDF +obsts <- MeanDims(obs$data, c('lat', 'lon'), na.rm = T) +print(quantile(obsts, c(0.1,0.3,0.6,0.9), na.rm = T)) +expts <- MeanDims(exp$data, c('lat', 'lon'), na.rm = T) +exp.qmts <- MeanDims(exp.qm$data, c('lat', 'lon'), na.rm = T) +empty <- array(NA, c(dataset = 2, member = 225, sdate = 26, ftime = 31, monthly = 8)) +data <- abind(expts, exp.qmts, along = 1) +names(dim(data)) <- names(dim(expts)) +data <- abind(data, empty, along = 2) +names(dim(data)) <- names(dim(expts)) +fsts <- MeanDims(fs, c('lat', 'lon'), na.rm = T) +data <- abind(data, fsts, along = 1) +names(dim(data)) <- c('dataset', 'members', 'sdate', 'ftime', 'monthly') +agg_png("/esarchive/scratch/nperez/CSTools_manuscript/v20201201/FiguresPDF_11DecemberAll2.png", + width = 750, height = 650, units = 'px', res = 144) +PlotForecastPDF(data[,,1,11,2], tercile.limits = c(0.58, 2.4), obs = obsts[,,1,11,2], + extreme.limits = c(0.06, 7.26), color.set = 'hydro', add.ensmemb = 'no', + var.name = "Precipitation (mm)", title = "Forecasts issued on Nov 1993 for 11th December 1993", + fcst.names = c("ECMWFS5C3S", "Bias Corrected", "Downscaled")) +dev.off() + + + + diff --git a/inst/doc/launch_UseCase2_PrecipitationDownscaling.sh b/inst/doc/launch_UseCase2_PrecipitationDownscaling.sh new file mode 100644 index 00000000..55c3e848 --- /dev/null +++ b/inst/doc/launch_UseCase2_PrecipitationDownscaling.sh @@ -0,0 +1,10 @@ +#!/bin/bash +#BSUB -W 2:00 +#BSUB -n 16 +#BSUB -M 7000 +#BSUB -J RainFARM_Downsc +#BSUB -o /esarchive/scratch/nperez/CSTools_manuscript/v20210603/lsf-%J.out +#BSUB -e /esarchive/scratch/nperez/CSTools_manuscript/v20210603/lsf-%J.err + +module load R +Rscript /esarchive/scratch/nperez/git/cstools/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R -- GitLab From 42fb98ad109305d50b53e82309f84aad223fac53 Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 8 Jun 2021 18:03:22 +0200 Subject: [PATCH 02/11] fixing plots --- inst/doc/UseCase1_WindEvent_March2018.R | 36 ++++++++++++---- ...e2_PrecipitationDownscaling_RainFARM_RF4.R | 41 +++++++++++-------- 2 files changed, 50 insertions(+), 27 deletions(-) diff --git a/inst/doc/UseCase1_WindEvent_March2018.R b/inst/doc/UseCase1_WindEvent_March2018.R index 657d1d06..719f5bf3 100644 --- a/inst/doc/UseCase1_WindEvent_March2018.R +++ b/inst/doc/UseCase1_WindEvent_March2018.R @@ -1,7 +1,22 @@ +#!/usr/bin/env Rscript +rm(list=ls()); gc(); + +### Creation date: 3rd June 2021 +# Author: N. Pérez-Zanón +# Refer CSTools package and manuscript when using it. +# ---------------------------------------- +# Wind speed bias adjustment for assessment of an extreme event +# The System5-ECMWF downloaded from C3S seasonal forecasts initialized +# in December 2017, January 2018 and February 2018 +# This code includes the bias adjustent and the results visualization +# ---------------------------------------- library(CSTools) library(s2dv) library(ragg) +library(multiApply) +output_dir <- "/esarchive/scratch/nperez/CSTools_manuscript/Wind/" + exp_path <- list(name = "ECMWFS5", path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_f6h/$VAR_NAME$_$START_DATE$.nc") @@ -15,6 +30,7 @@ months_in_advance <- c('02', '01', '12') wind_fsct_BC <- list() wind_ref_terciles <- NULL wind_ref_extremes <- NULL +wind_thres_latlon <- NULL # Observations March 2018 wind_obs <- CSTools::CST_Load(var = 'windagl100', obs = list(obs_path), sdates = '20180301', nmember = 1, @@ -64,12 +80,12 @@ for (mm in 1:3) { str(wind_ref$Dates) dim(wind_ref$data) print(wind_ref$Dates$start) + wind_ref_terciles <- rbind(wind_ref_terciles, quantile(MeanDims(wind_ref$data, c('lat', 'lon')), c(0.3, 0.6))) wind_ref_extremes <- rbind(wind_ref_extremes, quantile(MeanDims(wind_ref$data, c('lat', 'lon')), c(0.1, 0.9))) source("/esarchive/scratch/nperez/git/cstools/R/CST_BiasCorrection.R") - library(multiApply) wind_fsct <- CST_BiasCorrection(exp = wind_hcst, obs = wind_ref, exp_cor = wind_fcst) @@ -119,10 +135,10 @@ for (mm in 1:3) { } return(y) })$output1 - + wind_thres_latlon <- abind::abind(wind_thres_latlon, thres, along = 4) source("/esarchive/scratch/nperez/git/cstools/R/PlotCombinedMap.R") source("/esarchive/scratch/nperez/git/cstools/R/PlotMostLikelyQuantileMap.R") - agg_png(paste0("/esarchive/scratch/nperez/CSTools_manuscript/Wind/MostLikely_", mm, "_obstercile.png"), + agg_png(paste0(output_dir, "MostLikely_", mm, "_obstercile.png"), width = 1000, height = 1000, units = 'px', res = 144) PlotMostLikelyQuantileMap(probs = Mean_PB, lon = wind_fsct$lon, lat = wind_fsct$lat, intylat = 2, intxlon = 2, @@ -149,7 +165,7 @@ visual <- data.frame(dec = as.vector(MeanDims(wind_fsct_BC[[3]]$data, c('lat', ' print("IS DATA LOADED") print("Wait") -agg_png("/esarchive/scratch/nperez/CSTools_manuscript/Wind/PlotForecast_IP.png", +agg_png(paste0(output_dir, "PlotForecast_IP.png"), width = 1000, height = 1000, units = 'px',res = 144) PlotForecastPDF(visual, tercile.limits = wind_ref_terciles, extreme.limits = wind_ref_extremes, @@ -158,8 +174,10 @@ PlotForecastPDF(visual, tercile.limits = wind_ref_terciles, dev.off() # Plotting observed terciles: - -wind_obs_obstercile <- Apply(list(wind_obs$data), target_dims = NULL, +names(dim(wind_thres_latlon)) <- c('thres', 'lat', 'lon', 'sdate') +wind_thres_latlon <- ClimProjDiags::Subset(wind_thres_latlon, indices = 1, along = 'sdate') +wind_obs_obstercile <- Apply(list(wind_obs$data, wind_thres_latlon), + target_dims = list(NULL, 'thres'), fun = function(x, tercile) { if (x <= tercile[1]) { res <- 1 @@ -169,9 +187,9 @@ wind_obs_obstercile <- Apply(list(wind_obs$data), target_dims = NULL, res <- 2 } return(res) - }, tercile = wind_ref_terciles[1,])$output1 + })$output1 -agg_png("/esarchive/scratch/nperez/CSTools_manuscript/Wind/MostLikely_Observed_obstercile.png", +agg_png(paste0(output_dir, "MostLikely_Observed_obstercile.png"), width = 1000, height = 1000, units = 'px', res = 144) s2dv::PlotEquiMap(wind_obs_obstercile, @@ -183,7 +201,7 @@ s2dv::PlotEquiMap(wind_obs_obstercile, toptitle = c("Observed terciles March 2018")) dev.off() - +# All gridpoints are above normal observed tercile. print("DONE") diff --git a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R index 19f59fa5..6917ae6d 100644 --- a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R +++ b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R @@ -3,6 +3,8 @@ rm(list=ls()); gc(); ### Creation date: 3rd June 2021 # Author: N. Pérez-Zanón +# Adapted from the original version Terzago et al. 2020 +# https://drive.google.com/file/d/1qp2gbtKdBl4XmsyOeaEhFENwpeUuJwkf/view # Refer CSTools package and manuscript when using it. # ---------------------------------------- # This code is divided in 4 steps plus visualization of the results @@ -40,7 +42,7 @@ era5$data <- era5$data * 24 * 3600 * 1000 # with or without this line -> same re era5 <- CST_SplitDim(era5, split_dim = 'sdate', indices = rep(1:12, 26)) slope <- CST_RFSlope(era5, time_dim = c('sdate', 'ftime'), kmin = 5) save(slope, file = paste0(dir_output, 'Slope.RDS'), version = 2) - +slope_plot <- slope # -------------------------------------------- # STEP 2: @@ -66,7 +68,7 @@ exp.qm <- CST_QuantileMapping(exp, obs, method = "QUANT", wet.day = FALSE, sample_dims = c('member', 'sdate', 'ftime'), ncores = 4) -save(exp.qm, file = paste0(dir_output, 'ExpQM.RDS', version = 2) +save(exp.qm, file = paste0(dir_output, 'ExpQM.RDS'), version = 2) # -------------------------------------------- # STEP 3: @@ -91,6 +93,7 @@ save(exp.qm, file = paste0(dir_output, 'ExpQM.RDS', version = 2) # Datasets = 'WorldClim') #weight <- CST_RFWeights(wc_month, lon = 4:11, lat = 49:42, nf = 4) #save(weight, file = paste0(dir_output, 'weights.RDS')) + load(paste0(dir_output, 'weights.RDS')) # -------------------------------------------- @@ -108,42 +111,44 @@ newfs <- CST_MergeDims(fs, merge_dims = c("ftime", "monthly"), na.rm = TRUE) newfs$Dates[[1]] <- exp$Dates[[1]] -CST_SaveExp(newfs, destination = output_dir) +CST_SaveExp(newfs, destination = dir_output) # -------------------------------------------- # Visualization # -------------------------------------------- -agg_png(paste0(output_dir, "EXP_11dec.png"), - width = 500, height = 600, units = 'px',res = 144) +library(s2dv) +agg_png(paste0(dir_output, "EXP_11dec.png"), + width = 1000, height = 1000, units = 'px',res = 144) PlotEquiMap(exp$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, filled.continents = FALSE, bar_limits = c(0,40), intylat = 2, intxlon = 2, title_scale = 0.7, toptitle = 'ECMWF-S5C3S', units = 'precipitation (mm)') dev.off() -agg_png(paste0(output_dir, "EXPQM_11dec.png"), - width = 500, height = 600, units = 'px',res = 144) +agg_png(paste0(dir_output, "EXPQM_11dec.png"), + width = 1000, height = 1000, units = 'px',res = 144) PlotEquiMap(exp.qm$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, filled.continents = FALSE, bar_limits = c(0,40), intylat = 2, intxlon = 2, title_scale = 0.7, - toptitle = 'Bias Corrected', units = 'precipitation (mm)') + toptitle = 'Bias Adjusted', units = 'precipitation (mm)') dev.off() -agg_png(paste0(output_dir, "Down_11dec.png"), - width = 500, height = 600, units = 'px',res = 144) +agg_png(paste0(dir_output, "Down_11dec.png"), + width = 1000, height = 1000, units = 'px',res = 144) PlotEquiMap(fs$data[1,1,1,11,,,2],lon = fs$lon, lat = fs$lat, filled.continents = FALSE, bar_limits = c(0,40), intylat = 2, intxlon = 2, title_scale = 0.7, toptitle = 'Downsacaled', units = 'precipitation (mm)') dev.off() -agg_png(paste0(output_dir, "WeightsNov.png"), - width = 500, height = 600, units = 'px',res = 144) -PlotEquiMap(weight$data[,,11], lon = weight$lon, lat = weight$lat, +agg_png(paste0(dir_output, "WeightsDec.png"), + width = 1000, height = 1000, units = 'px',res = 144) +PlotEquiMap(weight$data[,,12], lon = weight$lon, lat = weight$lat, filled.continents = FALSE, intylat = 2, intxlon = 2, toptitle = 'December Weights WorldClim2') dev.off() -agg_png(paste0(output_dir, "Slope.png"), - width = 500, height = 600, units = 'px',res = 144) -plot(1:12, slope[3:12,1:2] type = 'b', main = 'slope', pch = 3) +agg_png(paste0(dir_output, "Slope.png"), + width = 1000, height = 1000, units = 'px',res = 144) +plot(1:12, slope_plot, type = 'b', main = 'slope', pch = 3) +line(12, slope_plot[12], type = 'p', col = 'red', pch = 3) dev.off() # Plot ForecastPDF @@ -159,12 +164,12 @@ names(dim(data)) <- names(dim(expts)) fsts <- MeanDims(fs, c('lat', 'lon'), na.rm = T) data <- abind(data, fsts, along = 1) names(dim(data)) <- c('dataset', 'members', 'sdate', 'ftime', 'monthly') -agg_png("/esarchive/scratch/nperez/CSTools_manuscript/v20201201/FiguresPDF_11DecemberAll2.png", +agg_png(paste0(dir_output, "/FiguresPDF_11DecemberAll2.png", width = 750, height = 650, units = 'px', res = 144) PlotForecastPDF(data[,,1,11,2], tercile.limits = c(0.58, 2.4), obs = obsts[,,1,11,2], extreme.limits = c(0.06, 7.26), color.set = 'hydro', add.ensmemb = 'no', var.name = "Precipitation (mm)", title = "Forecasts issued on Nov 1993 for 11th December 1993", - fcst.names = c("ECMWFS5C3S", "Bias Corrected", "Downscaled")) + fcst.names = c("ECMWFS5C3S", "Bias Adjusted", "Downscaled")) dev.off() -- GitLab From e0ba23a0b2d1bd91cf20e8a33962c7a264dfba94 Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 11 Jun 2021 11:20:18 +0200 Subject: [PATCH 03/11] Use cases and extra_string parameter in SaveExp --- R/CST_SaveExp.R | 28 ++- ..._PrecipitationDownscaling_RainFARM_RF100.R | 195 ++++++++++++++++++ ...e2_PrecipitationDownscaling_RainFARM_RF4.R | 63 +++--- man/CST_SaveExp.Rd | 4 +- man/SaveExp.Rd | 7 +- 5 files changed, 259 insertions(+), 38 deletions(-) create mode 100644 inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R diff --git a/R/CST_SaveExp.R b/R/CST_SaveExp.R index 9c689ff7..f9290b18 100644 --- a/R/CST_SaveExp.R +++ b/R/CST_SaveExp.R @@ -13,6 +13,7 @@ #'folder tree: destination/experiment/variable/. By default the function #'creates and saves the data into the folder "CST_Data" in the working #'directory. +#'@param extra_string a character string to be include as part of the file name, for instance, to identify member or realization. It would be added to the file name between underscore characters. #' #'@seealso \code{\link{CST_Load}}, \code{\link{as.s2dv_cube}} and \code{\link{s2dv_cube}} #' @@ -29,7 +30,7 @@ #'} #' #'@export -CST_SaveExp <- function(data, destination = "./CST_Data") { +CST_SaveExp <- function(data, destination = "./CST_Data", extra_string = NULL) { if (!is.character(destination) & length(destination) > 1) { stop("Parameter 'destination' must be a character string of one element ", "indicating the name of the file (including the folder if needed) ", @@ -55,7 +56,8 @@ CST_SaveExp <- function(data, destination = "./CST_Data") { SaveExp(data = data$data, lon = data$lon, lat = data$lat, Dataset = names(data$Datasets), var_name = var_name, units = units, cdo_grid_name = cdo_grid_name, projection = projection, - startdates = sdates, Dates = time_values, destination) + startdates = sdates, Dates = time_values, destination, + extra_string = extra_string) } #'Save an experiment in a format compatible with CST_Load #'@description This function is created for compatibility with CST_Load/Load for saving post-processed datasets such as those calibrated of downscaled with CSTools functions @@ -73,8 +75,9 @@ CST_SaveExp <- function(data, destination = "./CST_Data") { #'@param cdo_grid_name a character string indicating the name of the grid e.g.: 'r360x181' #'@param projection a character string indicating the projection name #'@param destination a character string indicating the path where to store the NetCDF files +#'@param extra_string a character string to be include as part of the file name, for instance, to identify member or realization. #' -#'@return the function creates as many files as sdates per dataset. Each file could contain multiple members +#'@return the function creates as many files as sdates per dataset. Each file could contain multiple members. It would be added to the file name between underscore characters. #' The path will be created with the name of the variable and each Datasets. #' #'@import ncdf4 @@ -102,7 +105,8 @@ CST_SaveExp <- function(data, destination = "./CST_Data") { #'} #'@export SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, - cdo_grid_name, projection, destination) { + cdo_grid_name, projection, destination, + extra_string = NULL) { dimname <- names(dim(data)) if (any(dimname == "ftime")) { dimname[which(dimname == "ftime")] <- "time" @@ -136,7 +140,11 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, stop("Element 'data' in parameter 'data' has more than one 'sdate'", " dimension.") } - + if (!is.null(extra_string)) { + if (!is.character(extra_string)) { + stop("Parameter 'extra_string' must be a character string.") + } + } dataset_pos <- which(dimname == "dataset" | dimname == "dat") dims <- dim(data) if (length(dataset_pos) == 0) { @@ -227,7 +235,7 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, NULL, 'time'), fun = .saveExp, var_name = var_name, units = units, dims_var = dims_var, cdo_grid_name = cdo_grid_name, projection = projection, - destination = path) + destination = path, extra_string = extra_string) } } @@ -252,7 +260,7 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, #Dates <- as.Date(c("1900-11-01", "1900-12-01", "1901-01-01", "1901-02-01", "1901-03-01")) #.saveExp(data, sdate, Dates, var_name, units, dims_var, cdo_grid_name = 'r360x181', projection = 'none', destination) .saveExp <- function(data, sdate, Dates, var_name, units, dims_var, - cdo_grid_name, projection, destination) { + cdo_grid_name, projection, destination, extra_string) { dim_names <- names(dim(data)) if (any(dim_names != c('longitude', 'latitude', 'member', 'time'))) { data <- Reorder(data, c('longitude', 'latitude', 'member', 'time')) @@ -266,7 +274,11 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, datanc <- ncvar_def(name = var_name, units = units, dim = dims_var, missval = -99999) - file_name <- paste0(var_name, "_", sdate, ".nc") + if (is.null(extra_string)) { + file_name <- paste0(var_name, "_", sdate, ".nc") + } else { + file_name <- paste0(var_name, "_", extra_string, "_", sdate, ".nc") + } full_filename <- file.path(destination, file_name) file_nc <- nc_create(full_filename, datanc) ncvar_put(file_nc, datanc, data) diff --git a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R new file mode 100644 index 00000000..7b057180 --- /dev/null +++ b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R @@ -0,0 +1,195 @@ +#!/usr/bin/env Rscript +rm(list=ls()); gc(); + +### Creation date: 9th June 2021 +# Author: N. Pérez-Zanón +# Adapted from the original version Terzago et al. 2020 +# https://drive.google.com/file/d/1qp2gbtKdBl4XmsyOeaEhFENwpeUuJwkf/view +# Refer CSTools package and manuscript when using it. +# ---------------------------------------- +# This code is using divided in 4 steps plus visualization of the results +# Taking advantage of +# STEP 1: Compute Slope values for RainFARM downscaling method +# STEP 2: Apply Quantile Mapping Correction to Seasonal Precipitation Forecast +# STEP 3: Compute Weights for RainFARM downscaling method +# STEP 4: Apply RainFARM downscaling method +# ---------------------------------------- +# Note: STEP 3 requires a machine with internet connexion. +# In this file, the lines are commented since they have been run and the +# result saved on disk, then the result is loaded. +# ---------------------------------------- +# +# Load required libraries and setup output directory: +library(CSTools) +library(ClimProjDiags) +library(zeallot) +library(ragg) +dir_output <- '/esarchive/scratch/nperez/CSTools_manuscript/v20210603/' #slash end + +# -------------------------------------------- +# STEP 1: +# -------------------------------------------- +era5 <- list(name = 'era5', + path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc') +years <- unlist(lapply(1993:2018, function(x){ + paste0(x, sprintf("%02d",1:12), '01')})) +era5 <- CST_Load(var = 'prlr', + exp = list(era5), + sdates = years, nmember = 1, + storefreq = "daily", sampleperiod = 1, + latmin = 37.5, latmax = 53.25, lonmin = 2.5, lonmax = 18.25, + output = 'lonlat', nprocs = 1) +era5$data <- era5$data * 24 * 3600 * 1000 # with or without this line -> same result +era5 <- CST_SplitDim(era5, split_dim = 'sdate', indices = rep(1:12, 26)) +slope <- CST_RFSlope(era5, time_dim = c('sdate', 'ftime'), kmin = 5) +save(slope, file = paste0(dir_output, 'Slope.RDS'), version = 2) +slope_plot <- slope + +# -------------------------------------------- +# STEP 2: +# -------------------------------------------- +StartDates <- paste0(1993:2018, '1101') +exp <- list(name = 'ecmwfS5', + path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc") +obs <- list(name = 'era5', + path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc') +#obs <- list(name = 'era5', path = '/esarchive/scratch/nperez/ERA5/daily_total_prlr_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc') + +c(exp, obs) %<-% CST_Load(var = 'prlr', exp = list(exp), obs = list(obs), + sdates = StartDates, nmember = 25, + storefreq = "daily", sampleperiod = 1, + latmin = 42, latmax = 49, lonmin = 4, lonmax = 11, + output = 'lonlat', nprocs = 1) +exp <- CST_SplitDim(exp, split_dim = c('ftime')) +obs <- CST_SplitDim(obs, split_dim = c('ftime')) +exp$data <- exp$data * 24 * 3600 * 1000 +obs$data <- obs$data * 24 * 3600 * 1000 +exp$data[which(exp$data < 0)] <- 0 +exp.qm <- CST_QuantileMapping(exp, obs, method = "QUANT", + wet.day = FALSE, + sample_dims = c('member', 'sdate', 'ftime'), + ncores = 4) +save(exp.qm, file = paste0(dir_output, 'ExpQM.RDS'), version = 2) + +# -------------------------------------------- +# STEP 3: +# -------------------------------------------- +#library(raster); library(s2dv); library(CSTools) +#worldclim <- getData("worldclim", var = "prec", res = 0.5, lon = 5, lat = 45) +#wc_month <- lapply(1:12, FUN = function(x) { +# res <- crop(worldclim[[x]], +# extent(3.5, 11.5, 41.5, 49.5)) +# res <- as.array(res) +# names(dim(res)) <- c('lat', 'lon', 'month') +# return(res) +# }) +#xy <- xyFromCell(crop(worldclim[[1]], extent(3.5, 11.5, 41.5, 49.5)), +# 1:length(crop(worldclim[[1]], extent(3.5, 11.5, 41.5, 49.5)))) +#lons <- unique(xy[,1]) +#lats <- unique(xy[,2]) +#wc_month <- unlist(wc_month) +#dim(wc_month) <- c(lat = length(lats), lon = length(lons), monthly = 12) +#wc_month <- Reorder(wc_month, c('lon', 'lat', 'monthly')) +#wc_month <- s2dv_cube(data = wc_month, lon = lons, lat = lats, +# Datasets = 'WorldClim') +#weight <- CST_RFWeights(wc_month, lon = 4:11, lat = 49:42, nf = 100) +#save(weight, file = paste0(dir_output, 'weightsRF100.RDS')) + +load(paste0(dir_output, 'weightsRF100.RDS')) + +# -------------------------------------------- +# Visualization +# -------------------------------------------- +agg_png(paste0(dir_output, "RF100_WeightsDec.png"), + width = 1000, height = 1100, units = 'px',res = 144) +PlotEquiMap(weight$data[,,12], lon = weight$lon, lat = weight$lat, + filled.continents = FALSE, title_scale = 1, + intylat = 2, intxlon = 2, + toptitle = 'December Weights RF 100') +dev.off() + +# -------------------------------------------- +# STEP 4: +# -------------------------------------------- + +weights <- Subset(weight$data, along = 'monthly', indices = c(11, 12, 1:6)) +slope <- Subset(slope, along = 'monthly', indices = c(11, 12, 1:6), + drop = 'non-selected') +k = 1 # To create the member ID +#----- +# To be removed when CSTools 4.0.1 is published: +source("/gpfs/scratch/bsc32/bsc32339/CSTools_manuscript/CST_SaveExp.R") +library(multiApply) +library(ncdf4) +library(s2dv) +#----- +for (realizations in 1:10) { + for (member in 1:25) { + result <- data # to store the data + result$data <- NULL + for (month in 1:8) { + data <- exp.qm # to take the correct piece of data + data$data <- data$data[1, member, , , , , month] + fs <- CST_RainFARM(data, nf = 100, + weights = weights, slope = slope[month], + kmin = 1, nens = 1, verbose = TRUE, + nprocs = 8, + drop_realization = TRUE) + print(dim(fs$data)) + result$data <- abind::abind(result$data, fs$data, along = 5) + if (month == 2 & member == 1 & realization == 1) { + # ---------------------------- + # Visualization: + # ---------------------------- + agg_png(paste0(dir_output, "RF100_Down_11dec.png"), + width = 1000, height = 1100, units = 'px',res = 144) + PlotEquiMap(fs$data[1,11,,],lon = fs$lon, lat = fs$lat, + filled.continents = FALSE, bar_limits = c(0,40), + intylat = 2, intxlon = 2, title_scale = 1, + triangle_ends = c(TRUE, FALSE), + toptitle = 'Downsacaled RF 100', units = 'precipitation (mm)') + dev.off() + # PlotForecastPDF + library(abind) + obsts <- MeanDims(obs$data, c('lat', 'lon'), na.rm = T) + print(quantile(obsts, c(0.1, 0.3, 0.6, 0.9), na.rm = T)) + expts <- MeanDims(exp$data, c('lat', 'lon'), na.rm = T) + exp.qmts <- MeanDims(exp.qm$data, c('lat', 'lon'), na.rm = T) + empty <- array(NA, c(dataset = 2, member = 225, sdate = 26, + ftime = 31, monthly = 8)) + data <- abind(expts, exp.qmts, along = 1) + names(dim(data)) <- names(dim(expts)) + data <- abind(data, empty, along = 2) + names(dim(data)) <- names(dim(expts)) + fsts <- MeanDims(fs$data, c('lat', 'lon'), na.rm = T) + print(dim(fsts)) + print(dim(data)) + data <- abind(data, fsts, along = 1) + names(dim(data)) <- c('dataset', 'members', 'sdate', 'ftime', 'monthly') + agg_png(paste0(dir_output, "/FiguresPDF_RF100_11DecemberAll2.png"), + width = 750, height = 650, units = 'px', res = 144) + PlotForecastPDF(data[,,1,11,2], tercile.limits = c(0.58, 2.4), + obs = obsts[,,1,11,2], + extreme.limits = c(0.06, 7.26), + color.set = 'hydro', add.ensmemb = 'no', + var.name = "Precipitation (mm)", + title = "Forecasts issued on Nov 1993 for 11th December 1993", + fcst.names = c("ECMWFS5C3S", "Bias Adjusted", "Downscaled")) + dev.off() + } + result$lon <- fs$lon + result$lat <- fs$lat + result <- CST_MergeDims(result, merge_dims = c("ftime", "monthly"), + na.rm = TRUE) + result$Dataset <- paste0('RF100_ECMWFC3S_QM_member_', member, '_real_', + realizations) + result$Dates[[1]] <- exp$Dates[[1]] + CST_SaveExp(result, destination = dir_output, + extra_string = paste0('member', k)) + gc() + k = k + 1 + rm(list= list('fs', 'result', 'data')) + } +} + + diff --git a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R index 6917ae6d..44b1a94e 100644 --- a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R +++ b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R @@ -92,16 +92,17 @@ save(exp.qm, file = paste0(dir_output, 'ExpQM.RDS'), version = 2) #wc_month <- s2dv_cube(data = wc_month, lon = lons, lat = lats, # Datasets = 'WorldClim') #weight <- CST_RFWeights(wc_month, lon = 4:11, lat = 49:42, nf = 4) -#save(weight, file = paste0(dir_output, 'weights.RDS')) +#save(weight, file = paste0(dir_output, 'weightsRF4.RDS')) -load(paste0(dir_output, 'weights.RDS')) +load(paste0(dir_output, 'weightsRF4.RDS')) # -------------------------------------------- # STEP 4: # -------------------------------------------- -weights <- Subset(weight$data, along = 'monthly', indices = c(11,12,1:1:6)) -slope <- Subset(slope, along = 'monthly', indices = c(11,12,1:1:6), drop = 'non-selected') +weights <- Subset(weight$data, along = 'monthly', indices = c(11, 12, 1:6)) +slope <- Subset(slope, along = 'monthly', indices = c(11, 12, 1:6), + drop = 'non-selected') fs <- CST_RainFARM(exp.qm, nf = 4, weights = weights, slope = slope, kmin = 1, nens = 10, verbose = TRUE, @@ -118,42 +119,47 @@ CST_SaveExp(newfs, destination = dir_output) # -------------------------------------------- library(s2dv) agg_png(paste0(dir_output, "EXP_11dec.png"), - width = 1000, height = 1000, units = 'px',res = 144) + width = 1000, height = 1100, units = 'px',res = 144) PlotEquiMap(exp$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, filled.continents = FALSE, bar_limits = c(0,40), - intylat = 2, intxlon = 2, title_scale = 0.7, + intylat = 2, intxlon = 2, title_scale = 1, + triangle_ends = c(TRUE, FALSE), toptitle = 'ECMWF-S5C3S', units = 'precipitation (mm)') dev.off() agg_png(paste0(dir_output, "EXPQM_11dec.png"), - width = 1000, height = 1000, units = 'px',res = 144) + width = 1000, height = 1100, units = 'px',res = 144) PlotEquiMap(exp.qm$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, filled.continents = FALSE, bar_limits = c(0,40), - intylat = 2, intxlon = 2, title_scale = 0.7, + intylat = 2, intxlon = 2, title_scale = 1, + triangle_ends = c(TRUE, FALSE), toptitle = 'Bias Adjusted', units = 'precipitation (mm)') dev.off() -agg_png(paste0(dir_output, "Down_11dec.png"), - width = 1000, height = 1000, units = 'px',res = 144) +agg_png(paste0(dir_output, "RF4_Down_11dec.png"), + width = 1000, height = 1100, units = 'px',res = 144) PlotEquiMap(fs$data[1,1,1,11,,,2],lon = fs$lon, lat = fs$lat, filled.continents = FALSE, bar_limits = c(0,40), - intylat = 2, intxlon = 2, title_scale = 0.7, - toptitle = 'Downsacaled', units = 'precipitation (mm)') + intylat = 2, intxlon = 2, title_scale = 1, + triangle_ends = c(TRUE, FALSE), + toptitle = 'Downsacaled RF 4', units = 'precipitation (mm)') dev.off() -agg_png(paste0(dir_output, "WeightsDec.png"), - width = 1000, height = 1000, units = 'px',res = 144) +agg_png(paste0(dir_output, "RF4_WeightsDec.png"), + width = 1000, height = 1100, units = 'px',res = 144) PlotEquiMap(weight$data[,,12], lon = weight$lon, lat = weight$lat, - filled.continents = FALSE, + filled.continents = FALSE, title_scale = 1, intylat = 2, intxlon = 2, - toptitle = 'December Weights WorldClim2') + toptitle = 'December Weights RF 4') dev.off() agg_png(paste0(dir_output, "Slope.png"), - width = 1000, height = 1000, units = 'px',res = 144) -plot(1:12, slope_plot, type = 'b', main = 'slope', pch = 3) -line(12, slope_plot[12], type = 'p', col = 'red', pch = 3) + width = 700, height = 700, units = 'px',res = 144) +plot(1:12, slope_plot, type = 'b', main = 'Slope', pch = 16, xlab = 'month', + ylab = 'Slope', bty = 'n') +lines(12, slope_plot[12], type = 'p', col = 'red', pch = 16) dev.off() # Plot ForecastPDF +library(abind) obsts <- MeanDims(obs$data, c('lat', 'lon'), na.rm = T) -print(quantile(obsts, c(0.1,0.3,0.6,0.9), na.rm = T)) +print(quantile(obsts, c(0.1, 0.3, 0.6, 0.9), na.rm = T)) expts <- MeanDims(exp$data, c('lat', 'lon'), na.rm = T) exp.qmts <- MeanDims(exp.qm$data, c('lat', 'lon'), na.rm = T) empty <- array(NA, c(dataset = 2, member = 225, sdate = 26, ftime = 31, monthly = 8)) @@ -161,15 +167,18 @@ data <- abind(expts, exp.qmts, along = 1) names(dim(data)) <- names(dim(expts)) data <- abind(data, empty, along = 2) names(dim(data)) <- names(dim(expts)) -fsts <- MeanDims(fs, c('lat', 'lon'), na.rm = T) +fsts <- MeanDims(fs$data, c('lat', 'lon'), na.rm = T) data <- abind(data, fsts, along = 1) names(dim(data)) <- c('dataset', 'members', 'sdate', 'ftime', 'monthly') -agg_png(paste0(dir_output, "/FiguresPDF_11DecemberAll2.png", - width = 750, height = 650, units = 'px', res = 144) -PlotForecastPDF(data[,,1,11,2], tercile.limits = c(0.58, 2.4), obs = obsts[,,1,11,2], - extreme.limits = c(0.06, 7.26), color.set = 'hydro', add.ensmemb = 'no', - var.name = "Precipitation (mm)", title = "Forecasts issued on Nov 1993 for 11th December 1993", - fcst.names = c("ECMWFS5C3S", "Bias Adjusted", "Downscaled")) +agg_png(paste0(dir_output, "/FiguresPDF_RF4_11DecemberAll2.png"), + width = 1000, height = 700, units = 'px', res = 144) +PlotForecastPDF(data[,,1,11,2], tercile.limits = c(0.67, 2.5), + #obs = obsts[,,1,11,2], + extreme.limits = c(0.09, 7.3), color.set = 'hydro', + #add.ensmemb = 'no', + var.name = "Precipitation (mm)", + title = "Forecasts issued on Nov 1993 for 11th December 1993", + fcst.names = c("ECMWFS5C3S", "Bias Adjusted", "Downscaled RF 4")) dev.off() diff --git a/man/CST_SaveExp.Rd b/man/CST_SaveExp.Rd index ddd9164e..92c97282 100644 --- a/man/CST_SaveExp.Rd +++ b/man/CST_SaveExp.Rd @@ -5,7 +5,7 @@ \title{Save CSTools objects of class 's2dv_cube' containing experiments or observed data in NetCDF format} \usage{ -CST_SaveExp(data, destination = "./CST_Data") +CST_SaveExp(data, destination = "./CST_Data", extra_string = NULL) } \arguments{ \item{data}{an object of class \code{s2dv_cube}.} @@ -15,6 +15,8 @@ to save the data. NetCDF file for each starting date are saved into the folder tree: destination/experiment/variable/. By default the function creates and saves the data into the folder "CST_Data" in the working directory.} + +\item{extra_string}{a character string to be include as part of the file name, for instance, to identify member or realization. It would be added to the file name between underscore characters.} } \description{ This function allows to divide and save a object of class diff --git a/man/SaveExp.Rd b/man/SaveExp.Rd index 40ace2db..345974b2 100644 --- a/man/SaveExp.Rd +++ b/man/SaveExp.Rd @@ -15,7 +15,8 @@ SaveExp( Dates, cdo_grid_name, projection, - destination + destination, + extra_string = NULL ) } \arguments{ @@ -40,9 +41,11 @@ SaveExp( \item{projection}{a character string indicating the projection name} \item{destination}{a character string indicating the path where to store the NetCDF files} + +\item{extra_string}{a character string to be include as part of the file name, for instance, to identify member or realization.} } \value{ -the function creates as many files as sdates per dataset. Each file could contain multiple members +the function creates as many files as sdates per dataset. Each file could contain multiple members. It would be added to the file name between underscore characters. The path will be created with the name of the variable and each Datasets. } \description{ -- GitLab From 8824d09546961f4f36ef427013d8248f2955bdf1 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 14 Jun 2021 17:08:23 +0200 Subject: [PATCH 04/11] remove commented lines --- ...Case2_PrecipitationDownscaling_RainFARM_RF4.R | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R index 44b1a94e..5af7770c 100644 --- a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R +++ b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R @@ -69,7 +69,7 @@ exp.qm <- CST_QuantileMapping(exp, obs, method = "QUANT", sample_dims = c('member', 'sdate', 'ftime'), ncores = 4) save(exp.qm, file = paste0(dir_output, 'ExpQM.RDS'), version = 2) - +save(exp, file = paste0(dir_output, 'Exp.RDS'), version = 2) # -------------------------------------------- # STEP 3: # -------------------------------------------- @@ -112,7 +112,7 @@ newfs <- CST_MergeDims(fs, merge_dims = c("ftime", "monthly"), na.rm = TRUE) newfs$Dates[[1]] <- exp$Dates[[1]] -CST_SaveExp(newfs, destination = dir_output) +CST_SaveExp(newfs, destination = paste0(dir_output, 'RF4/')) # -------------------------------------------- # Visualization @@ -124,7 +124,7 @@ PlotEquiMap(exp$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, filled.continents = FALSE, bar_limits = c(0,40), intylat = 2, intxlon = 2, title_scale = 1, triangle_ends = c(TRUE, FALSE), - toptitle = 'ECMWF-S5C3S', units = 'precipitation (mm)') + toptitle = 'SEAS5', units = 'precipitation (mm)') dev.off() agg_png(paste0(dir_output, "EXPQM_11dec.png"), width = 1000, height = 1100, units = 'px',res = 144) @@ -140,14 +140,14 @@ PlotEquiMap(fs$data[1,1,1,11,,,2],lon = fs$lon, lat = fs$lat, filled.continents = FALSE, bar_limits = c(0,40), intylat = 2, intxlon = 2, title_scale = 1, triangle_ends = c(TRUE, FALSE), - toptitle = 'Downsacaled RF 4', units = 'precipitation (mm)') + toptitle = 'Downscaled nf 4', units = 'precipitation (mm)') dev.off() agg_png(paste0(dir_output, "RF4_WeightsDec.png"), width = 1000, height = 1100, units = 'px',res = 144) PlotEquiMap(weight$data[,,12], lon = weight$lon, lat = weight$lat, filled.continents = FALSE, title_scale = 1, intylat = 2, intxlon = 2, - toptitle = 'December Weights RF 4') + toptitle = 'December Weights nf 4') dev.off() agg_png(paste0(dir_output, "Slope.png"), width = 700, height = 700, units = 'px',res = 144) @@ -171,14 +171,12 @@ fsts <- MeanDims(fs$data, c('lat', 'lon'), na.rm = T) data <- abind(data, fsts, along = 1) names(dim(data)) <- c('dataset', 'members', 'sdate', 'ftime', 'monthly') agg_png(paste0(dir_output, "/FiguresPDF_RF4_11DecemberAll2.png"), - width = 1000, height = 700, units = 'px', res = 144) + width = 1100, height = 700, units = 'px', res = 144) PlotForecastPDF(data[,,1,11,2], tercile.limits = c(0.67, 2.5), - #obs = obsts[,,1,11,2], extreme.limits = c(0.09, 7.3), color.set = 'hydro', - #add.ensmemb = 'no', var.name = "Precipitation (mm)", title = "Forecasts issued on Nov 1993 for 11th December 1993", - fcst.names = c("ECMWFS5C3S", "Bias Adjusted", "Downscaled RF 4")) + fcst.names = c("SEAS5", "Bias Adjusted", "Downscaled nf 4")) dev.off() -- GitLab From becb13727e7179a6556182b84722cd9666177870 Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 18 Jun 2021 17:52:24 +0200 Subject: [PATCH 05/11] Fix in PlotCombinedMap to scale bar titles for Use Case 1 --- DESCRIPTION | 2 +- R/PlotCombinedMap.R | 4 +- ...e2_PrecipitationDownscaling_RainFARM_RF4.R | 52 +++++++++---------- man/PlotCombinedMap.Rd | 3 ++ 4 files changed, 33 insertions(+), 28 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index aecd89b3..97a80c72 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -84,4 +84,4 @@ VignetteBuilder: knitr License: Apache License 2.0 Encoding: UTF-8 LazyData: true -RoxygenNote: 7.0.2 +RoxygenNote: 7.0.1 diff --git a/R/PlotCombinedMap.R b/R/PlotCombinedMap.R index 7169e4c4..8af4a14d 100644 --- a/R/PlotCombinedMap.R +++ b/R/PlotCombinedMap.R @@ -22,6 +22,7 @@ #' layers via the parameter 'dot_symbol'. #'@param bar_titles Optional vector of character strings providing the titles to be shown on top of each of the colour bars. #'@param legend_scale Scale factor for the size of the colour bar labels. Takes 1 by default. +#'@param cex_bar_titles Scale factor for the sizes of the bar titles. Takes 1.5 by default. #'@param fileout File where to save the plot. If not specified (default) a graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp and tiff #'@param width File width, in the units specified in the parameter size_units (inches by default). Takes 8 by default. #'@param height File height, in the units specified in the parameter size_units (inches by default). Takes 5 by default. @@ -69,6 +70,7 @@ PlotCombinedMap <- function(maps, lon, lat, mask = NULL, col_mask = 'grey', dots = NULL, bar_titles = NULL, legend_scale = 1, + cex_bar_titles = 1.5, fileout = NULL, width = 8, height = 5, size_units = 'in', res = 100, ...) { @@ -422,7 +424,7 @@ PlotCombinedMap <- function(maps, lon, lat, draw_separators = TRUE, extra_margin = c(2, 0, 2, 0), label_scale = legend_scale * 1.5) if (!is.null(bar_titles)) { - mtext(bar_titles[[k]], 3, line = -3, cex = 1.5) + mtext(bar_titles[[k]], 3, line = -3, cex = cex_bar_titles) } } diff --git a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R index 5af7770c..a4c4ed7a 100644 --- a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R +++ b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R @@ -24,7 +24,6 @@ library(ClimProjDiags) library(zeallot) library(ragg) dir_output <- '/esarchive/scratch/nperez/CSTools_manuscript/v20210603/' #slash end - # -------------------------------------------- # STEP 1: # -------------------------------------------- @@ -119,60 +118,61 @@ CST_SaveExp(newfs, destination = paste0(dir_output, 'RF4/')) # -------------------------------------------- library(s2dv) agg_png(paste0(dir_output, "EXP_11dec.png"), - width = 1000, height = 1100, units = 'px',res = 144) + width = 800, height = 900, units = 'px',res = 144) PlotEquiMap(exp$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, filled.continents = FALSE, bar_limits = c(0,40), - intylat = 2, intxlon = 2, title_scale = 1, - triangle_ends = c(TRUE, FALSE), + intylat = 2, intxlon = 2, title_scale = 0.8, + triangle_ends = c(TRUE, FALSE), degree_sym = TRUE, units_scale = 1.8, toptitle = 'SEAS5', units = 'precipitation (mm)') dev.off() agg_png(paste0(dir_output, "EXPQM_11dec.png"), - width = 1000, height = 1100, units = 'px',res = 144) + width = 800, height = 900, units = 'px',res = 144) PlotEquiMap(exp.qm$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, filled.continents = FALSE, bar_limits = c(0,40), - intylat = 2, intxlon = 2, title_scale = 1, - triangle_ends = c(TRUE, FALSE), + intylat = 2, intxlon = 2, title_scale = 0.8, + triangle_ends = c(TRUE, FALSE), degree_sym = TRUE, units_scale = 1.8, toptitle = 'Bias Adjusted', units = 'precipitation (mm)') dev.off() agg_png(paste0(dir_output, "RF4_Down_11dec.png"), - width = 1000, height = 1100, units = 'px',res = 144) + width = 800, height = 900, units = 'px',res = 144) PlotEquiMap(fs$data[1,1,1,11,,,2],lon = fs$lon, lat = fs$lat, filled.continents = FALSE, bar_limits = c(0,40), - intylat = 2, intxlon = 2, title_scale = 1, - triangle_ends = c(TRUE, FALSE), + intylat = 2, intxlon = 2, title_scale = 0.8, + triangle_ends = c(TRUE, FALSE), degree_sym = TRUE, units_scale = 1.8, toptitle = 'Downscaled nf 4', units = 'precipitation (mm)') dev.off() agg_png(paste0(dir_output, "RF4_WeightsDec.png"), - width = 1000, height = 1100, units = 'px',res = 144) + width = 800, height = 900, units = 'px',res = 144) PlotEquiMap(weight$data[,,12], lon = weight$lon, lat = weight$lat, - filled.continents = FALSE, title_scale = 1, - intylat = 2, intxlon = 2, + filled.continents = FALSE, title_scale = 0.8, + intylat = 2, intxlon = 2, degree_sym = TRUE, toptitle = 'December Weights nf 4') dev.off() agg_png(paste0(dir_output, "Slope.png"), width = 700, height = 700, units = 'px',res = 144) plot(1:12, slope_plot, type = 'b', main = 'Slope', pch = 16, xlab = 'month', - ylab = 'Slope', bty = 'n') + ylab = 'Slope', bty = 'n', cex.main = 1.2, cex.lab = 1.2) lines(12, slope_plot[12], type = 'p', col = 'red', pch = 16) dev.off() # Plot ForecastPDF library(abind) -obsts <- MeanDims(obs$data, c('lat', 'lon'), na.rm = T) -print(quantile(obsts, c(0.1, 0.3, 0.6, 0.9), na.rm = T)) -expts <- MeanDims(exp$data, c('lat', 'lon'), na.rm = T) -exp.qmts <- MeanDims(exp.qm$data, c('lat', 'lon'), na.rm = T) -empty <- array(NA, c(dataset = 2, member = 225, sdate = 26, ftime = 31, monthly = 8)) -data <- abind(expts, exp.qmts, along = 1) -names(dim(data)) <- names(dim(expts)) +#obsts <- MeanDims(obs$data, c('lat', 'lon'), na.rm = T) +#print(quantile(obsts, c(0.1, 0.3, 0.6, 0.9), na.rm = T)) +#expts <- MeanDims(exp$data, c('lat', 'lon'), na.rm = T) +#exp.qmts <- MeanDims(exp.qm$data, c('lat', 'lon'), na.rm = T) +data <- rbind(exp$data[1, ,1,11,5,4,2], exp.qm$data[1,,1,11,5,4,2]) +empty <- array(NA, c(dataset = 2, members = 225)) +names(dim(data)) <- names(dim(empty)) data <- abind(data, empty, along = 2) -names(dim(data)) <- names(dim(expts)) -fsts <- MeanDims(fs$data, c('lat', 'lon'), na.rm = T) -data <- abind(data, fsts, along = 1) -names(dim(data)) <- c('dataset', 'members', 'sdate', 'ftime', 'monthly') +names(dim(data)) <- names(dim(empty)) +data <- abind(data, fs$data[1,,1,11,15,18 ,2], along = 1) +names(dim(data)) <- names(dim(empty)) +print(dim(data)) +print(quantile(obs$data[1,,,11,5,4,2])) agg_png(paste0(dir_output, "/FiguresPDF_RF4_11DecemberAll2.png"), width = 1100, height = 700, units = 'px', res = 144) -PlotForecastPDF(data[,,1,11,2], tercile.limits = c(0.67, 2.5), +PlotForecastPDF(data, tercile.limits = c(0.67, 2.5), extreme.limits = c(0.09, 7.3), color.set = 'hydro', var.name = "Precipitation (mm)", title = "Forecasts issued on Nov 1993 for 11th December 1993", diff --git a/man/PlotCombinedMap.Rd b/man/PlotCombinedMap.Rd index e631761e..3d6661e1 100644 --- a/man/PlotCombinedMap.Rd +++ b/man/PlotCombinedMap.Rd @@ -19,6 +19,7 @@ PlotCombinedMap( dots = NULL, bar_titles = NULL, legend_scale = 1, + cex_bar_titles = 1.5, fileout = NULL, width = 8, height = 5, @@ -61,6 +62,8 @@ layers via the parameter 'dot_symbol'.} \item{legend_scale}{Scale factor for the size of the colour bar labels. Takes 1 by default.} +\item{cex_bar_titles}{Scale factor for the sizes of the bar titles. Takes 1.5 by default.} + \item{fileout}{File where to save the plot. If not specified (default) a graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp and tiff} \item{width}{File width, in the units specified in the parameter size_units (inches by default). Takes 8 by default.} -- GitLab From 8420b6000e226e4084c00ff6f52d0a3544c57893 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 13 Sep 2021 11:50:57 +0200 Subject: [PATCH 06/11] Latests changes use cases --- inst/doc/UseCase1_WindEvent_March2018.R | 30 +++++++----- ...e2_PrecipitationDownscaling_RainFARM_RF4.R | 49 +++++++++++++++---- ...aunch_UseCase2_PrecipitationDownscaling.sh | 2 +- 3 files changed, 58 insertions(+), 23 deletions(-) diff --git a/inst/doc/UseCase1_WindEvent_March2018.R b/inst/doc/UseCase1_WindEvent_March2018.R index 719f5bf3..d3bbb936 100644 --- a/inst/doc/UseCase1_WindEvent_March2018.R +++ b/inst/doc/UseCase1_WindEvent_March2018.R @@ -11,11 +11,11 @@ rm(list=ls()); gc(); # This code includes the bias adjustent and the results visualization # ---------------------------------------- -library(CSTools) +#library(CSTools) library(s2dv) library(ragg) library(multiApply) -output_dir <- "/esarchive/scratch/nperez/CSTools_manuscript/Wind/" +output_dir <- "/esarchive/scratch/nperez/CSTools_manuscript/v20210603/" exp_path <- list(name = "ECMWFS5", @@ -138,15 +138,17 @@ for (mm in 1:3) { wind_thres_latlon <- abind::abind(wind_thres_latlon, thres, along = 4) source("/esarchive/scratch/nperez/git/cstools/R/PlotCombinedMap.R") source("/esarchive/scratch/nperez/git/cstools/R/PlotMostLikelyQuantileMap.R") - agg_png(paste0(output_dir, "MostLikely_", mm, "_obstercile.png"), - width = 1000, height = 1000, units = 'px', res = 144) - PlotMostLikelyQuantileMap(probs = Mean_PB, lon = wind_fsct$lon, lat = wind_fsct$lat, + agg_png(paste0(output_dir, "Wind_MostLikely_", mm, "_obstercile.png"), + width = 1050, height = 1000, units = 'px', res = 144) + PlotMostLikelyQuantileMap(probs = Mean_PB, lon = wind_fsct$lon, + lat = wind_fsct$lat, sizetit = 1.5, intylat = 2, intxlon = 2, coast_width = 1.5, legend_scale = 0.8, - cat_dim = 'bin', dot_size = 2, + cat_dim = 'bin', dot_size = 2.5, + axes_label_scale = 1.6, degree_sym = TRUE, dots = filtered_obs_terciles[,,1,1,1,1], toptitle = c(paste0("Initialized on ", - month.name[as.numeric(months_in_advance[mm])]))) + month.name[as.numeric(months_in_advance[mm])]))) dev.off() } @@ -154,7 +156,7 @@ visual <- data.frame(dec = as.vector(MeanDims(wind_fsct_BC[[3]]$data, c('lat', ' jan = as.vector(MeanDims(wind_fsct_BC[[2]]$data, c('lat', 'lon'))), feb = as.vector(MeanDims(wind_fsct_BC[[1]]$data, c('lat', 'lon')))) - wind_obs_areave <- CST_Load(var = 'windagl100', obs = list(obs_path), + wind_obs_areave <- CSTools::CST_Load(var = 'windagl100', obs = list(obs_path), sdates = '20180301', nmember = 1, leadtimemin = 1, leadtimemax = 1, storefreq = "monthly", sampleperiod = 1, @@ -165,12 +167,14 @@ visual <- data.frame(dec = as.vector(MeanDims(wind_fsct_BC[[3]]$data, c('lat', ' print("IS DATA LOADED") print("Wait") -agg_png(paste0(output_dir, "PlotForecast_IP.png"), - width = 1000, height = 1000, units = 'px',res = 144) -PlotForecastPDF(visual, tercile.limits = wind_ref_terciles, +agg_png(paste0(output_dir, "Wind_PlotForecast_IP.png"), + width = 1000, height = 1000, units = 'px',res = 150) +CSTools::PlotForecastPDF(visual, tercile.limits = wind_ref_terciles, extreme.limits = wind_ref_extremes, - var.name = "Wind Speed 100 m (m/s)", title = "Bias Corrected forecasts at IP", - fcst.names = c("December", "January", "February"), obs = as.vector(wind_obs_areave$data)) + var.name = "Wind Speed 100 m (m/s)", + title = "Bias Corrected forecasts at IP", + fcst.names = c("December", "January", "February"), + obs = as.vector(wind_obs_areave$data)) dev.off() # Plotting observed terciles: diff --git a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R index a4c4ed7a..ee2df194 100644 --- a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R +++ b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R @@ -98,7 +98,7 @@ load(paste0(dir_output, 'weightsRF4.RDS')) # -------------------------------------------- # STEP 4: # -------------------------------------------- - +Rprof() weights <- Subset(weight$data, along = 'monthly', indices = c(11, 12, 1:6)) slope <- Subset(slope, along = 'monthly', indices = c(11, 12, 1:6), drop = 'non-selected') @@ -113,6 +113,8 @@ newfs <- CST_MergeDims(fs, merge_dims = c("ftime", "monthly"), newfs$Dates[[1]] <- exp$Dates[[1]] CST_SaveExp(newfs, destination = paste0(dir_output, 'RF4/')) +Rprof(NULL) +profile.info <- summaryRprof(paste0(dir_output, "Rprof.out")) # -------------------------------------------- # Visualization # -------------------------------------------- @@ -122,7 +124,8 @@ agg_png(paste0(dir_output, "EXP_11dec.png"), PlotEquiMap(exp$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, filled.continents = FALSE, bar_limits = c(0,40), intylat = 2, intxlon = 2, title_scale = 0.8, - triangle_ends = c(TRUE, FALSE), degree_sym = TRUE, units_scale = 1.8, + bar_label_scale = 1.3, axes_label_scale = 1.2, + triangle_ends = c(TRUE, FALSE), degree_sym = TRUE, units_scale = 1.5, toptitle = 'SEAS5', units = 'precipitation (mm)') dev.off() agg_png(paste0(dir_output, "EXPQM_11dec.png"), @@ -130,7 +133,8 @@ agg_png(paste0(dir_output, "EXPQM_11dec.png"), PlotEquiMap(exp.qm$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, filled.continents = FALSE, bar_limits = c(0,40), intylat = 2, intxlon = 2, title_scale = 0.8, - triangle_ends = c(TRUE, FALSE), degree_sym = TRUE, units_scale = 1.8, + bar_label_scale = 1.3, axes_label_scale = 1.2, + triangle_ends = c(TRUE, FALSE), degree_sym = TRUE, units_scale = 1.5, toptitle = 'Bias Adjusted', units = 'precipitation (mm)') dev.off() agg_png(paste0(dir_output, "RF4_Down_11dec.png"), @@ -138,21 +142,23 @@ agg_png(paste0(dir_output, "RF4_Down_11dec.png"), PlotEquiMap(fs$data[1,1,1,11,,,2],lon = fs$lon, lat = fs$lat, filled.continents = FALSE, bar_limits = c(0,40), intylat = 2, intxlon = 2, title_scale = 0.8, - triangle_ends = c(TRUE, FALSE), degree_sym = TRUE, units_scale = 1.8, + bar_label_scale = 1.3, axes_label_scale = 1.2, + triangle_ends = c(TRUE, TRUE), degree_sym = TRUE, units_scale = 1.5, toptitle = 'Downscaled nf 4', units = 'precipitation (mm)') dev.off() agg_png(paste0(dir_output, "RF4_WeightsDec.png"), width = 800, height = 900, units = 'px',res = 144) PlotEquiMap(weight$data[,,12], lon = weight$lon, lat = weight$lat, filled.continents = FALSE, title_scale = 0.8, + bar_label_scale = 1.3, axes_label_scale = 1.2, intylat = 2, intxlon = 2, degree_sym = TRUE, toptitle = 'December Weights nf 4') dev.off() agg_png(paste0(dir_output, "Slope.png"), width = 700, height = 700, units = 'px',res = 144) plot(1:12, slope_plot, type = 'b', main = 'Slope', pch = 16, xlab = 'month', - ylab = 'Slope', bty = 'n', cex.main = 1.2, cex.lab = 1.2) -lines(12, slope_plot[12], type = 'p', col = 'red', pch = 16) + ylab = 'Slope', bty = 'n', cex.main = 1.5, cex.lab = 1.3, cex = 1.5) +lines(12, slope_plot[12], type = 'p', col = 'red', pch = 16, cex = 1.3) dev.off() # Plot ForecastPDF @@ -161,6 +167,8 @@ library(abind) #print(quantile(obsts, c(0.1, 0.3, 0.6, 0.9), na.rm = T)) #expts <- MeanDims(exp$data, c('lat', 'lon'), na.rm = T) #exp.qmts <- MeanDims(exp.qm$data, c('lat', 'lon'), na.rm = T) +print("Quantiles gridpoint") +print(quantile(as.vector(exp.qm$data[1,,,11,5,4,2]), c(0.1,0.3,0.6,0.9))) data <- rbind(exp$data[1, ,1,11,5,4,2], exp.qm$data[1,,1,11,5,4,2]) empty <- array(NA, c(dataset = 2, members = 225)) names(dim(data)) <- names(dim(empty)) @@ -170,9 +178,30 @@ data <- abind(data, fs$data[1,,1,11,15,18 ,2], along = 1) names(dim(data)) <- names(dim(empty)) print(dim(data)) print(quantile(obs$data[1,,,11,5,4,2])) -agg_png(paste0(dir_output, "/FiguresPDF_RF4_11DecemberAll2.png"), - width = 1100, height = 700, units = 'px', res = 144) -PlotForecastPDF(data, tercile.limits = c(0.67, 2.5), +agg_png(paste0(dir_output, "/FiguresPDF_RF4_11December_GridPoint.png"), + width = 1000, height = 1000, units = 'px', res = 144) +PlotForecastPDF(data, tercile.limits = c(0.02, 1.17), + extreme.limits = c(0.0155555, 7.75), color.set = 'hydro', + var.name = "Precipitation (mm)", add.ensmemb = 'no', + title = "Forecasts issued on Nov 1993 for 11th December 1993", + fcst.names = c("SEAS5", "Bias Adjusted", "Downscaled nf 4")) +dev.off() + +obsts <- MeanDims(obs$data, c('lat', 'lon'), na.rm = T) +print(quantile(obsts, c(0.1, 0.3, 0.6, 0.9), na.rm = T)) +expts <- MeanDims(exp$data, c('lat', 'lon'), na.rm = T) +exp.qmts <- MeanDims(exp.qm$data, c('lat', 'lon'), na.rm = T) +empty <- array(NA, c(dataset = 2, member = 225, sdate = 26, ftime = 31, monthly = 8)) +data <- abind(expts, exp.qmts, along = 1) +names(dim(data)) <- names(dim(expts)) +data <- abind(data, empty, along = 2) +names(dim(data)) <- names(dim(expts)) +fsts <- MeanDims(fs$data, c('lat', 'lon'), na.rm = T) +data <- abind(data, fsts, along = 1) +names(dim(data)) <- c('dataset', 'members', 'sdate', 'ftime', 'monthly') +agg_png(paste0(dir_output, "/FiguresPDF_RF4_11December_Areave.png"), + width = 1400, height = 800, units = 'px', res = 144) +PlotForecastPDF(data[,,1,11,2], tercile.limits = c(0.67, 2.5), extreme.limits = c(0.09, 7.3), color.set = 'hydro', var.name = "Precipitation (mm)", title = "Forecasts issued on Nov 1993 for 11th December 1993", @@ -182,3 +211,5 @@ dev.off() + + diff --git a/inst/doc/launch_UseCase2_PrecipitationDownscaling.sh b/inst/doc/launch_UseCase2_PrecipitationDownscaling.sh index 55c3e848..2e1dd939 100644 --- a/inst/doc/launch_UseCase2_PrecipitationDownscaling.sh +++ b/inst/doc/launch_UseCase2_PrecipitationDownscaling.sh @@ -1,5 +1,5 @@ #!/bin/bash -#BSUB -W 2:00 +#BSUB -W 6:00 #BSUB -n 16 #BSUB -M 7000 #BSUB -J RainFARM_Downsc -- GitLab From b92f0291ea913d4e92a4975095b9f82559432443 Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 1 Oct 2021 12:37:17 +0200 Subject: [PATCH 07/11] Update --- ..._PrecipitationDownscaling_RainFARM_RF100.R | 37 +------------------ ..._UseCase2_PrecipitationDownscaling_RF4.sh} | 0 2 files changed, 1 insertion(+), 36 deletions(-) rename inst/doc/{launch_UseCase2_PrecipitationDownscaling.sh => launch_UseCase2_PrecipitationDownscaling_RF4.sh} (100%) diff --git a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R index 7b057180..146382dc 100644 --- a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R +++ b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R @@ -116,16 +116,9 @@ weights <- Subset(weight$data, along = 'monthly', indices = c(11, 12, 1:6)) slope <- Subset(slope, along = 'monthly', indices = c(11, 12, 1:6), drop = 'non-selected') k = 1 # To create the member ID -#----- -# To be removed when CSTools 4.0.1 is published: -source("/gpfs/scratch/bsc32/bsc32339/CSTools_manuscript/CST_SaveExp.R") -library(multiApply) -library(ncdf4) -library(s2dv) -#----- for (realizations in 1:10) { for (member in 1:25) { - result <- data # to store the data + result <- exp.qm # to store the data result$data <- NULL for (month in 1:8) { data <- exp.qm # to take the correct piece of data @@ -135,7 +128,6 @@ for (realizations in 1:10) { kmin = 1, nens = 1, verbose = TRUE, nprocs = 8, drop_realization = TRUE) - print(dim(fs$data)) result$data <- abind::abind(result$data, fs$data, along = 5) if (month == 2 & member == 1 & realization == 1) { # ---------------------------- @@ -149,33 +141,6 @@ for (realizations in 1:10) { triangle_ends = c(TRUE, FALSE), toptitle = 'Downsacaled RF 100', units = 'precipitation (mm)') dev.off() - # PlotForecastPDF - library(abind) - obsts <- MeanDims(obs$data, c('lat', 'lon'), na.rm = T) - print(quantile(obsts, c(0.1, 0.3, 0.6, 0.9), na.rm = T)) - expts <- MeanDims(exp$data, c('lat', 'lon'), na.rm = T) - exp.qmts <- MeanDims(exp.qm$data, c('lat', 'lon'), na.rm = T) - empty <- array(NA, c(dataset = 2, member = 225, sdate = 26, - ftime = 31, monthly = 8)) - data <- abind(expts, exp.qmts, along = 1) - names(dim(data)) <- names(dim(expts)) - data <- abind(data, empty, along = 2) - names(dim(data)) <- names(dim(expts)) - fsts <- MeanDims(fs$data, c('lat', 'lon'), na.rm = T) - print(dim(fsts)) - print(dim(data)) - data <- abind(data, fsts, along = 1) - names(dim(data)) <- c('dataset', 'members', 'sdate', 'ftime', 'monthly') - agg_png(paste0(dir_output, "/FiguresPDF_RF100_11DecemberAll2.png"), - width = 750, height = 650, units = 'px', res = 144) - PlotForecastPDF(data[,,1,11,2], tercile.limits = c(0.58, 2.4), - obs = obsts[,,1,11,2], - extreme.limits = c(0.06, 7.26), - color.set = 'hydro', add.ensmemb = 'no', - var.name = "Precipitation (mm)", - title = "Forecasts issued on Nov 1993 for 11th December 1993", - fcst.names = c("ECMWFS5C3S", "Bias Adjusted", "Downscaled")) - dev.off() } result$lon <- fs$lon result$lat <- fs$lat diff --git a/inst/doc/launch_UseCase2_PrecipitationDownscaling.sh b/inst/doc/launch_UseCase2_PrecipitationDownscaling_RF4.sh similarity index 100% rename from inst/doc/launch_UseCase2_PrecipitationDownscaling.sh rename to inst/doc/launch_UseCase2_PrecipitationDownscaling_RF4.sh -- GitLab From d2b117c82082dd6fd2d36cdd0fd4fcda926d1df8 Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 1 Oct 2021 12:41:51 +0200 Subject: [PATCH 08/11] Update code use case 3 --- .../UseCase3_data_preparation_SCHEME_model.R | 90 +++++++++++++++++-- 1 file changed, 83 insertions(+), 7 deletions(-) diff --git a/inst/doc/UseCase3_data_preparation_SCHEME_model.R b/inst/doc/UseCase3_data_preparation_SCHEME_model.R index 0ff104db..ada24ef2 100644 --- a/inst/doc/UseCase3_data_preparation_SCHEME_model.R +++ b/inst/doc/UseCase3_data_preparation_SCHEME_model.R @@ -1,3 +1,6 @@ +# Author: Bert Van Schaeybroeck +# Use Case 3: Seasonal forecasts for a river flow +# ----------------------------------------------- rm(list = ls()) library(CSTools) library(s2dverification) @@ -18,6 +21,7 @@ domain.high.res <- "greece_high_res" domain.low.res <- "greece_low_res" sdate.mon.to.use <- 5 #Month of startdate for ECMWF Sys5 possibilities are 5 (May) or 11 (November) sdate.day.to.use <- 1 #Day of startdate for ECMWF Sys5 only possibility is 1 +make.plot.msk <- F #mask to indicate if figures need to be made based on the output of the Analog function. #LOCAL PARAMETERS (to be adjusted for each working system) #--------------------------------------------------------- @@ -210,6 +214,10 @@ exp.msl.eur.merge <- CST_MergeDims( merge_dims = c("sdate", "ftime"), rename_dim = "sdate") +obs.msl.eur.merge.an <- CST_Anomaly(exp = obs.msl.eur.merge, dim_anom = 3) +exp.msl.eur.merge.an <- CST_Anomaly(exp = exp.msl.eur.merge, dim_anom = 3) + + #Load observational and forecast set of variable that needs to be calibrated and downscaled: file.to.load <- paste0(dir.rdata, "data_all.RData") if(file.exists(file.to.load)){ @@ -333,8 +341,8 @@ lon.low.res <- as.vector(cal.merge$lon) lat.low.res <- as.vector(cal.merge$lat) lon.high.res <- as.vector(obs.high.res$lon) lat.high.res <- as.vector(obs.high.res$lat) -lon.eur <- as.vector(obs.msl.eur.merge$lon) -lat.eur <- as.vector(obs.msl.eur.merge$lat) +lon.eur <- as.vector(obs.msl.eur.merge.an$lon) +lat.eur <- as.vector(obs.msl.eur.merge.an$lat) #amount of lead times in months. For ECMWF Sys5 it is 7: amt.lead.mon <- as.numeric(dim(cal.merge$data)["monthly"]) @@ -387,9 +395,9 @@ i.mbr.obs <- 1 for(i.mbr in seq(1, amt.mbr)){ for(i.mon in seq(1, amt.lead.mon)){ for(i.sdate in seq(1, amt.sdate)){ - i.mbr <- 1 - i.mon = 1 - i.sdate = 24 + #i.mbr <- 1 + #i.mon = 1 + #i.sdate = 24 date.to.use <- sub.time[ i.sdate, i.mon] date.an.lst <- sub.time[ , i.mon] cat("i.mbr = ", i.mbr, ", i.mon =", i.mon, ", i.sdate = ", @@ -400,10 +408,10 @@ for(i.mbr in seq(1, amt.mbr)){ exp.low.res.tmp <- exp.merge$data[i.dataset, i.mbr, i.sdate, , , i.mon] cal.low.res.tmp <- cal.merge$data[i.dataset, i.mbr, i.sdate, , , i.mon] #Extract the large-scale pressure field of that day - exp.msl.eur.tmp <- exp.msl.eur.merge$data[i.dataset, i.mbr, i.sdate, , , i.mon] + exp.msl.eur.tmp <- exp.msl.eur.merge.an$data[i.dataset, i.mbr, i.sdate, , , i.mon] #Extract all observations that will be used to find analogs - obs.msl.eur.tmp <- obs.msl.eur.merge$data[i.dataset, i.mbr.obs, , , , i.mon]#-i.sdate + obs.msl.eur.tmp <- obs.msl.eur.merge.an$data[i.dataset, i.mbr.obs, , , , i.mon]#-i.sdate obs.low.res.tmp <- obs.low.res.merge$data[i.dataset, i.mbr.obs, , , , i.mon] #-i.sdate obs.high.res.tmp <- obs.high.res.merge$data[i.dataset, i.mbr.obs, , , , i.mon] #-i.sdate names(dim(obs.high.res.tmp)) <- c("time", "lat", "lon") @@ -428,9 +436,77 @@ for(i.mbr in seq(1, amt.mbr)){ excludeTime = date.to.use, AnalogsInfo = T, nAnalogs = 1000) + + + if(make.plot.msk){ + corr.date <- as.character(res$dates[1]) #select the date of the most + corr.dex <- which(date.an.lst == corr.date) + + #The following figure shows the uncalibrated raw model field (analogous to Fig. 9a) + file.fig <- paste0("mbr_", i.mbr, "_mon_", i.mon, + "_sdate_", date.to.use, "_exp.low.res.pdf") + pdf(file = file.fig) + PlotEquiMap( + exp.low.res.tmp[ , ], + lon = obs.low.res.merge$lon, + lat = obs.low.res.merge$lat, + filled.continents = F, + intylat = 2, + intxlon = 2, + title_scale = 0.7, #bar_limits = c(0, 60), + units = "precipitation (mm)") + dev.off() + + #The following figure includes the calibrated model field (analogous to Fig. 9b) + file.fig <- paste0("mbr_", i.mbr, "_mon_", i.mon, + "_sdate_", date.to.use, "_cal.low.res.pdf") + pdf(file = file.fig) + PlotEquiMap( + cal.low.res.tmp, + lon = obs.low.res.merge$lon, + lat = obs.low.res.merge$lat, + filled.continents = F, + intylat = 2, + intxlon = 2, + title_scale = 0.7, #bar_limits = c(0, 60), + units = "precipitation (mm)") + + #The following figure includes the analog upscaled field (analogous to Fig. 9c) + file.fig <- paste0("mbr_", i.mbr, "_mon_", i.mon, + "_sdate_", date.to.use, "_obs.low.res.pdf") + pdf(file = file.fig) + PlotEquiMap( + obs.low.res.tmp[corr.dex, , ], + lon = obs.low.res.merge$lon, + lat = obs.low.res.merge$lat, + filled.continents = F, + intylat = 2, + intxlon = 2, + title_scale = 0.7, #bar_limits = c(0, 60), + units = "precipitation (mm)") + dev.off() + + #The following figure includes the analog field (analogous to Fig. 9d) + file.fig <- paste0("mbr_", i.mbr, "_mon_", i.mon, + "_sdate_", date.to.use, "_obs.high.res.pdf") + pdf(file = file.fig) + PlotEquiMap( + obs.high.res.tmp[corr.dex, , ], + lon = obs.high.res.merge$lon, + lat = obs.high.res.merge$lat, + filled.continents = F, + intylat = 2, + intxlon = 2, + title_scale = 0.7, #bar_limits = c(0, 60), + units = "precipitation (mm)") + dev.off() + + + } } } } } + -- GitLab From 69cc3b077a54fc7c84289e22437ae8b6a4eb57c8 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 4 Oct 2021 10:36:54 +0200 Subject: [PATCH 09/11] remove tests from the package --- .Rbuildignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.Rbuildignore b/.Rbuildignore index b2d8e5fc..fa596e70 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -5,4 +5,4 @@ ./.nc$ .*^(?!data)\.RData$ .*\.gitlab-ci.yml$ -#^tests$ +^tests$ -- GitLab From 7305bc8ad2013f171fb5477e73b46600320f2dc6 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 4 Oct 2021 10:52:39 +0200 Subject: [PATCH 10/11] news --- NEWS.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 154c6307..6d66f0a2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,11 +1,13 @@ ### CSTools 4.0.1 -<<<<<<< HEAD **Submission date to CRAN: XX-06-2021** - New features: + Dynamical Bias Correction method: `CST_ProxiesAttractors` and `CST_DynBiasCorrection` (optionally `Predictability`) + CST_BiasCorrection and BiasCorrection allows to calibrate a forecast given the calibration in the hindcast by using parameter 'exp_cor'. + + Use cases + + CST_SaveExp includes parameter extra_string + + PlotCombinedMap includes parameter cex_bar_titles - Fixes: + Calibration retains correlation absolute value -- GitLab From c972dba35a75ba30660edbe5e920a27b5eb3c5d7 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 4 Oct 2021 10:55:04 +0200 Subject: [PATCH 11/11] roxygen2 version installed in R 3.4.2 --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 318b382f..52999f2a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -85,4 +85,4 @@ VignetteBuilder: knitr License: Apache License 2.0 Encoding: UTF-8 LazyData: true -RoxygenNote: 7.0.1 +RoxygenNote: 7.0.2 -- GitLab