Repository: Deltares/Wflow.jl
Pull request #74 (2020-12-18 18:35:54)

drop a few dependencies
BenchmarkTools was only needed for testing, so it becomes a test dependency only.
The second commit removes the CSV and DataFrames dependencies.
This speeds up loading a precompiled Wflow from 4.3 to 3.2 seconds on my machine, julia 1.6.



Patch:
diff --git a/Project.toml b/Project.toml
index c5478888e..2d7f6dace 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,11 +5,9 @@ version = "0.1.0"
 
 [deps]
 BasicModelInterface = "59605e27-edc0-445a-b93d-c09a3a50b330"
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CFTime = "179af706-886a-5703-950a-314cd64e0468"
-CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 FieldMetadata = "bf96fef3-21d2-5d20-8afa-0e7d4c32a885"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
@@ -22,10 +20,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
 BasicModelInterface = "0.1"
-BenchmarkTools = "0.5"
 CFTime = "0.1"
-CSV = "0.6, 0.7, 0.8"
-DataFrames = "0.21, 0.22"
 FieldMetadata = "0.3"
 LightGraphs = "1.3"
 NCDatasets = "0.10.4, 0.11"
@@ -36,9 +31,10 @@ UnPack = "1"
 julia = "1.5"
 
 [extras]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Polynomials", "Random", "Test"]
+test = ["BenchmarkTools", "Polynomials", "Random", "Test"]
diff --git a/docs/src/lateral/gwf.md b/docs/src/lateral/gwf.md
index ea7bfde39..ce69bc752 100644
--- a/docs/src/lateral/gwf.md
+++ b/docs/src/lateral/gwf.md
@@ -96,7 +96,6 @@ head [L] of the head boundary and  ``\phi`` the head of the aquifer cell.
 A volumetric well rate [L``^3`` T``^{-1}``] can be specified as a boundary condition.
 
 
-!!! note 
-    For an unconfined aquifer the boundary fluxes are checked, in case of a dry aquifer cell
-    a negative flux is not allowed.
+!!! note For an unconfined aquifer the boundary fluxes are checked, in case of a dry aquifer
+    cell a negative flux is not allowed.
 
diff --git a/src/Wflow.jl b/src/Wflow.jl
index a90a27536..e8f24d70f 100644
--- a/src/Wflow.jl
+++ b/src/Wflow.jl
@@ -7,12 +7,11 @@ using NCDatasets
 using StaticArrays
 using Statistics
 using UnPack
-using CSV
 using Random
-using DataFrames
 using BasicModelInterface
 using FieldMetadata
 using Parameters
+using DelimitedFiles
 
 @metadata get_units "mm"
 
diff --git a/src/io.jl b/src/io.jl
index 75d9d0408..eb56f277d 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -909,3 +909,28 @@ function rewind!(clock)
     clock.time -= clock.Δt
     return clock
 end
+
+"Read a rating curve from CSV into a NamedTuple of vectors"
+function read_sh_csv(path)
+    data, header = readdlm(path, ',', Float64, header = true)
+    names = vec(uppercase.(header))
+    idx_h = findfirst(==("H"), names)
+    idx_s = findfirst(==("S"), names)
+
+    if isnothing(idx_h) || isnothing(idx_s)
+        error("$path needs to provide H and S columns, got $names")
+    end
+
+    return (H = data[:, idx_h], S = data[:, idx_s])
+end
+
+"Read a specific storage curve from CSV into a NamedTuple of vectors"
+function read_hq_csv(path)
+    data = readdlm(path, ',', Float64, skipstart=1)
+    # Q is a matrix with 365 columns, one for each day in the year
+    return (H=data[:, 1], Q=data[:, 2:end])
+end
+
+# these represent the type of the rating curve and specific storage data
+const SH = NamedTuple{(:H, :S), Tuple{Vector{Float64}, Vector{Float64}}}
+const HQ = NamedTuple{(:H, :Q), Tuple{Vector{Float64}, Matrix{Float64}}}
diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl
index 515e4103b..44137b190 100644
--- a/src/reservoir_lake.jl
+++ b/src/reservoir_lake.jl
@@ -194,8 +194,8 @@ end
     outflowfunc::Vector{Int} | "-"          # type of lake rating curve, 1: Q = f(H) from lake data and interpolation, 2: General Q = b(H - H₀)ᵉ, 3: Case of Puls Approach Q = b(H - H₀)²
     b::Vector{T} | "m3/2 s-1 (if e=3/2)"    # rating curve coefficient
     e::Vector{T} | "-"                      # rating curve exponent
