Merge pull request #1242 from ESMValGroup/version2_julia_rainfarm

RainFARM in pure Julia

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

esmvaltool/diag_scripts/rainfarm/rainfarm.jl

@@ -0,0 +1,153 @@
+# #############################################################################
+# 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-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
+#
+# ############################################################################
+
+import YAML
+using RainFARM
+using Printf
+
+function provenance_record(infile)
+    xprov = Dict( "ancestors" => infile,
+        "authors" => ["vonhardenberg_jost", "arnone_enrico"],
+        "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(joinpath(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 = joinpath(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) 
+auxiliary_data_dir = get(settings, "auxiliary_data_dir", "") 
+
+if weights_climo isa Bool # Compatibility with old standard
+    weights_climo = ""
+end
+
+# 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
+    (infilename, ) = splitext(basename(infile))
+    outfilename = joinpath(work_dir, infilename * "_downscaled")
+
+    println(diag_base, ": calling RainFARM for ", infilename)
+
+    (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
+
+    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,
+                    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

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

esmvaltool/install/R/r_requirements.txt

@@ -6,7 +6,6 @@ dotCall64
 functional
 ggplot2
 gridExtra
-JuliaCall
 lintr
 logging
 mapproj

esmvaltool/recipes/recipe_rainfarm.yml

@@ -43,10 +43,10 @@ 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
         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