From b1d6f6ce40e3b8c19859557fe386cfc543ca9fb2 Mon Sep 17 00:00:00 2001 From: jhardenberg Date: Sat, 10 Aug 2019 23:30:44 +0200 Subject: [PATCH 1/6] first working version --- esmvaltool/diag_scripts/rainfarm/rainfarm.jl | 137 +++++++++++++++++++ esmvaltool/diag_scripts/shared/external.jl | 32 +++++ esmvaltool/recipes/recipe_rainfarm.yml | 2 +- 3 files changed, 170 insertions(+), 1 deletion(-) create mode 100644 esmvaltool/diag_scripts/rainfarm/rainfarm.jl create mode 100644 esmvaltool/diag_scripts/shared/external.jl diff --git a/esmvaltool/diag_scripts/rainfarm/rainfarm.jl b/esmvaltool/diag_scripts/rainfarm/rainfarm.jl new file mode 100644 index 0000000000..3271325b26 --- /dev/null +++ b/esmvaltool/diag_scripts/rainfarm/rainfarm.jl @@ -0,0 +1,137 @@ +# ############################################################################# +# rainfarm.jl +# Authors: J. von Hardenberg (ISAC-CNR, Italy) +# E. Arnone (ISAC-CNR, Italy) +# ############################################################################# +# Description +# ESMValTool diagnostic calling the RainFARM library written in Julia (by von Hardenberg, ISAC-CNR, Italy). +# RainFARM is a stochastic precipitation downscaling method, further adapted for climate downscaling. +# +# Required +# CDO +# Julia language: https://julialang.org +# RainFARM Julia library: https://github.com/jhardenberg/RainFARM.jl +# +# Optional +# +# Caveats +# +# Modification history +# 20190810 hard_jo: rewritten in pure Julia, no R +# 20181210 hard_jo: cleanup and using juliacall +# 20180508-A_arnone_e: Conversion to v2.0 +# 20170908-A_arnone_e: 1st github version +# +# ############################################################################ + +import YAML +using RainFARM +using Printf + +function provenance_record(infile) + xprov = Dict( "ancestors" => infile, + "authors" => ["hard_jo", "arno_en"], + "references" => ["donofrio14jh", "rebora06jhm", + "terzago18nhess"], + "projects" => ["c3s-magic"], + "caption" => "RainFARM precipitation downscaling", + "statistics" => ["other"], + "realms" => ["atmos"], + "themes" => ["phys"], + "domains" => ["reg"] + ) + return(xprov) +end + +let +diag_scripts_dir = ENV["diag_scripts"] +include(string(diag_scripts_dir, "/shared/external.jl")) + +settings = YAML.load_file(ARGS[1]) + +metadata = YAML.load_file(settings["input_files"][1]) +climofiles = collect(keys(metadata)) +climolist = metadata[climofiles[1]] +varname = climolist["short_name"] +diag_base = climolist["diagnostic"] + +println(diag_base, ": starting routine") +println(diag_base, ": creating work and plot directories") +work_dir = settings["work_dir"] +run_dir = settings["run_dir"] +mkpath(work_dir) +mkpath(run_dir) +cd(run_dir) + +# setup provenance file and list +provenance_file = string(run_dir, "/diagnostic_provenance.yml") +provenance = Dict() + +# Reading parameters from the settings +nf = get(settings, "nf", 2) +slope = get(settings, "slope", 0.0) +kmin = get(settings, "kmin", 1) +nens = get(settings, "nens", 1) +weights_climo = get(settings, "weights_climo", "") +conserv_glob = get(settings, "conserv_glob", false) +conserv_smooth = get(settings, "conserv_smooth", true) + +# WEIGHTS CLIMO MISSING +ww = 1. + + # Conservation options +if (conserv_glob) + println("Conserving global field") +elseif (conserv_smooth) + println("Smooth conservation") +else + println("Box conservation") +end + +for (infile, value) in metadata + println("SLOPE", settings["slope"]) + println("SLOPE", slope) + (infilename, ) = splitext(basename(infile)) + outfilename = string(work_dir, "/", infilename, "_downscaled") + + println(diag_base, ": calling RainFARM") + + (pr, lon_mat, lat_mat) = read_netcdf2d(infile, varname) + + # Ensure grid is square and with even dims + #nmin <- min(dim(pr)[1], dim(pr)[2]) + #nmin <- floor(nmin / 2) * 2 + #pr <- pr[1:nmin, 1:nmin, ] + #if (is.vector(lon_mat)) { + # lon_mat <- lon_mat[1:nmin] + # lat_mat <- lat_mat[1:nmin] + #} else { + # lon_mat <- lon_mat[1:nmin, 1:nmin] + # lat_mat <- lat_mat[1:nmin, 1:nmin] + #} + + (lon_f, lat_f) = lon_lat_fine(lon_mat, lat_mat, nf); + + # Automatic spectral slope + if (slope == 0.) + (fxp, ftp)=fft3d(pr) + slope =fitslopex(fxp, kmin=kmin) + println("Computed spatial spectral slope: ", slope) + else + println("Fixed spatial spectral slope: ", slope) + end + + for iens=1:nens + println("Realization ", iens) + rd=rainfarm(pr, slope, nf, ww, fglob = conserv_glob, + fsmooth = conserv_smooth, verbose=true) + fname=@sprintf("%s_%04d.nc", outfilename, iens); + write_netcdf2d(fname, rd, lon_f, lat_f, varname, infile) + xprov = provenance_record(infile) + provenance[fname] = xprov + end +end + +# Write provenance file +create_yaml(provenance, provenance_file) +end diff --git a/esmvaltool/diag_scripts/shared/external.jl b/esmvaltool/diag_scripts/shared/external.jl new file mode 100644 index 0000000000..cae963ef84 --- /dev/null +++ b/esmvaltool/diag_scripts/shared/external.jl @@ -0,0 +1,32 @@ +""" + create_yaml(filename, dict) +Write a dictionary to a YAML file. +""" +function create_yaml(dict::Dict{Any,Any}, filename::AbstractString) + os = open(filename, "w") + for (key, value) in dict + println(os, "? ",key) + indent = ": " + print_yaml(os, value, indent) + end + close(os) +end + +function print_yaml(os::IOStream, obj::Dict{String,Any}, indent::String) + for (key, value) in obj + print(os, indent, key, ": ") + print_yaml(os, value, " ") + indent = " " + end +end + +function print_yaml(os::IOStream, obj::String, indent::String) + println(os, obj) +end + +function print_yaml(os::IOStream, obj::Array{String}, indent::String) + println(os) + for i = 1:length(obj) + println(os, indent, "- ", obj[i]) + end +end diff --git a/esmvaltool/recipes/recipe_rainfarm.yml b/esmvaltool/recipes/recipe_rainfarm.yml index 07fb7aec4b..f771015a6d 100644 --- a/esmvaltool/recipes/recipe_rainfarm.yml +++ b/esmvaltool/recipes/recipe_rainfarm.yml @@ -42,7 +42,7 @@ diagnostics: mip: day scripts: rainfarm: - script: rainfarm/rainfarm.R + script: rainfarm/rainfarm.jl slope: 1.7 # spatial spectral slope (set to 0 to compute from large scales) nens: 2 # number of ensemble members to be calculated nf: 8 # subdivisions for downscaling From dcef9ce4e23a3a00d3ef612cc06d2769d29bd9de Mon Sep 17 00:00:00 2001 From: jhardenberg Date: Sun, 11 Aug 2019 10:52:51 +0200 Subject: [PATCH 2/6] adapt to square grid --- esmvaltool/diag_scripts/rainfarm/rainfarm.jl | 110 +++++++++---------- 1 file changed, 54 insertions(+), 56 deletions(-) diff --git a/esmvaltool/diag_scripts/rainfarm/rainfarm.jl b/esmvaltool/diag_scripts/rainfarm/rainfarm.jl index 3271325b26..90c716151d 100644 --- a/esmvaltool/diag_scripts/rainfarm/rainfarm.jl +++ b/esmvaltool/diag_scripts/rainfarm/rainfarm.jl @@ -29,18 +29,18 @@ using RainFARM using Printf function provenance_record(infile) - xprov = Dict( "ancestors" => infile, - "authors" => ["hard_jo", "arno_en"], - "references" => ["donofrio14jh", "rebora06jhm", - "terzago18nhess"], - "projects" => ["c3s-magic"], - "caption" => "RainFARM precipitation downscaling", - "statistics" => ["other"], - "realms" => ["atmos"], - "themes" => ["phys"], - "domains" => ["reg"] - ) - return(xprov) + xprov = Dict( "ancestors" => infile, + "authors" => ["hard_jo", "arno_en"], + "references" => ["donofrio14jh", "rebora06jhm", + "terzago18nhess"], + "projects" => ["c3s-magic"], + "caption" => "RainFARM precipitation downscaling", + "statistics" => ["other"], + "realms" => ["atmos"], + "themes" => ["phys"], + "domains" => ["reg"] + ) + return(xprov) end let @@ -81,55 +81,53 @@ ww = 1. # Conservation options if (conserv_glob) - println("Conserving global field") + println("Conserving global field") elseif (conserv_smooth) - println("Smooth conservation") + println("Smooth conservation") else - println("Box conservation") + println("Box conservation") end for (infile, value) in metadata - println("SLOPE", settings["slope"]) - println("SLOPE", slope) - (infilename, ) = splitext(basename(infile)) - outfilename = string(work_dir, "/", infilename, "_downscaled") - - println(diag_base, ": calling RainFARM") - - (pr, lon_mat, lat_mat) = read_netcdf2d(infile, varname) - - # Ensure grid is square and with even dims - #nmin <- min(dim(pr)[1], dim(pr)[2]) - #nmin <- floor(nmin / 2) * 2 - #pr <- pr[1:nmin, 1:nmin, ] - #if (is.vector(lon_mat)) { - # lon_mat <- lon_mat[1:nmin] - # lat_mat <- lat_mat[1:nmin] - #} else { - # lon_mat <- lon_mat[1:nmin, 1:nmin] - # lat_mat <- lat_mat[1:nmin, 1:nmin] - #} - - (lon_f, lat_f) = lon_lat_fine(lon_mat, lat_mat, nf); - - # Automatic spectral slope - if (slope == 0.) - (fxp, ftp)=fft3d(pr) - slope =fitslopex(fxp, kmin=kmin) - println("Computed spatial spectral slope: ", slope) - else - println("Fixed spatial spectral slope: ", slope) - end - - for iens=1:nens - println("Realization ", iens) - rd=rainfarm(pr, slope, nf, ww, fglob = conserv_glob, - fsmooth = conserv_smooth, verbose=true) - fname=@sprintf("%s_%04d.nc", outfilename, iens); - write_netcdf2d(fname, rd, lon_f, lat_f, varname, infile) - xprov = provenance_record(infile) - provenance[fname] = xprov - end + (infilename, ) = splitext(basename(infile)) + outfilename = string(work_dir, "/", infilename, "_downscaled") + + println(diag_base, ": calling RainFARM") + + (pr, lon_mat, lat_mat) = read_netcdf2d(infile, varname) + + # Ensure grid is square and with even dims + nmin = min(size(pr)[1], size(pr)[2]) + nmin = floor(Int, nmin / 2) * 2 + pr = pr[1:nmin, 1:nmin, :] + if (ndims(lon_mat) == 1) + lon_mat = lon_mat[1:nmin] + lat_mat = lat_mat[1:nmin] + else + lon_mat = lon_mat[1:nmin, 1:nmin] + lat_mat = lat_mat[1:nmin, 1:nmin] + end + + (lon_f, lat_f) = lon_lat_fine(lon_mat, lat_mat, nf); + + # Automatic spectral slope + if (slope == 0.) + (fxp, ftp)=fft3d(pr) + slope =fitslopex(fxp, kmin=kmin) + println("Computed spatial spectral slope: ", slope) + else + println("Fixed spatial spectral slope: ", slope) + end + + for iens=1:nens + println("Realization ", iens) + rd=rainfarm(pr, slope, nf, ww, fglob = conserv_glob, + fsmooth = conserv_smooth, verbose=true) + fname=@sprintf("%s_%04d.nc", outfilename, iens); + write_netcdf2d(fname, rd, lon_f, lat_f, varname, infile) + xprov = provenance_record(infile) + provenance[fname] = xprov + end end # Write provenance file From ca167e46ae6a7b8c7b4e378386a35b4de338a46d Mon Sep 17 00:00:00 2001 From: jhardenberg Date: Mon, 12 Aug 2019 12:20:48 +0200 Subject: [PATCH 3/6] added weights --- esmvaltool/diag_scripts/rainfarm/rainfarm.jl | 30 +++++++++++++++----- esmvaltool/recipes/recipe_rainfarm.yml | 2 +- 2 files changed, 24 insertions(+), 8 deletions(-) diff --git a/esmvaltool/diag_scripts/rainfarm/rainfarm.jl b/esmvaltool/diag_scripts/rainfarm/rainfarm.jl index 90c716151d..ce41d8ccf9 100644 --- a/esmvaltool/diag_scripts/rainfarm/rainfarm.jl +++ b/esmvaltool/diag_scripts/rainfarm/rainfarm.jl @@ -45,7 +45,7 @@ end let diag_scripts_dir = ENV["diag_scripts"] -include(string(diag_scripts_dir, "/shared/external.jl")) +include(joinpath(diag_scripts_dir, "shared/external.jl")) settings = YAML.load_file(ARGS[1]) @@ -64,7 +64,7 @@ mkpath(run_dir) cd(run_dir) # setup provenance file and list -provenance_file = string(run_dir, "/diagnostic_provenance.yml") +provenance_file = joinpath(run_dir, "diagnostic_provenance.yml") provenance = Dict() # Reading parameters from the settings @@ -75,11 +75,13 @@ nens = get(settings, "nens", 1) weights_climo = get(settings, "weights_climo", "") conserv_glob = get(settings, "conserv_glob", false) conserv_smooth = get(settings, "conserv_smooth", true) +auxiliary_data_dir = get(settings, "auxiliary_data_dir", "") -# WEIGHTS CLIMO MISSING -ww = 1. +if weights_climo isa Bool # Compatibility with old standard + weights_climo = "" +end - # Conservation options +# Conservation options if (conserv_glob) println("Conserving global field") elseif (conserv_smooth) @@ -90,9 +92,9 @@ end for (infile, value) in metadata (infilename, ) = splitext(basename(infile)) - outfilename = string(work_dir, "/", infilename, "_downscaled") + outfilename = joinpath(work_dir, infilename * "_downscaled") - println(diag_base, ": calling RainFARM") + println(diag_base, ": calling RainFARM for ", infilename) (pr, lon_mat, lat_mat) = read_netcdf2d(infile, varname) @@ -119,6 +121,20 @@ for (infile, value) in metadata println("Fixed spatial spectral slope: ", slope) end + if weights_climo != "" + if weights_climo[1] != '/' + weights_climo = joinpath(auxiliary_data_dir, weights_climo) + end + println("Using external climatology for weights: ", weights_climo) + fileweights = joinpath(work_dir, infilename * "_w.nc") + + ww = rfweights(weights_climo, infile, nf, + weightsfn = fileweights, varname = varname, + fsmooth = conserv_smooth) + else + ww = 1. + end + for iens=1:nens println("Realization ", iens) rd=rainfarm(pr, slope, nf, ww, fglob = conserv_glob, diff --git a/esmvaltool/recipes/recipe_rainfarm.yml b/esmvaltool/recipes/recipe_rainfarm.yml index f771015a6d..b2fed1bb4e 100644 --- a/esmvaltool/recipes/recipe_rainfarm.yml +++ b/esmvaltool/recipes/recipe_rainfarm.yml @@ -48,4 +48,4 @@ diagnostics: nf: 8 # subdivisions for downscaling conserv_glob: false # conserve precipitation over full domain (choose either glob or smooth, glob has priority) conserv_smooth: true # conserve precipitation using convolution (if neither is chosen box conservation is used) - weights_climo: false # orographic weights: set to false or full path to a fine-scale precipitation climatology file + weights_climo: false # orographic weights: set to false (or omit) or path to a fine-scale precipitation climatology file From 6219b1d95fa09964108c5e5d0a06ac115dec5d90 Mon Sep 17 00:00:00 2001 From: jhardenberg Date: Mon, 12 Aug 2019 12:21:48 +0200 Subject: [PATCH 4/6] JuliaCall not needed anymore --- esmvaltool/install/R/r_requirements.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/esmvaltool/install/R/r_requirements.txt b/esmvaltool/install/R/r_requirements.txt index 8510439047..bc3a3ad9f6 100644 --- a/esmvaltool/install/R/r_requirements.txt +++ b/esmvaltool/install/R/r_requirements.txt @@ -6,7 +6,6 @@ dotCall64 functional ggplot2 gridExtra -JuliaCall lintr logging mapproj From ab7eb596ee2d0fbd72a8ee3efca9479d3adf55f2 Mon Sep 17 00:00:00 2001 From: jhardenberg Date: Mon, 12 Aug 2019 12:22:28 +0200 Subject: [PATCH 5/6] updated docs --- doc/sphinx/source/recipes/recipe_rainfarm.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/sphinx/source/recipes/recipe_rainfarm.rst b/doc/sphinx/source/recipes/recipe_rainfarm.rst index 695a4a4313..ce8e687f9f 100644 --- a/doc/sphinx/source/recipes/recipe_rainfarm.rst +++ b/doc/sphinx/source/recipes/recipe_rainfarm.rst @@ -19,7 +19,7 @@ Recipes are stored in recipes/ Diagnostics are stored in diag_scripts/rainfarm/ -* rainfarm.R +* rainfarm.jl User settings @@ -32,7 +32,7 @@ User settings * nf: number of subdivisions for downscaling (e.g. 8 will produce output fields with linear resolution increased by a factor 8) * conserv_glob: logical, if to conserve precipitation over full domain * conserv_smooth: logical, if to conserve precipitation using convolution (if neither conserv_glob or conserv_smooth is chosen, box conservation is used) -* weights_climo: set to false if no orographic weights are to be used, else set it to the full path to a fine-scale precipitation climatology file. The file is expected to be in NetCDF format and should contain at least one precipitation field. If several fields at different times are provided, a climatology is derived by time averaging. Suitable climatology files could be for example a fine-scale precipitation climatology from a high-resolution regional climate model (see e.g. Terzago et al. 2018), a local high-resolution gridded climatology from observations, or a reconstruction such as those which can be downloaded from the WORLDCLIM (http://www.worldclim.org) or CHELSA (http://chelsa-climate.org) websites. The latter data will need to be converted to NetCDF format before being used (see for example the GDAL tools (https://www.gdal.org). +* weights_climo: set to false or omit if no orographic weights are to be used, else set it to the path to a fine-scale precipitation climatology file. If a relative file path is used, `auxiliary_data_dir` will be searched for this file. The file is expected to be in NetCDF format and should contain at least one precipitation field. If several fields at different times are provided, a climatology is derived by time averaging. Suitable climatology files could be for example a fine-scale precipitation climatology from a high-resolution regional climate model (see e.g. Terzago et al. 2018), a local high-resolution gridded climatology from observations, or a reconstruction such as those which can be downloaded from the WORLDCLIM (http://www.worldclim.org) or CHELSA (http://chelsa-climate.org) websites. The latter data will need to be converted to NetCDF format before being used (see for example the GDAL tools (https://www.gdal.org). Variables From 873ebe3786244fda3ac42382bc8e4938745e8863 Mon Sep 17 00:00:00 2001 From: Mattia Righi Date: Fri, 6 Sep 2019 12:16:21 +0200 Subject: [PATCH 6/6] Change author tags after new standard --- esmvaltool/diag_scripts/rainfarm/rainfarm.jl | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/esmvaltool/diag_scripts/rainfarm/rainfarm.jl b/esmvaltool/diag_scripts/rainfarm/rainfarm.jl index ce41d8ccf9..37e3992af6 100644 --- a/esmvaltool/diag_scripts/rainfarm/rainfarm.jl +++ b/esmvaltool/diag_scripts/rainfarm/rainfarm.jl @@ -4,8 +4,10 @@ # E. Arnone (ISAC-CNR, Italy) # ############################################################################# # Description -# ESMValTool diagnostic calling the RainFARM library written in Julia (by von Hardenberg, ISAC-CNR, Italy). -# RainFARM is a stochastic precipitation downscaling method, further adapted for climate downscaling. +# ESMValTool diagnostic calling the RainFARM library written in Julia +# (by von Hardenberg, ISAC-CNR, Italy). +# RainFARM is a stochastic precipitation downscaling method, further adapted +# for climate downscaling. # # Required # CDO @@ -17,10 +19,10 @@ # Caveats # # Modification history -# 20190810 hard_jo: rewritten in pure Julia, no R -# 20181210 hard_jo: cleanup and using juliacall -# 20180508-A_arnone_e: Conversion to v2.0 -# 20170908-A_arnone_e: 1st github version +# 20190810-vonhardenberg_jost: rewritten in pure Julia, no R +# 20181210-vonhardenberg_jost: cleanup and using juliacall +# 20180508-arnone_enrico: Conversion to v2.0 +# 20170908-arnone_enrico: 1st github version # # ############################################################################ @@ -30,7 +32,7 @@ using Printf function provenance_record(infile) xprov = Dict( "ancestors" => infile, - "authors" => ["hard_jo", "arno_en"], + "authors" => ["vonhardenberg_jost", "arnone_enrico"], "references" => ["donofrio14jh", "rebora06jhm", "terzago18nhess"], "projects" => ["c3s-magic"],