-    sh::Vector{DataFrame}                   # data for storage curve
-    hq::Vector{DataFrame}                   # data for rating curve
+    sh::Vector{Union{SH, Missing}}          # data for storage curve
+    hq::Vector{Union{HQ, Missing}}          # data for rating curve
     waterlevel::Vector{T} | "m"             # waterlevel H [m] of lake
     inflow::Vector{T} | "m3"                # inflow to the lake [m³]
     storage::Vector{T} | "m3"               # storage lake [m³]
@@ -316,9 +316,11 @@ function initialize_natural_lake(config, static_path, nc, inds_riv, nriv, pits)
     # length of all lake cells. To do that we need to introduce a mapping.
     n_lakes = length(inds_lake)
 
-    sh = Vector{DataFrame}(undef, n_lakes)
-    hq = Vector{DataFrame}(undef, n_lakes)
+    sh = Vector{Union{SH, Missing}}(missing, n_lakes)
+    hq = Vector{Union{HQ, Missing}}(missing, n_lakes)
     for i = 1:n_lakes
+        lakeloc = lakelocs[i]
+
         if linked_lakelocs[i] > 0
             linked_lakelocs[i] = i
         else
@@ -326,23 +328,13 @@ function initialize_natural_lake(config, static_path, nc, inds_riv, nriv, pits)
         end
 
         if lake_storfunc[i] == 2
-            sh[i] = CSV.read(
-                joinpath(static_path, "lake_sh_$(lakelocs[i])"),
-                DataFrame,
-                type = Float64,
-            )
-        else
-            sh[i] = DataFrame()
+            csv_path = joinpath(static_path, "lake_sh_$lakeloc.csv")
+            sh[i] = read_sh_csv(csv_path)
         end
 
         if lake_outflowfunc[i] == 1
-            hq[i] = CSV.read(
-                joinpath(static_path, "lake_hq_$(lakelocs[i])"),
-                DataFrame,
-                type = Float64,
-            )
-        else
-            hq[i] = DataFrame()
+            csv_path = joinpath(static_path, "lake_hq_$lakeloc.csv")
+            hq[i] = read_hq_csv(csv_path)
         end
 
         if lake_outflowfunc[i] == 3 && lake_storfunc[i] != 1
@@ -431,7 +423,7 @@ function update(lake::NaturalLake, i, inflow, doy, timestepsecs)
     lo = lake.lowerlake_ind[i]
     has_lowerlake = lo != 0
 
-    col = max(doy + 1, 366)
+    col = max(doy, 365)
     ### Modified Puls Approach (Burek et al., 2013, LISFLOOD) ###
     # outflowfunc = 3
     # Calculate lake factor and SI parameter
@@ -467,7 +459,7 @@ function update(lake::NaturalLake, i, inflow, doy, timestepsecs)
 
         if lake.outflowfunc[i] == 1
             outflow =
-                interpolate_linear(lake.waterlevel[i], lake.hq[i][!, 1], lake.hq[i][!, col])
+                interpolate_linear(lake.waterlevel[i], lake.hq[i].H, lake.hq[i].Q[:, col])
         else
             if diff_wl >= 0.0
                 outflow = max(


Test patch:
diff --git a/test/reservoir_lake.jl b/test/reservoir_lake.jl
index be831a555..864248d8a 100644
--- a/test/reservoir_lake.jl
+++ b/test/reservoir_lake.jl
@@ -37,9 +37,9 @@ end
         inflow = [0.0],
         b = [0.22],
         e = [2.0],
-        sh = [DataFrame()],
-        hq = [DataFrame()],
-        storage = Wflow.initialize_storage([1], [180510409.0], [18.5], [DataFrame()]),
+        sh = [missing],
+        hq = [missing],
+        storage = Wflow.initialize_storage([1], [180510409.0], [18.5], [missing]),
         waterlevel = [18.5],
         precipitation = [20.0],
         evaporation = [3.2],
@@ -57,10 +57,13 @@ end
 
 datadir = joinpath(@__DIR__, "data")
 sh = [
-    CSV.read(joinpath(datadir, "lake_sh_1.csv"), DataFrame, type = Float64),
-    CSV.read(joinpath(datadir, "lake_sh_2.csv"), DataFrame, type = Float64),
+    Wflow.read_sh_csv(joinpath(datadir, "lake_sh_1.csv")),
+    Wflow.read_sh_csv(joinpath(datadir, "lake_sh_2.csv")),
 ]
 @testset "linked lakes (HBV)" begin
+    @test keys(sh[1]) == (:H, :S)
+    @test typeof(values(sh[1])) == Tuple{Vector{Float64}, Vector{Float64}}
+
     lake = Wflow.NaturalLake{Float64}(
         lowerlake_ind = [2, 0],
         area = [472461536.0, 60851088.0],
@@ -73,8 +76,8 @@ sh = [
         e = [1.5, 1.5],
         sh = sh,
         hq = [
-            DataFrame(),
-            CSV.read(joinpath(datadir, "lake_hq_2.csv"), DataFrame, type = Float64),
+            missing,
+            Wflow.read_hq_csv(joinpath(datadir, "lake_hq_2.csv")),
         ],
         waterlevel = [395.03027, 394.87833],
         precipitation = [10.0, 10.0],
diff --git a/test/run_hbv.jl b/test/run_hbv.jl
index 91a969a05..ec1496c55 100644
--- a/test/run_hbv.jl
+++ b/test/run_hbv.jl
@@ -10,8 +10,7 @@ model = Wflow.update(model)
 # test if the first timestep was written to the CSV file
 flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
 @testset "CSV output" begin
-    csv = CSV.File(model.writer.csv_path)
-    row = csv[1]
+    row = csv_first_row(model.writer.csv_path)
 
     @test row.time == DateTime("2000-01-01T00:00:00")
     @test row.Q ≈ 517.9761666607884
diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index 7aa282473..592318ac9 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -11,8 +11,7 @@ model = Wflow.update(model)
 # test if the first timestep was written to the CSV file
 flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
 @testset "CSV output" begin
-    csv = CSV.File(model.writer.csv_path)
-    row = csv[1]
+    row = csv_first_row(model.writer.csv_path)
 
     @test row.time == DateTime("2000-01-01T00:00:00")
     @test row.Q ≈ 7.815587720999557
diff --git a/test/run_sbm_gwf.jl b/test/run_sbm_gwf.jl
index e639742bd..b3b65cd21 100644
--- a/test/run_sbm_gwf.jl
+++ b/test/run_sbm_gwf.jl
@@ -10,8 +10,7 @@ model = Wflow.update_sbm_gwf(model)
 # test if the first timestep was written to the CSV file
 flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
 @testset "CSV output" begin
-    csv = CSV.File(model.writer.csv_path)
-    row = csv[1]
+    row = csv_first_row(model.writer.csv_path)
 
     @test row.time == DateTime("2000-06-01T00:00:00")
     @test row.Q_av ≈ 0.016260519762890745
diff --git a/test/runtests.jl b/test/runtests.jl
index ba441eac4..8b3a0276a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,12 +7,11 @@ using Statistics
 using Test
 using UnPack
 using Wflow
-using CSV
-using DataFrames
 using BenchmarkTools
 using Base.MathConstants: eulergamma
 using BasicModelInterface
 import Polynomials
+using DelimitedFiles
 
 const BMI = BasicModelInterface
 
diff --git a/test/testing_utils.jl b/test/testing_utils.jl
index 01e47d457..0f27b7c7f 100644
--- a/test/testing_utils.jl
+++ b/test/testing_utils.jl
@@ -136,3 +136,23 @@ function expint(n::Integer, z)
         return Eᵢ
     end
 end
+
+"Return the first row of a Wflow output CSV file as a NamedTuple"
+function csv_first_row(path)
+    # silly function to avoid needing CSV.jl as a test dependency
+    header, dataline = open(path) do io
+        header = readline(io)
+        dataline = readline(io)
+        (header, dataline)
+    end
+
+    names = Tuple(Symbol.(split(header, ',')))
+    ncol = length(names)
+    # this assumes the first column is a time, the rest a float
+    types = Tuple{DateTime, fill(Float64, ncol-1)...}
+
+    parts = split(dataline, ',')
+    values = parse.(Float64, parts[2:end])
+    row = NamedTuple{names, types}((DateTime(parts[1]), values...))
+    return row
+end