Repository: Deltares/Wflow.jl
Pull request #128 (2021-05-04 20:52:43)

use a single internal dimension order with increasing axes
Previously we accepted x and y dimensions in either order, and adapted to that internally. However, some issues arose with other variants present in and between the netCDF inputs that made it clear that we need to support more variants. Some benchmarking showed that moving away from in place netCDF loading hardly affected performance. So to be both more flexible, and reduce some internal code complexity, we treat the data at entry to always follow the dimension order (x, y, layer, time), and to have increasing axes for each dimension.


Patch:
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index a70173c56..3e15391dd 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -13,8 +13,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
-          - '1'
+          - '1.6'
         os:
           - ubuntu-latest
           - macOS-latest
diff --git a/Project.toml b/Project.toml
index ab895503d..520ac3b67 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Wflow"
 uuid = "d48b7d99-76e7-47ae-b1d5-ff0c1cf9a818"
 authors = ["Deltares and contributors"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 BasicModelInterface = "59605e27-edc0-445a-b93d-c09a3a50b330"
@@ -30,14 +30,15 @@ ProgressLogging = "0.1"
 StaticArrays = "0.12, 1.0"
 TOML = "1.0"
 UnPack = "1"
-julia = "1.5"
+julia = "1.6"
 
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "Logging", "Polynomials", "Random", "Test"]
+test = ["BenchmarkTools", "Downloads", "Logging", "Polynomials", "Random", "Test"]
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 8827386da..f4d833af3 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -7,19 +7,25 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## Unreleased
 
+## v0.3.0 - 2021-05-10
+
 ### Changed
 - New deposition process for coarse sediment in the reservoirs with a new parameter 
   `restrapefficiency` in the sediment model.
 - New variables added to the `LandSediment` and `RiverSediment` structs in order to save 
   more output from the sediment model.
 - Added variables `volume` and `inwater` to `SurfaceFlow` struct, this is convenient for the
-  coupling with the water quality model Delwaq. 
+  coupling with the water quality model Delwaq.
+- River water level (`h`) and discharge (`q`) forced directly into the `RiverSediment`
+  struct (instead of using the `OverlandFlowSediment` struct first). 
+- Require Julia 1.6 or later.
 
 ### Added
 - Modify model parameters and forcing through the TOML file (see [Modify parameters](@ref)).
 - Run wflow\_sbm (SBM + kinematic wave) in two parts (until recharge and after subsurface
   flow) from BMI, including the option to switch off the lateral subsurface component of
   wflow\_sbm.
+- Support more netCDF dimension and axis order variants.
 
 ### Fixed
 - Corrected a bug in sediment deposition in the river (case when incoming sediment load is
diff --git a/docs/src/config.md b/docs/src/config.md
index 7fd15df2a..97d95ed73 100644
--- a/docs/src/config.md
+++ b/docs/src/config.md
@@ -236,7 +236,7 @@ parameter = "vertical.temperature"
 location = "temp_byindex"
 name = "temp_index"
 index.x = 100
-index.y = 50
+index.y = 264
 parameter = "vertical.temperature"
 ```
 
@@ -287,7 +287,7 @@ parameter = "vertical.temperature"
 [[csv.column]]
 header = "temp_byindex"
 index.x = 100
-index.y = 50
+index.y = 264
 parameter = "vertical.temperature"
 
 [[csv.column]]
diff --git a/src/Wflow.jl b/src/Wflow.jl
index ff9eb68bb..071f77bf2 100644
--- a/src/Wflow.jl
+++ b/src/Wflow.jl
@@ -36,7 +36,7 @@ function Clock(config)
 end
 
 function Clock(config, reader)
-    nctimes = ncread(reader.dataset, "time")
+    nctimes = reader.dataset["time"][:]
     # if the config file does not have a start or endtime, folow the NetCDF times
     # and add them to the config
     # if the timestep is not given, use the difference between NetCDF time 1 and 2
diff --git a/src/bmi.jl b/src/bmi.jl
index 890d1e78d..f127a3658 100644
--- a/src/bmi.jl
+++ b/src/bmi.jl
@@ -257,7 +257,7 @@ function BMI.get_grid_x(model::Model, grid::Int)
     @unpack dataset = reader
     sel = active_indices(model.network, symbols(grids[grid]))
     inds = [sel[i][1] for i in eachindex(sel)]
-    x_nc = "x" in keys(dataset.dim) ? ncread(dataset, "x") : ncread(dataset, "lon")
+    x_nc = read_x_axis(dataset)
     return x_nc[inds]
 end
 
@@ -266,7 +266,7 @@ function BMI.get_grid_y(model::Model, grid::Int)
     @unpack dataset = reader
     sel = active_indices(model.network, symbols(grids[grid]))
     inds = [sel[i][2] for i in eachindex(sel)]
-    y_nc = "y" in keys(dataset.dim) ? ncread(dataset, "y") : ncread(dataset, "lat")
+    y_nc = read_y_axis(dataset)
     return y_nc[inds]
 end
 
diff --git a/src/hbv_model.jl b/src/hbv_model.jl
index 2fc6bd85d..8aebdbf04 100644
--- a/src/hbv_model.jl
+++ b/src/hbv_model.jl
@@ -27,11 +27,6 @@ function initialize_hbv_model(config::Config)
     nc = NCDataset(static_path)
     dims = dimnames(nc[param(config, "input.subcatchment")])
 
-    # There is no need to permute the dimensions of the data, since the active indices are
-    # correctly calculated in both ways.
-    # The dimension order only needs to be known for interpreting the LDD directions
-    # and creating the coordinate maps.
-    dims_xy = dims[2] in ("y", "lat")
     subcatch_2d = ncread(nc, param(config, "input.subcatchment"); allow_missing = true)
     # indices based on catchment
     inds, rev_inds = active_indices(subcatch_2d, missing)
@@ -264,13 +259,9 @@ function initialize_hbv_model(config::Config)
     )
 
     # read x, y coordinates and calculate cell length [m]
-    y_nc = "y" in keys(nc.dim) ? ncread(nc, "y") : ncread(nc, "lat")
-    x_nc = "x" in keys(nc.dim) ? ncread(nc, "x") : ncread(nc, "lon")
-    if dims_xy
-        y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
-    else
-        y = repeat(y_nc, outer = (1, length(x_nc)))[inds]
-    end
+    y_nc = read_y_axis(nc)
+    x_nc = read_x_axis(nc)
+    y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
     cellength = abs(mean(diff(x_nc)))
 
     xl = fill(mv, n)
@@ -434,7 +425,6 @@ function initialize_hbv_model(config::Config)
         lake = nothing,
     )
 
-    pcr_dir = dims_xy ? permute_indices(pcrdir) : pcrdir
     graph = flowgraph(ldd, inds, pcr_dir)
 
     riverslope =
@@ -497,7 +487,7 @@ function initialize_hbv_model(config::Config)
         lake = isempty(lake) ? nothing : lake.reverse_indices,
     )
     writer =
-        prepare_writer(config, reader, modelmap, indices_reverse, x_nc, y_nc, dims_xy, nc)
+        prepare_writer(config, reader, modelmap, indices_reverse, x_nc, y_nc, nc)
     close(nc)
 
     # for each domain save the directed acyclic graph, the traversion order,
diff --git a/src/horizontal_process.jl b/src/horizontal_process.jl
index 1b93b2b2d..d724d5557 100644
--- a/src/horizontal_process.jl
+++ b/src/horizontal_process.jl
@@ -1,5 +1,5 @@
 "Convert a gridded drainage direction to a directed graph"
-function flowgraph(ldd::AbstractVector, inds::AbstractVector, pcrdir::AbstractVector)
+function flowgraph(ldd::AbstractVector, inds::AbstractVector, pcr_dir::AbstractVector)
     # prepare a directed graph to be filled
     n = length(inds)
     graph = DiGraph(n)
@@ -9,7 +9,7 @@ function flowgraph(ldd::AbstractVector, inds::AbstractVector, pcrdir::AbstractVe
         ldd_val = ldd[from_node]
         # skip pits to prevent cycles
         ldd_val == 5 && continue
-        to_index = from_index + pcrdir[ldd_val]
+        to_index = from_index + pcr_dir[ldd_val]
         # find the node id of the downstream cell
         to_node = searchsortedfirst(inds, to_index)
         add_edge!(graph, from_node, to_node)
diff --git a/src/io.jl b/src/io.jl
index 839333d10..d70b65e34 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -72,7 +72,7 @@ Base.iterate(config::Config, state) = iterate(Dict(config), state)
 "Extract NetCDF variable name `ncname` from `var` (type `String` or `Config`). If `var` has 
 type `Config`, either `scale` and `offset` are expected (with `ncname`) or a `value` (uniform
 value), these are stored as part of `NamedTuple` `modifier`"
-function ncvar_name_modifier(var)
+function ncvar_name_modifier(var; verbose=true)
     ncname = nothing
     modifier = (scale = 1.0, offset = 0.0, value = nothing)
     if isa(var, Config)
@@ -83,7 +83,9 @@ function ncvar_name_modifier(var)
             scale = param(var, "scale", 1.0)
             offset = param(var, "offset", 0.0)
             modifier = (scale = scale, offset = offset, value = nothing)
-            @info "NetCDF parameter $ncname is modified with scale $scale and offset $offset"
+            if verbose
+                @info "NetCDF parameter $ncname is modified with scale $scale and offset $offset"
+            end
         elseif haskey(var, "value")
             modifier = (scale = 1.0, offset = 0.0, value = param(var, "value"))
         end
@@ -94,9 +96,9 @@ function ncvar_name_modifier(var)
 end
 
 "Extract a NetCDF variable at a given time"
-function get_at!(
-    buffer,
-    var::NCDatasets.CFVariable,
+function get_at(
+    ds::NCDataset,
+    varname::AbstractString,
     times::AbstractVector{<:TimeType},
     t::TimeType,
 )
@@ -104,53 +106,18 @@ function get_at!(
     i = findfirst(>=(t), times)
     t < first(times) && throw(DomainError("time $t before dataset begin $(first(times))"))
     i === nothing && throw(DomainError("time $t after dataset end $(last(times))"))
-    return get_at!(buffer, var, i)
+    return get_at(ds, varname, i)
 end
 
-function get_at!(buffer, var::NCDatasets.CFVariable, i)
-    # assumes the dataset has 12 time steps, from January to December
-    dim = findfirst(==("time"), NCDatasets.dimnames(var))
-    # load in place, using a lower level NCDatasets function
-    if ndims(var) == 3
-        if dim == 1
-            NCDatasets.load!(var.var, buffer, i, :, :)
-        elseif dim == 3
-            NCDatasets.load!(var.var, buffer, :, :, i)
-        else
-            error("Time dimension expected at position 1 or 3 (number of dimensions is 3)")
-        end
-    elseif ndims(var) == 4
-        if dim == 1
-            len = size(var, 2)
-            if len == 1
-                NCDatasets.load!(var.var, buffer, i, 1, :, :)
-            else
-                dimname = NCDatasets.dimnames(var)[2]
-                error("Unsupported dimension $dimname of size $len")
-            end
-        elseif dim == 4
-            len = size(var, 3)
-            if len == 1
-                NCDatasets.load!(var.var, buffer, :, :, 1, i)
-            else
-                dimname = NCDatasets.dimnames(var)[3]
-                error("Unsupported dimension $dimname of size $len")
-            end
-        else
-            error("Time dimension expected at position 1 or 4 (number of dimensions is 4)")
-        end
-    else
-        error("Expected a 3 or 4 dimensional variable in the NetCDF")
-    end
-    return buffer
+function get_at(ds::NCDataset, varname::AbstractString, i)
+    return read_standardized(ds, varname, (x = :, y = :, time=i))
 end
 
 "Get dynamic NetCDF input for the given time"
 function update_forcing!(model)
     @unpack vertical, clock, reader, network, config = model
-    @unpack dataset, forcing_parameters, buffer = reader
-    sel = network.land.indices
-    nctimes = ncread(dataset, "time")
+    @unpack dataset, forcing_parameters = reader
+    nctimes = dataset["time"][:]
 
     do_reservoirs = get(config.model, "reservoirs", false)::Bool
     do_lakes = get(config.model, "lakes", false)::Bool
@@ -179,10 +146,10 @@ function update_forcing!(model)
     # load from NetCDF into the model according to the mapping
     for (par, ncvar) in forcing_parameters
         time = convert(eltype(nctimes), clock.time)
-        buffer = get_at!(buffer, dataset[ncvar.name], nctimes, time)
+        data = get_at(dataset, ncvar.name, nctimes, time)
 
         if ncvar.scale != 1.0 || ncvar.offset != 0.0
-            buffer .= buffer .* ncvar.scale .+ ncvar.offset
+            data .= data .* ncvar.scale .+ ncvar.offset
         end
 
         # calculate the mean precipitation and evaporation over the lakes and reservoirs
@@ -191,22 +158,29 @@ function update_forcing!(model)
         if par in mover_params
             if do_reservoirs
                 for (i, sel_reservoir) in enumerate(sel_reservoirs)
-                    avg = mean(buffer[sel_reservoir])
-                    buffer[sel_reservoir] .= 0
+                    avg = mean(data[sel_reservoir])
+                    data[sel_reservoir] .= 0
                     param_res[par][i] = avg
                 end
             end
             if do_lakes
                 for (i, sel_lake) in enumerate(sel_lakes)
-                    avg = mean(buffer[sel_lake])
-                    buffer[sel_lake] .= 0
+                    avg = mean(data[sel_lake])
+                    data[sel_lake] .= 0
                     param_lake[par][i] = avg
                 end
             end
         end
 
         param_vector = param(model, par)
-        param_vector .= buffer[sel]
+        sel = active_indices(network, par)
+        data_sel = data[sel]
+        if any(ismissing, data_sel)
+            print(par)
+            msg = "Forcing data has missing values on active model cells for $(ncvar.name)"
+            throw(ArgumentError(msg))
+        end
+        param_vector .= data_sel
     end
 
     return model
@@ -230,7 +204,7 @@ monthday_passed(curr, avail) = (curr[1] >= avail[1]) && (curr[2] >= avail[2])
 "Get cyclic NetCDF input for the given time"
 function update_cyclic!(model)
     @unpack vertical, clock, reader, network, config = model
-    @unpack cyclic_dataset, cyclic_times, cyclic_parameters, buffer = reader
+    @unpack cyclic_dataset, cyclic_times, cyclic_parameters = reader
     sel = network.land.indices
 
     # pick up the data that is valid for the past 24 hours
@@ -244,9 +218,9 @@ function update_cyclic!(model)
 
         # load from NetCDF into the model according to the mapping
         for (par, ncvar) in cyclic_parameters
-            buffer = get_at!(buffer, cyclic_dataset[ncvar.name], i)
+            data = get_at(cyclic_dataset, ncvar.name, i)
             param_vector = param(model, par)
-            param_vector .= buffer[sel]
+            param_vector .= data[sel]
             if ncvar.scale != 1.0 || ncvar.offset != 0.0
                 param_vector .= param_vector .* ncvar.scale .+ ncvar.offset
             end
@@ -452,24 +426,12 @@ function add_time(ds, time)
     return i
 end
 
-function checkdims(dims)
-    # TODO check if the x y ordering is equal to the staticmaps NetCDF
-    # with ensembles (e.g. Delft-FEWS) dims == 4
-    @assert length(dims) == 3 || length(dims) == 4
-    @assert "time" in dims
-    @assert ("x" in dims) || ("lon" in dims)
-    @assert ("y" in dims) || ("lat" in dims)
-    @assert dims[2] != "time"
-    return dims
-end
-
-struct NCReader{T}
+struct NCReader
     dataset::NCDataset
     cyclic_dataset::Union{NCDataset,Nothing}
     cyclic_times::Vector{Tuple{Int,Int}}
     forcing_parameters::Dict{Tuple{Symbol,Vararg{Symbol}},NamedTuple}
     cyclic_parameters::Dict{Tuple{Symbol,Vararg{Symbol}},NamedTuple}
-    buffer::Matrix{T}
 end
 
 struct Writer
@@ -490,29 +452,10 @@ end
 
 function prepare_reader(path, cyclic_path, config)
     dataset = NCDataset(path)
-    ncvar1, _ = ncvar_name_modifier(param(config, "input." * first(config.input.forcing)))
+    # set verbose to false to avoid logging modifications twice
+    ncvar1, _ = ncvar_name_modifier(param(config, "input." * first(config.input.forcing)); verbose=false)
     var = dataset[ncvar1].var
 
-    scale_factor = get(var.attrib, "scale_factor", nothing)
-    add_offset = get(var.attrib, "add_offset", nothing)
-    # TODO support scale_factor and add_offset with in place loading
-    # TODO check other forcing parameters as well
-    if !(isnothing(scale_factor) || isone(scale_factor))
-        error(
-            "scale_factor in NetCDF forcing not supported, found $scale_factor in $ncvar1",
-        )
-    end
-    if !(isnothing(add_offset) || iszero(add_offset))
-        error("add_offset in NetCDF forcing not supported, found $add_offset in $ncvar1")
-    end
-
-    T = eltype(var)
-    dims = dimnames(var)
-    checkdims(dims)
-    timelast = last(dims) == "time"
-    lateral_size = timelast ? size(var)[1:2] : size(var)[2:3]
-    buffer = zeros(T, lateral_size)
-
     # check for cyclic parameters
     do_cyclic = haskey(config.input, "cyclic")
 
@@ -565,7 +508,6 @@ function prepare_reader(path, cyclic_path, config)
         cyclic_times,
         forcing_parameters,
         cyclic_parameters,
-        buffer,
     )
 end
 
@@ -728,7 +670,6 @@ function prepare_writer(
     rev_inds,
     x_nc,
     y_nc,
-    dims_xy,
     nc_static;
     maxlayers = nothing,
 )
@@ -799,7 +740,7 @@ function prepare_writer(
         for var in config.netcdf.variable
             parameter = var["parameter"]
             reducer_func =
-                get_reducer_func(var, rev_inds, x_nc, y_nc, dims_xy, config, nc_static)
+                get_reducer_func(var, rev_inds, x_nc, y_nc, config, nc_static)
             push!(nc_scalar, (parameter = parameter, reducer = reducer_func))
         end
     else
@@ -826,7 +767,7 @@ function prepare_writer(
         for col in config.csv.column
             parameter = col["parameter"]
             reducer_func =
-                get_reducer_func(col, rev_inds, x_nc, y_nc, dims_xy, config, nc_static)
+                get_reducer_func(col, rev_inds, x_nc, y_nc, config, nc_static)
             push!(csv_cols, (parameter = parameter, reducer = reducer_func))
         end
     else
@@ -871,10 +812,10 @@ end
 "Write a new timestep with grid data to a NetCDF file"
 function write_netcdf_timestep(model, dataset, parameters)
     @unpack vertical, clock, reader, network = model
-    @unpack buffer = reader
 
     time_index = add_time(dataset, clock.time)
 
+    buffer = zeros(Union{Float, Missing}, size(model.network.land.reverse_indices))
     for (key, val) in parameters
         @unpack par, vector = val
         sel = active_indices(network, par)
@@ -882,14 +823,14 @@ function write_netcdf_timestep(model, dataset, parameters)
         elemtype = eltype(vector)
         if elemtype <: AbstractFloat
             # ensure no other information is written
-            fill!(buffer, NaN)
+            fill!(buffer, missing)
             buffer[sel] .= vector
             dataset[key][:, :, time_index] = buffer
         elseif elemtype <: SVector
             nlayer = length(first(vector))
             for i = 1:nlayer
                 # ensure no other information is written
-                fill!(buffer, NaN)
+                fill!(buffer, missing)
                 buffer[sel] .= getindex.(vector, i)
                 dataset[key][:, :, i, time_index] = buffer
             end
@@ -1001,7 +942,7 @@ function reducerfunction(reducer::AbstractString)
 end
 
 "Get a reducer function based on output settings for scalar data defined in a dictionary"
-function reducer(col, rev_inds, x_nc, y_nc, dims_xy, config, dataset)
+function reducer(col, rev_inds, x_nc, y_nc, config, dataset)
     if haskey(col, "map")
         # assumes the parameter in "map" has a 2D input map, with
         # integers indicating the points or zones that are to be aggregated
@@ -1041,13 +982,10 @@ function reducer(col, rev_inds, x_nc, y_nc, dims_xy, config, dataset)
             return x -> getindex(x, index)
         elseif index isa Dict
             # index into the 2D input/output arrays
-            # the first always corresponds to the y dimension, then the x dimension
+            # the first always corresponds to the x dimension, then the y dimension
             # this is 1-based
-            i = index["y"]::Int
-            j = index["x"]::Int
-            if dims_xy
-                i, j = j, i
-            end
+            i = index["x"]::Int
+            j = index["y"]::Int
             ind = rev_inds[i, j]
             iszero(ind) && error("inactive loc specified for output")
             return A -> getindex(A, ind)
@@ -1060,7 +998,7 @@ function reducer(col, rev_inds, x_nc, y_nc, dims_xy, config, dataset)
         # find the closest cell center index
         _, iy = findmin(abs.(y_nc .- y))
         _, ix = findmin(abs.(x_nc .- x))
-        I = dims_xy ? CartesianIndex(ix, iy) : CartesianIndex(iy, ix)
+        I = CartesianIndex(ix, iy)
         i = rev_inds[I]
         iszero(i) && error("inactive coordinate specified for output")
         return A -> getindex(A, i)
@@ -1142,3 +1080,218 @@ end
 # these represent the type of the rating curve and specific storage data
 const SH = NamedTuple{(:H, :S),Tuple{Vector{Float},Vector{Float}}}
 const HQ = NamedTuple{(:H, :Q),Tuple{Vector{Float},Matrix{Float}}}
+
+is_increasing(v) = last(v) > first(v)
+
+"""
+    nc_dim_name(dims::Vector{Symbol}, name::Symbol)
+    nc_dim_name(ds::NCDataset, name::Symbol)
+
+Given a NetCDF dataset or list of dimensions, and an internal dimension name, return the
+corresponding NetCDF dimension name. Certain common alternatives are supported, e.g. :lon or
+:longitude instead of :x.
+"""
+function nc_dim_name(dims::Vector{Symbol}, name::Symbol)
+    # direct naming
+    name in dims && return name
+    # list of common alternative names
+    if name == :x
+        for candidate in (:lon, :longitude)
+            candidate in dims && return candidate
+        end
+        error("No x dimension candidate found")
+    elseif name == :y
+        for candidate in (:lat, :latitude)
+            candidate in dims && return candidate
+        end
+        error("No y dimension candidate found")
+    else
+        error("Unknown dimension $name")
+    end
+end
+
+nc_dim_name(ds::NCDataset, name::Symbol) = nc_dim_name(Symbol.(keys(ds.dim)), name)
+
+"""
+    nc_dim(ds::NCDataset, name::Symbol)
+
+Return the dimension coordinate, based on the internal name (:x, :y, :layer, :time),
+which will map to the correct NetCDF name using `nc_dim_name`.
+"""
+nc_dim(ds::NCDataset, name) = ds[nc_dim_name(ds, name)]
+
+
+"""
+    internal_dim_name(name::Symbol)
+
+Given a NetCDF dimension name string, return the corresponding internal dimension name.
+"""
+function internal_dim_name(name::Symbol)
+    if name in (:x, :lon, :longitude)
+        return :x
+    elseif name in (:y, :lat, :latitude)
+        return :y
+    elseif name in (:time, :layer)
+        return name
+    else
+        error("Unknown dimension $name")
+    end
+end
+
+"""
+    read_dims(A::NCDatasets.CFVariable, dim_sel)
+
+Return the data of a NetCDF data variable as an Array. Only dimensions in `dim_sel`, a
+NamedTuple like (x=:, y=:, time=1). Other dimensions that may be present need to be size 1,
+otherwise an error is thrown.
+
+`dim_sel` keys should be the internal dimension names; :x, :y, :time, :layer.
+"""
+function read_dims(A::NCDatasets.CFVariable, dim_sel::NamedTuple)
+    dimsizes = dimsize(A)
+    indexer = []
+    data_dim_order = Symbol[]
+    # need to iterate in this order, to get the indices in order
+    for (dim_name, dim_size) in pairs(dimsizes)
+        dim = internal_dim_name(dim_name)
+        # each dim should either be in dim_sel or be size 1
+        if dim in keys(dim_sel)
+            idx = dim_sel[dim]
+            push!(indexer, idx)
+            # Would be nice to generalize this to anything that drops the dimension
+            # when used to index. Not sure how we could support `begin` or `end` here.
+            if !(idx isa Int)
+                push!(data_dim_order, dim)
+            end
+        elseif dim_size == 1
+            push!(indexer, 1)
+        else
+            throw(ArgumentError("""NetCDF dimension $dim_name has length $dim_size.
+                Only extra dimensions of length 1 are supported."""))
+        end
+    end
+    data = A[indexer...]
+    order = Tuple(data_dim_order)
+    @assert ndims(data) == length(order) "dimensions not properly recorded"
+    return data, order
+end
+
+"""
+    dim_directions(ds, dim_names)
+
+Given a NCDataset and a list of internal dimension names, return a NamedTuple, which maps
+the dimension name to `true` if it is increasing, or `false` otherwise.
+
+For the layer dimension, we allow coordinate arrays to be missing, in which case we
+consider it increasing, going from the top layer (1) to deeper layers. This is to keep
+accepting data that we have accepted before.
+"""
+function dim_directions(ds::NCDataset, dim_names)
+    pairs = Pair{Symbol,Bool}[]
+    for d in dim_names
+        if d == :layer && !(haskey(ds, "layer"))
+            inc = true
+        else
+            inc = is_increasing(nc_dim(ds, d))
+        end
+        push!(pairs, d => inc)
+    end
+    return NamedTuple(pairs)
+end
+
+"""
+    permute_data(data, dim_names)
+
+Given an Array of data, and a list of its dimension names, return a permuted version
+of the data such that the dimension order will follow (:x, :y, :layer). No permutation
+is done if this is not needed.
+"""
+function permute_data(data, dim_names)
+    @assert ndims(data) == length(dim_names)
+    desired_order = (:x, :y, :layer)
+    @assert dim_names ⊆ desired_order
+    if length(dim_names) == 2
+        if first(dim_names) == :x
+            return data, dim_names
+        else
+            return permutedims(data), reverse(dim_names)
+        end
+    elseif length(dim_names) == 3
+        permutation = [findfirst(==(d), dim_names) for d in desired_order]
+        if permutation == (1, 2, 3)
+            return data, dim_names
+        else
+            return permutedims(data, permutation), dim_names[permutation]
+        end
+    else
+        error("Unsupported number of dimensions")
+    end
+end
+
+"""
+    reverse_data!(data, dims_increasing)
+
+Reverse the data as needed, such that it is increasing along all dimensions.
+"""
+function reverse_data!(data, dims_increasing)
+    # for the reverse call it is important that the dims_increasing tuple is ordered in the
+    # desired internal ordering, just like the data is after permutation
+    if length(dims_increasing) == 2
+        dims_increasing_ordered = (x = dims_increasing.x, y = dims_increasing.y)
+    elseif length(dims_increasing) == 3
+        dims_increasing_ordered =
+            (x = dims_increasing.x, y = dims_increasing.y, layer = dims_increasing.layer)
+    else
+        error("Unsupported number of dimensions")
+    end
+    # the non increasing dimensions should be reversed, to make them all increasing
+    dims = Tuple(findall(.!values(dims_increasing_ordered)))
+    return reverse!(data; dims)
+end
+
+"""
+    read_standardized(ds::NCDataset, varname::AbstractString, dim_names)
+
+Read the dimensions listed in dim_names from a variable with name `varname` from a NetCDF
+dataset `ds`. `dim_sel` should be a NamedTuple like (x=:, y=:, time=1), which will return
+a 2 dimensional array with x and y axes, representing the first index in the time dimension.
+"""
+function read_standardized(ds::NCDataset, varname::AbstractString, dim_sel::NamedTuple)
+    data, data_dim_order = read_dims(ds[varname], dim_sel)
+    data, new_dim_order = permute_data(data, data_dim_order)
+    dims_increasing = dim_directions(ds, new_dim_order)
+    data = reverse_data!(data, dims_increasing)
+    return data
+end
+
+"""
+    read_x_axis(ds::NCDataset)
+
+Return the x coordinate Vector{Float64}, whether it is called x, lon or longitude.
+Also sorts the vector to be increasing, to match `read_standardized`.
+"""
+function read_x_axis(ds::NCDataset)
+    candidates = ("x", "lon", "longitude")
+    for candidate in candidates
+        if haskey(ds, candidate)
+            return sort!(Float64.(ds[candidate][:]))
+        end
+    end
+    error("no x axis found in $(path(ds))")
+end
+
+"""
+    read_y_axis(ds::NCDataset)
+
+Return the y coordinate Vector{Float64}, whether it is called y, lat or latitude.
+Also sorts the vector to be increasing, to match `read_standardized`.
+"""
+function read_y_axis(ds::NCDataset)
+    candidates = ("y", "lat", "latitude")
+    for candidate in candidates
+        if haskey(ds, candidate)
+            return sort!(Float64.(ds[candidate][:]))
+        end
+    end
+    error("no y axis found in $(path(ds))")
+end
diff --git a/src/sbm.jl b/src/sbm.jl
index a823f7b99..169f4581d 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -436,7 +436,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         sel = inds,
         defaults = 10.0,
         type = Float,
-        dimname = "layer",
+        dimname = :layer,
     )
     if size(c, 1) != maxlayers
         parname = param(config, "input.vertical.c")
@@ -449,7 +449,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         sel = inds,
         defaults = 1.0,
         type = Float,
-        dimname = "layer",
+        dimname = :layer,
     )
     if size(kvfrac, 1) != maxlayers
         parname = param(config, "input.vertical.khfrac")
diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
index 9a871385f..77706ddcf 100644
--- a/src/sbm_gwf_model.jl
+++ b/src/sbm_gwf_model.jl
@@ -36,13 +36,6 @@ function initialize_sbm_gwf_model(config::Config)
     kw_land_tstep = get(config.model, "kw_land_tstep", 0)
 
     nc = NCDataset(static_path)
-    dims = dimnames(nc[param(config, "input.subcatchment")])
-
-    # There is no need to permute the dimensions of the data, since the active indices are
-    # correctly calculated in both ways.
-    # The dimension order only needs to be known for interpreting the LDD directions
-    # and creating the coordinate maps.
-    dims_xy = dims[2] in ("y", "lat")
 
     subcatch_2d = ncread(nc, param(config, "input.subcatchment"); allow_missing = true)
     # indices based on catchment
@@ -61,13 +54,9 @@ function initialize_sbm_gwf_model(config::Config)
 
     altitude = ncread(nc, param(config, "input.altitude"); sel = inds, type = Float)
     # read x, y coordinates and calculate cell length [m]
-    y_nc = "y" in keys(nc.dim) ? ncread(nc, "y") : ncread(nc, "lat")
-    x_nc = "x" in keys(nc.dim) ? ncread(nc, "x") : ncread(nc, "lon")
-    if dims_xy
-        y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
-    else
-        y = repeat(y_nc, outer = (1, length(x_nc)))[inds]
-    end
+    y_nc = read_y_axis(nc)
+    x_nc = read_x_axis(nc)
+    y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
     cellength = abs(mean(diff(x_nc)))
 
     xl = fill(mv, n)
@@ -108,6 +97,7 @@ function initialize_sbm_gwf_model(config::Config)
     βₗ = ncread(nc, param(config, "input.lateral.land.slope"); sel = inds, type = Float)
     clamp!(βₗ, 0.00001, Inf)
     ldd_2d = ncread(nc, param(config, "input.ldd"); allow_missing = true)
+
     ldd = ldd_2d[inds]
     dl = fill(mv, n)
     dw = fill(mv, n)
@@ -158,7 +148,6 @@ function initialize_sbm_gwf_model(config::Config)
         lake = nothing,
     )
 
-    pcr_dir = dims_xy ? permute_indices(pcrdir) : pcrdir
     graph = flowgraph(ldd, inds, pcr_dir)
 
     # river flow (kinematic wave)
@@ -364,7 +353,6 @@ function initialize_sbm_gwf_model(config::Config)
         indices_reverse,
         x_nc,
         y_nc,
-        dims_xy,
         nc,
         maxlayers = sbm.maxlayers,
     )
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index 6c8a461b7..6362579b1 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -28,12 +28,6 @@ function initialize_sbm_model(config::Config)
     nc = NCDataset(static_path)
     dims = dimnames(nc[param(config, "input.subcatchment")])
 
-    # There is no need to permute the dimensions of the data, since the active indices are
-    # correctly calculated in both ways.
-    # The dimension order only needs to be known for interpreting the LDD directions
-    # and creating the coordinate maps.
-    dims_xy = dims[2] in ("y", "lat")
-
     subcatch_2d = ncread(nc, param(config, "input.subcatchment"); allow_missing = true)
     # indices based on catchment
     inds, rev_inds = active_indices(subcatch_2d, missing)
@@ -50,13 +44,9 @@ function initialize_sbm_model(config::Config)
     riverlength = riverlength_2d[inds]
 
     # read x, y coordinates and calculate cell length [m]
-    y_nc = "y" in keys(nc.dim) ? ncread(nc, "y") : ncread(nc, "lat")
-    x_nc = "x" in keys(nc.dim) ? ncread(nc, "x") : ncread(nc, "lon")
-    if dims_xy
-        y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
-    else
-        y = repeat(y_nc, outer = (1, length(x_nc)))[inds]
-    end
+    y_nc = read_y_axis(nc)
+    x_nc = read_x_axis(nc)
+    y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
     cellength = abs(mean(diff(x_nc)))
 
     xl = fill(mv, n)
@@ -197,7 +187,6 @@ function initialize_sbm_model(config::Config)
         lake = nothing,
     )
 
-    pcr_dir = dims_xy ? permute_indices(pcrdir) : pcrdir
     graph = flowgraph(ldd, inds, pcr_dir)
 
 
@@ -267,7 +256,6 @@ function initialize_sbm_model(config::Config)
         indices_reverse,
         x_nc,
         y_nc,
-        dims_xy,
         nc,
         maxlayers = sbm.maxlayers,
     )
diff --git a/src/sediment.jl b/src/sediment.jl
index 8da72e716..56cbeba35 100644
--- a/src/sediment.jl
+++ b/src/sediment.jl
@@ -667,9 +667,6 @@ end
     n::Int | "-"
     # Filter with river cells
     rivcell::Vector{T} | "-"
-    # Forcing for the river
-    h_riv::Vector{T} | "m"
-    q_riv::Vector{T} | "m3 s-1"
     # Total eroded soil [ton]
     soilloss::Vector{T} | "t"
     # Eroded soil per particle class [ton]
diff --git a/src/sediment_model.jl b/src/sediment_model.jl
index f8e17e69c..9ba50d6e5 100644
--- a/src/sediment_model.jl
+++ b/src/sediment_model.jl
@@ -23,12 +23,6 @@ function initialize_sediment_model(config::Config)
     nc = NCDataset(static_path)
     dims = dimnames(nc[param(config, "input.subcatchment")])
 
-    # There is no need to permute the dimensions of the data, since the active indices are
-    # correctly calculated in both ways.
-    # The dimension order only needs to be known for interpreting the LDD directions
-    # and creating the coordinate maps.
-    dims_xy = dims[2] in ("y", "lat")
-
     subcatch_2d = ncread(nc, param(config, "input.subcatchment"); allow_missing = true)
     # indices based on catchment
     inds, rev_inds = active_indices(subcatch_2d, missing)
@@ -52,13 +46,9 @@ function initialize_sediment_model(config::Config)
     lake = ()
 
     # read x, y coordinates and calculate cell length [m]
-    y_nc = "y" in keys(nc.dim) ? ncread(nc, "y") : ncread(nc, "lat")
-    x_nc = "x" in keys(nc.dim) ? ncread(nc, "x") : ncread(nc, "lon")
-    if dims_xy
-        y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
-    else
-        y = repeat(y_nc, outer = (1, length(x_nc)))[inds]
-    end
+    y_nc = read_y_axis(nc)
+    x_nc = read_x_axis(nc)
+    y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
     cellength = abs(mean(diff(x_nc)))
 
     xl = fill(mv, n)
@@ -82,8 +72,6 @@ function initialize_sediment_model(config::Config)
     ols = OverlandFlowSediment{Float}(
         n = n,
         rivcell = rivcell,
-        h_riv = fill(mv, n),
-        q_riv = fill(mv, n),
         soilloss = fill(mv, n),
         erosclay = fill(mv, n),
         erossilt = fill(mv, n),
@@ -110,7 +98,6 @@ function initialize_sediment_model(config::Config)
         inlandlagg = fill(mv, n),
     )
 
-    pcr_dir = dims_xy ? permute_indices(pcrdir) : pcrdir
     graph = flowgraph(ldd, inds, pcr_dir)
 
     # River processes
@@ -136,7 +123,7 @@ function initialize_sediment_model(config::Config)
         lake = isempty(lake) ? nothing : lake.reverse_indices,
     )
     writer =
-        prepare_writer(config, reader, modelmap, indices_reverse, x_nc, y_nc, dims_xy, nc)
+        prepare_writer(config, reader, modelmap, indices_reverse, x_nc, y_nc, nc)
     close(nc)
 
     # for each domain save the directed acyclic graph, the traversion order,
@@ -210,12 +197,7 @@ function update(model::Model{N,L,V,R,W}) where {N,L,V<:LandSediment,R,W}
     do_river = get(config.model, "runrivermodel", false)::Bool
 
     if do_river
-        # Forcing come from lateral.land instead of netcdf directly
-        # Fix update_forcing in io to allow forcing for river (problem wih sel)
         inds_riv = network.index_river
-        lateral.river.h_riv .= lateral.land.h_riv[inds_riv]
-        lateral.river.q_riv .= lateral.land.q_riv[inds_riv]
-
         lateral.river.inlandclay .= lateral.land.inlandclay[inds_riv]
         lateral.river.inlandsilt .= lateral.land.inlandsilt[inds_riv]
         lateral.river.inlandsand .= lateral.land.inlandsand[inds_riv]
diff --git a/src/utils.jl b/src/utils.jl
index 20a639c77..3acb0c0ad 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,15 +1,15 @@
 
 "Map from PCRaster LDD value to a CartesianIndex"
-const pcrdir = [
-    CartesianIndex(1, -1),  # 1
-    CartesianIndex(1, 0),  # 2
-    CartesianIndex(1, 1),  # 3
-    CartesianIndex(0, -1),  # 4
+const pcr_dir = [
+    CartesianIndex(-1, -1),  # 1
+    CartesianIndex(0, -1),  # 2
+    CartesianIndex(1, -1),  # 3
+    CartesianIndex(-1, 0),  # 4
     CartesianIndex(0, 0),  # 5
-    CartesianIndex(0, 1),  # 6
-    CartesianIndex(-1, -1),  # 7
-    CartesianIndex(-1, 0),  # 8
-    CartesianIndex(-1, 1),  # 9
+    CartesianIndex(1, 0),  # 6
+    CartesianIndex(-1, 1),  # 7
+    CartesianIndex(0, 1),  # 8
+    CartesianIndex(1, 1),  # 9
 ]
 
 const mv = Float(NaN)
@@ -17,9 +17,6 @@ const mv = Float(NaN)
 # timestep that the parameter units are defined in
 const basetimestep = Second(Day(1))
 
-# in case the input data is permuted the other way
-permute_indices(inds) = [CartesianIndex(i[2], i[1]) for i in inds]
-
 "Set at indices pit values (default = 5) in a gridded local drainage direction vector"
 function set_pit_ldd(pits_2d, ldd, indices; pit = 5)
     pits = pits_2d[indices]
@@ -112,7 +109,8 @@ function set_states(instate_path, model, state_ncnames; type = nothing)
             dims = length(dimnames(ds[ncname]))
             # 4 dims, for example (x,y,layer,time) where dim layer is an SVector for soil layers
             if dims == 4
-                A = permutedims(ds[ncname][sel, :, 1])
+                A = read_standardized(ds, ncname, (x = :, y = :, layer = :, time=1))
+                A = permutedims(A[sel, :])
                 # note that this array is allowed to have missings, since not every vertical
                 # column is `maxlayers` layers deep
                 A = replace!(A, missing => NaN)
@@ -126,7 +124,8 @@ function set_states(instate_path, model, state_ncnames; type = nothing)
                 param(model, state) .= svectorscopy(A, Val{size(A)[1]}())
                 # 3 dims (x,y,time)
             elseif dims == 3
-                A = ds[ncname][sel, 1]
+                A = read_standardized(ds, ncname, (x = :, y = :, time=1))
+                A = A[sel]
                 A = nomissing(A)
                 # Convert to desired type if needed
                 if !isnothing(type)
@@ -174,6 +173,13 @@ function ncread(
     dimname = nothing,
 )
 
+    if isnothing(dimname)
+        dim_sel = (x = :, y = :)
+    else
+        @assert dimname == :layer
+        dim_sel = (x = :, y = :, layer = :)
+    end
+
     if isnothing(var)
         @assert !isnothing(defaults)
         if !isnothing(type)
@@ -182,7 +188,7 @@ function ncread(
         if isnothing(dimname)
             return Base.fill(defaults, length(sel))
         else
-            return Base.fill(defaults, (nc.dim[dimname], length(sel)))
+            return Base.fill(defaults, (nc.dim[String(dimname)], length(sel)))
         end
     end
 
@@ -191,17 +197,17 @@ function ncread(
     # NetCDF var is read directly.
     var, mod = ncvar_name_modifier(var)
     if mod.scale != 1.0 || mod.offset != 0.0
-        A = nc[var][:] .* mod.scale .+ mod.offset
+        A = read_standardized(nc, var, dim_sel) .* mod.scale .+ mod.offset
     elseif !isnothing(mod.value)
         if isnothing(dimname)
             return Base.fill(mod.value, length(sel))
         else
-            return Base.fill(mod.value, (nc.dim[dimname], length(sel)))
+            return Base.fill(mod.value, (nc.dim[String(dimname)], length(sel)))
         end
     else
         # Read the entire variable into memory, applying scale, offset and
         # set fill_values to missing.
-        A = nc[var][:]
+        A = read_standardized(nc, var, dim_sel)
     end
 
     # Take out only the active cells
@@ -209,12 +215,7 @@ function ncread(
         if isnothing(dimname)
             A = A[sel]
         else
-            dim = findfirst(==(dimname), dimnames(nc[var]))
-            if dim == 3
-                A = permutedims(A[sel, :])
-            elseif dim == 1
-                A = A[:, sel]
-            end
+            A = permutedims(A[sel, :])
         end
     end
 


Test patch:
diff --git a/test/bmi.jl b/test/bmi.jl
index a4a427e62..5e9c17814 100644
--- a/test/bmi.jl
+++ b/test/bmi.jl
@@ -50,7 +50,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
             @test BMI.get_current_time(model) == 9.467712e8
             @test mean(BMI.get_value(model, "vertical.zi")) ≈ 276.252648379751
             @test BMI.get_value_at_indices(model, "lateral.river.q", [1, 100, 5617]) ≈
-                  [0.0019494398840044726, 0.0552464309192915, 2.662842827356121]
+                  [0.9784867989663215, 1.6359646870951168, 0.027160430280557993]
             BMI.set_value(model, "vertical.zi", fill(300.0, 50070))
             @test mean(BMI.get_value(model, "vertical.zi")) == 300.0
             BMI.set_value_at_indices(model, "vertical.zi", [1], [250.0])
@@ -64,10 +64,10 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
             @test BMI.get_grid_size(model, 2) == 5655
             @test BMI.get_grid_size(model, 4) == 50070
             @test BMI.get_grid_shape(model, 4) == 50070
-            @test minimum(BMI.get_grid_x(model, 4)) ≈ 5.4291666666666725
-            @test maximum(BMI.get_grid_x(model, 4)) ≈ 7.845833333333333
-            @test BMI.get_grid_x(model, 0) ≈ [5.920833333333338, 5.7625000000000055]
-            @test BMI.get_grid_y(model, 0) ≈ [49.9125, 48.920834]
+            @test minimum(BMI.get_grid_x(model, 4)) ≈ 5.4291666666666725f0
+            @test maximum(BMI.get_grid_x(model, 4)) ≈ 7.845833333333333f0
+            @test BMI.get_grid_x(model, 0) ≈ [5.7625f0, 5.920834f0]
+            @test BMI.get_grid_y(model, 0) ≈ [48.920834f0, 49.9125f0]
             @test BMI.get_grid_node_count(model, 0) == 2
         end
 
@@ -89,10 +89,10 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
 
         @testset "recharge part of SBM" begin
             sbm = model.vertical
-            @test sbm.interception[1] ≈ 0.6469898223876953f0
-            @test sbm.ustorelayerdepth[1][1] ≈ 0.8873711853394088f0
-            @test sbm.snow[1] ≈ 0.6029989752244306f0
-            @test sbm.recharge[5] ≈ -0.011238714493537794f0
+            @test sbm.interception[1] ≈ 0.6329902410507202f0
+            @test sbm.ustorelayerdepth[1][1] ≈ 0.0f0
+            @test sbm.snow[1] ≈ 3.4530311597151853f0
+            @test sbm.recharge[5] ≈ 0.0f0
             @test sbm.zi[5] ≈ 300.0f0
         end
 
@@ -105,10 +105,10 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
         @testset "SBM after subsurface flow" begin
             sbm = model.vertical
             sub = model.lateral.subsurface
-            @test sbm.interception[1] ≈ 0.6469898223876953f0
-            @test sbm.ustorelayerdepth[1][1] ≈ 0.8873711853394088f0
-            @test sbm.snow[1] ≈ 0.6029989752244306f0
-            @test sbm.recharge[5] ≈ -0.011238714493537794f0
+            @test sbm.interception[1] ≈ 0.6329902410507202f0
+            @test sbm.ustorelayerdepth[1][1] ≈ 0.0f0
+            @test sbm.snow[1] ≈ 3.4530311597151853f0
+            @test sbm.recharge[5] ≈ 0.0f0
             @test sbm.zi[5] ≈ 250.0f0
             @test sub.zi[5] ≈ 250.0f0
             @test sub.exfiltwater[1] ≈ 0.01f0
diff --git a/test/hbv_config.toml b/test/hbv_config.toml
index a421b0763..a60e10198 100644
--- a/test/hbv_config.toml
+++ b/test/hbv_config.toml
@@ -122,7 +122,7 @@ parameter = "vertical.temperature"
 [[csv.column]]
 header = "temp_byindex"
 index.x = 88
-index.y = 22
+index.y = 95
 parameter = "vertical.temperature"
 
 [[csv.column]]
diff --git a/test/horizontal_process.jl b/test/horizontal_process.jl
index 7c584d8eb..91942d303 100644
--- a/test/horizontal_process.jl
+++ b/test/horizontal_process.jl
@@ -3,23 +3,25 @@ const ldd_mv = 255
 
 # read the staticmaps into memory
 nc = NCDataset(staticmaps_rhine_path)
-ldd_2d = Wflow.ncread(nc, "ldd")
+# helper function to get the axis order and directionality right
+read_right(nc, var) = reverse(permutedims(nc[var][:]); dims=2)
+ldd_2d = read_right(nc, "ldd")
 
 inds, _ = Wflow.active_indices(ldd_2d, ldd_mv)
 n = length(inds)
 
 # take out only the active cells
 ldd = ldd_2d[inds]
-slope = Wflow.ncread(nc, "slope"; sel = inds)
-N = Wflow.ncread(nc, "N"; sel = inds)
-Qold = Wflow.ncread(nc, "Qold"; sel = inds)
-Bw = Wflow.ncread(nc, "Bw"; sel = inds)
-waterlevel = Wflow.ncread(nc, "waterlevel"; sel = inds)
-DCL = Wflow.ncread(nc, "DCL"; sel = inds)
+slope = read_right(nc, "slope")[inds]
+N = read_right(nc, "N")[inds]
+Qold = read_right(nc, "Qold")[inds]
+Bw = read_right(nc, "Bw")[inds]
+waterlevel = read_right(nc, "waterlevel")[inds]
+DCL = read_right(nc, "DCL")[inds]
 close(nc)
 
 # create the directed acyclic graph from the drainage direction array
-graph = Wflow.flowgraph(ldd, inds, Wflow.pcrdir)
+graph = Wflow.flowgraph(ldd, inds, Wflow.pcr_dir)
 # a topological sort is used for visiting nodes in order from upstream to downstream
 toposort = topological_sort_by_dfs(graph)
 sink = toposort[end]
@@ -38,8 +40,8 @@ Q = Wflow.kin_wave!(Q, graph, toposort, Qold, q, α, β, DCL, Δt_sec)
 
 @testset "flow rate" begin
     @test sum(Q) ≈ 2.957806043289641e6
-    @test Q[toposort[1]] ≈ 0.00021453683608501235
-    @test Q[toposort[n-100]] ≈ 4054.9507466731234
+    @test Q[toposort[1]] ≈ 0.007260052312634069f0
+    @test Q[toposort[n-100]] ≈ 3945.762718338739f0
     @test Q[sink] ≈ 4131.101474418251
 end
 
diff --git a/test/io.jl b/test/io.jl
index 73e3d0962..af3c98153 100644
--- a/test/io.jl
+++ b/test/io.jl
@@ -24,15 +24,6 @@ config = Wflow.Config(tomlpath)
     @test collect(keys(config.output)) == ["lateral", "vertical", "path"]
 end
 
-@testset "checkdims" begin
-    @test_throws AssertionError Wflow.checkdims(("z", "lat", "time"))
-    @test_throws AssertionError Wflow.checkdims(("z", "lat"))
-    @test Wflow.checkdims(("lon", "lat", "time")) == ("lon", "lat", "time")
-    @test Wflow.checkdims(("time", "lon", "lat")) == ("time", "lon", "lat")
-    @test Wflow.checkdims(("lat", "lon", "time")) == ("lat", "lon", "time")
-    @test Wflow.checkdims(("time", "lat", "lon")) == ("time", "lat", "lon")
-end
-
 @testset "Clock constructor" begin
     config = Wflow.Config(tomlpath)
 
@@ -153,7 +144,7 @@ end
     @test Wflow.monthday_passed((2, 2), (1, 1))  # month and day later
     @test !Wflow.monthday_passed((2, 1), (2, 2))  # day before
     @test !Wflow.monthday_passed((1, 2), (2, 2))  # month before
-    @test !Wflow.monthday_passed((1, 1), (2, 2))  # day and month before       
+    @test !Wflow.monthday_passed((1, 1), (2, 2))  # day and month before
 end
 
 # test reading and setting of warm states (reinit=false)
@@ -177,21 +168,21 @@ model = Wflow.initialize_sbm_model(config)
 end
 
 @testset "warm states" begin
-    @test Wflow.param(model, "lateral.river.reservoir.volume")[2] ≈ 2.7393418e7
-    @test Wflow.param(model, "vertical.satwaterdepth")[41120] ≈ 201.51429748535156
-    @test Wflow.param(model, "vertical.snow")[41120] ≈ 4.21874475479126
-    @test Wflow.param(model, "vertical.tsoil")[41120] ≈ -1.9285825490951538
-    @test Wflow.param(model, "vertical.ustorelayerdepth")[1][1] ≈ 16.73013687133789
-    @test Wflow.param(model, "vertical.snowwater")[41120] ≈ 0.42188167572021484
-    @test Wflow.param(model, "vertical.canopystorage")[1] ≈ 0.0
-    @test Wflow.param(model, "vertical.zi")[1] ≈ 380.67793082060416
-    @test Wflow.param(model, "lateral.subsurface.ssf")[39308] ≈ 1.1614302208e11
-    @test Wflow.param(model, "lateral.river.q")[5501] ≈ 111.46229553222656
-    @test Wflow.param(model, "lateral.river.h")[5501] ≈ 8.555977821350098
-    @test Wflow.param(model, "lateral.river.volume")[5501] ≈ 249163.33754433296
-    @test Wflow.param(model, "lateral.land.q")[39626] ≈ 1.0575231313705444
-    @test Wflow.param(model, "lateral.land.h")[39626] ≈ 0.03456519544124603
-    @test Wflow.param(model, "lateral.land.volume")[39626] ≈ 19379.23274404319
+    @test Wflow.param(model, "lateral.river.reservoir.volume")[1] ≈ 2.7393418e7
+    @test Wflow.param(model, "vertical.satwaterdepth")[9115] ≈ 201.51429748535156
+    @test Wflow.param(model, "vertical.snow")[9115] ≈ 4.21874475479126
+    @test Wflow.param(model, "vertical.tsoil")[9115] ≈ -1.9285825490951538
+    @test Wflow.param(model, "vertical.ustorelayerdepth")[50069][1] ≈ 16.73013687133789
+    @test Wflow.param(model, "vertical.snowwater")[9115] ≈ 0.42188167572021484
+    @test Wflow.param(model, "vertical.canopystorage")[50069] ≈ 0.0
+    @test Wflow.param(model, "vertical.zi")[50069] ≈ 380.67793082060416
+    @test Wflow.param(model, "lateral.subsurface.ssf")[10606] ≈ 1.1614302208e11
+    @test Wflow.param(model, "lateral.river.q")[149] ≈ 111.46229553222656
+    @test Wflow.param(model, "lateral.river.h")[149] ≈ 8.555977821350098
+    @test Wflow.param(model, "lateral.river.volume")[149] ≈ 249163.33754433296
+    @test Wflow.param(model, "lateral.land.q")[10584] ≈ 1.0575231313705444
+    @test Wflow.param(model, "lateral.land.h")[10584] ≈ 0.03456519544124603
+    @test Wflow.param(model, "lateral.land.volume")[10584] ≈ 19379.23274404319
 end
 
 @testset "reducer" begin
@@ -210,7 +201,7 @@ end
     # test if the reverse index reverses the index
     linear_index = 100
     cartesian_index = indices[linear_index]
-    @test cartesian_index === CartesianIndex(115, 6)
+    @test cartesian_index === CartesianIndex(178, 7)
     @test reverse_indices[cartesian_index] === linear_index
 end
 
@@ -218,7 +209,7 @@ end
     @unpack vertical = model
     @test vertical.cfmax[1] ≈ 3.7565300464630127
     @test vertical.soilthickness[1] ≈ 2000.0
-    @test vertical.precipitation[100] ≈ 2.197660446166992
+    @test vertical.precipitation[49951] ≈ 2.197660446166992
 end
 
 config.input.vertical.cfmax = Dict("value" => 2.0)
@@ -236,7 +227,7 @@ model = Wflow.initialize_sbm_model(config)
     @unpack vertical = model
     @test vertical.cfmax[1] == 2.0
     @test vertical.soilthickness[1] ≈ 2000.0 * 3.0 + 100.0
-    @test vertical.precipitation[100] ≈ 1.5 * 2.197660446166992
+    @test vertical.precipitation[49951] ≈ 1.5 * 2.197660446166992
 end
 
 Wflow.close_files(model, delete_output = false)
@@ -248,3 +239,58 @@ Wflow.close_files(model, delete_output = false)
     ds = Wflow.create_tracked_netcdf(path)
     close(ds)  # path is removed on process exit
 end
+
+@testset "NetCDF read variants" begin
+    NCDataset(staticmaps_moselle_path) do ds
+
+        @test Wflow.is_increasing(ds[:x])
+        @test !Wflow.is_increasing(ds[:y])
+
+        @test Wflow.nc_dim_name(ds, :time) == :time
+        @test Wflow.nc_dim_name([:longitude], :x) == :longitude
+        @test Wflow.nc_dim_name([:lat], :y) == :lat
+        @test_throws ErrorException Wflow.nc_dim_name(ds, :not_present)
+
+        x = collect(Wflow.nc_dim(ds, :x))
+        @test length(x) == 291
+        @test x isa Vector{Union{Float64,Missing}}
+
+        @test Wflow.internal_dim_name(:lon) == :x
+        @test Wflow.internal_dim_name(:latitude) == :y
+        @test Wflow.internal_dim_name(:time) == :time
+
+        @test_throws ArgumentError Wflow.read_dims(ds["c"], (x = :, y = :))
+        @test_throws ArgumentError Wflow.read_dims(ds["LAI"], (x = :, y = :))
+        data, data_dim_order = Wflow.read_dims(ds["wflow_dem"], (x = :, y = :))
+        @test data isa Matrix{Union{Float32,Missing}}
+        @test data[end, end] === missing
+        @test data[125, 1] ≈ 643.547f0
+        @test data_dim_order == (:x, :y)
+
+        @test Wflow.dim_directions(ds, (:x, :y)) === (x = true, y = false)
+        @test Wflow.dim_directions(ds, (:y, :x, :layer)) === (y = false, x = true, layer = true)
+
+        data, dims = Wflow.permute_data(zeros(1, 2, 3), (:layer, :y, :x))
+        @test size(data) == (3, 2, 1)
+        @test dims == (:x, :y, :layer)
+        data, dims = Wflow.permute_data(zeros(1, 2), (:x, :y))
+        @test size(data) == (1, 2)
+        @test dims == (:x, :y)
+        @test_throws AssertionError size(Wflow.permute_data(zeros(1, 2, 3), (:x, :y)))
+
+        data = collect(reshape(1:6, (2, 3)))
+        # flip y, which is the second dimension
+        @test Wflow.reverse_data!(data, (y = false, x = true))[1, :] == [5, 3, 1]
+        # and mutate it back, the NamedTuple order should not matter
+        @test Wflow.reverse_data!(data, (x = true, y = false))[1, :] == [1, 3, 5]
+        # flip both dimensions at the same time
+        data = Wflow.reverse_data!(data, (x = false, y = false))
+        @test data[1, :] == [6, 4, 2]
+        @test data[:, 1] == [6, 5]
+
+        data = Wflow.read_standardized(ds, "wflow_dem", (x = :, y = :))
+        # since in this case only the second dimension needs reversing, we can easily do it manually
+        manual_fix = reverse(ds["wflow_dem"]; dims = 2)
+        @test all(data .=== manual_fix)
+    end
+end
diff --git a/test/run_hbv.jl b/test/run_hbv.jl
index e595302a5..52b259c55 100644
--- a/test/run_hbv.jl
+++ b/test/run_hbv.jl
@@ -26,12 +26,12 @@ end
 
 @testset "first timestep" begin
     hbv = model.vertical
-    @test hbv.tt[1] ≈ 0.0
+    @test hbv.tt[4377] ≈ 0.0
     @test model.clock.iteration == 2
-    @test hbv.soilmoisture[1] ≈ 134.35299682617188f0
-    @test hbv.runoff[1] ≈ 7.406898120121746f0
-    @test hbv.soilevap[1] == 0.0
-    @test hbv.snow[1] == 0.0
+    @test hbv.soilmoisture[4377] ≈ 134.35299682617188f0
+    @test hbv.runoff[4377] ≈ 7.406898120121746f0
+    @test hbv.soilevap[4377] == 0.0
+    @test hbv.snow[4377] == 0.0
 end
 
 # run the second timestep
@@ -39,28 +39,27 @@ model = Wflow.update(model)
 
 @testset "second timestep" begin
     hbv = model.vertical
-    @test hbv.soilmoisture[1] ≈ 134.35299682617188f0
-    @test hbv.runoff[1] ≈ 4.3533463f0
-    @test hbv.soilevap[1] == 0.0
-    @test hbv.snow[1] == 0.0
+    @test hbv.soilmoisture[4377] ≈ 134.35299682617188f0
+    @test hbv.runoff[4377] ≈ 4.3533463f0
+    @test hbv.soilevap[4377] == 0.0
+    @test hbv.snow[4377] == 0.0
 end
 
 @testset "overland domain" begin
     q = model.lateral.land.q_av
     land = model.lateral.land
     @test sum(q) ≈ 3278.9423039643552f0
-    @test q[500] ≈ 0.2348646912841589f0
-    @test land.volume[500] ≈ 2057.7314802432425f0
-    @test land.inwater[500] ≈ 0.027351853491789667f0
-    @test q[1566] ≈ 0.030463805623037462f0
+    @test q[10354] ≈ 0.2348646912841589f0
+    @test land.volume[10354] ≈ 2057.7314802432425f0
+    @test land.inwater[10354] ≈ 0.027351853491789667f0
     @test q[network.land.order[end]] ≈ 0.296794455575309f0
 end
 
 @testset "river flow" begin
     q = model.lateral.river.q_av
     @test sum(q) ≈ 51943.24072462076f0
-    @test q[500] ≈ 5.693405508134953f0
-    @test q[88] ≈ 8.839306661545436f0
+    @test q[651] ≈ 5.693405508134953f0
+    @test q[1056] ≈ 8.839306661545436f0
     @test q[network.river.order[end]] ≈ 363.5647479049738f0
 end
 
diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index 81e201f59..5636b035f 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -15,7 +15,7 @@ flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
 
     @test row.time == DateTime("2000-01-01T00:00:00")
     @test row.Q ≈ 7.815587720999557f0
-    @test row.volume ≈ 4.364251626782245f7
+    @test row.volume ≈ 2.7783643457819197f7
     @test row.temp_bycoord ≈ 2.3279826641082764f0
     @test row.temp_byindex ≈ 2.3279826641082764f0
     @test row.Q_6336050 ≈ 0.02388490515704142f0
@@ -44,25 +44,25 @@ end
     ds = model.writer.dataset_scalar
     @test ds["time"][1] == DateTime("2000-01-01T00:00:00")
     @test ds["Q"][:] ≈ [
-        0.023884907f0,
-        0.012411069f0,
-        0.004848609f0,
-        0.011474802f0,
-        0.0005262529f0,
-        0.013467698f0,
-        0.003408281f0,
-        0.09773275f0,
-        0.0021476103f0,
-        0.0026493936f0,
-        0.0008708129f0,
-        0.0007291489f0,
-        0.0021553952f0,
-        0.002229833f0,
-        0.003104503f0,
-        3.3423893f0,
-        1.35827f0,
+        1.6809974f0,
         5.9423304f0,
-        1.6809973f0,
+        1.35827f0,
+        3.3423896f0,
+        0.0031045028f0,
+        0.0022298328f0,
+        0.0021553952f0,
+        0.0007291489f0,
+        0.00087081286f0,
+        0.0026493934f0,
+        0.0021476103f0,
+        0.09773275f0,
+        0.0034082807f0,
+        0.013467698f0,
+        0.00052625296f0,
+        0.011474802f0,
+        0.004848609f0,
+        0.01241107f0,
+        0.023884906f0,
     ]
     @test ds["Q_gauges"].attrib["cf_role"] == "timeseries_id"
     @test ds["temp_index"][:] ≈ [2.3279827f0]
@@ -73,14 +73,14 @@ end
 @testset "first timestep" begin
     sbm = model.vertical
 
-    @test sbm.tt[1] ≈ 1.2999999523162842f0
+    @test sbm.tt[50069] ≈ 1.2999999523162842f0
 
     @test model.clock.iteration == 2
 
-    @test sbm.θₛ[1] ≈ 0.48343977332115173f0
-    @test sbm.runoff[1] == 0.0
-    @test sbm.soilevap[1] == 0.0
-    @test sbm.snow[1] ≈ 0.6029989752244306f0
+    @test sbm.θₛ[50069] ≈ 0.48343977332115173f0
+    @test sbm.runoff[50069] == 0.0
+    @test sbm.soilevap[50069] == 0.0
+    @test sbm.snow[50069] ≈ 0.6029989752244306f0
 end
 
 # run the second timestep
@@ -88,16 +88,16 @@ model = Wflow.update(model)
 
 @testset "second timestep" begin
     sbm = model.vertical
-    @test sbm.θₛ[1] ≈ 0.48343977332115173f0
-    @test sbm.runoff[1] == 0.0
-    @test sbm.soilevap[1] ≈ 0.005865651540305367f0
-    @test sbm.snow[1] ≈ 0.009696763863612956f0
+    @test sbm.θₛ[50069] ≈ 0.48343977332115173f0
+    @test sbm.runoff[50069] == 0.0
+    @test sbm.soilevap[50069] ≈ 0.005865651540305367f0
+    @test sbm.snow[50069] ≈ 0.009696763863612956f0
 end
 
 @testset "subsurface flow" begin
     ssf = model.lateral.subsurface.ssf
     @test sum(ssf) ≈ 6.368140761295825f16
-    @test ssf[network.land.order[1]] ≈ 3.0449782003445332f13
+    @test ssf[network.land.order[1]] ≈ 5.324295565183691f11
     @test ssf[network.land.order[end-100]] ≈ 7.855716879739626f11
     @test ssf[network.land.order[end]] ≈ 2.161246841709492f11
 end
@@ -113,18 +113,18 @@ end
 @testset "river flow" begin
     q = model.lateral.river.q_av
     @test sum(q) ≈ 2807.884971209105f0
-    @test q[4061] ≈ 0.0016288040314320486f0
-    @test q[5617] ≈ 7.338169165884175f0
+    @test q[1622] ≈ 0.0016288040314320486f0
+    @test q[43] ≈ 7.338169165884175f0
     @test q[network.river.order[end]] ≈ 0.00610520650626283f0
 end
 
 @testset "reservoir simple" begin
     res = model.lateral.river.reservoir
-    @test res.outflow[2] ≈ 0.2174998592483153f0
-    @test res.inflow[2] ≈ 50.170880189190626f0
-    @test res.volume[2] ≈ 2.776162917050312f7
-    @test res.precipitation[2] ≈ 0.1765228509902954f0
-    @test res.evaporation[2] ≈ 0.5372688174247742f0
+    @test res.outflow[1] ≈ 0.2174998592483153f0
+    @test res.inflow[1] ≈ 50.170880189190626f0
+    @test res.volume[1] ≈ 2.776162917050312f7
+    @test res.precipitation[1] ≈ 0.1765228509902954f0
+    @test res.evaporation[1] ≈ 0.5372688174247742f0
 end
 
 # set these variables for comparison in "changed dynamic parameters"
@@ -161,8 +161,8 @@ end
 # downstream from pit is at river index 1739, CartesianIndex(142, 85)
 @testset "river flow at and downstream of pit" begin
     q = model.lateral.river.q_av
-    @test q[1765] ≈ 8.133435566601927f0
-    @test q[1739] ≈ 0.008996045895155833f0
+    @test q[3908] ≈ 8.133435566601927f0
+    @test q[3920] ≈ 0.008996045895155833f0
 end
 
 # test changing forcing and cyclic LAI parameter
diff --git a/test/run_sediment.jl b/test/run_sediment.jl
index eae48b0fc..e93d6b60e 100644
--- a/test/run_sediment.jl
+++ b/test/run_sediment.jl
@@ -17,7 +17,7 @@ model = Wflow.update(model)
     @test eros.dmsand[1] == 200.0f0
     @test eros.dmlagg[1] == 500.0f0
     @test mean(eros.interception) ≈ 0.4767846753916875f0
-    @test mean(eros.soilloss) ≈ 0.00846987746234039f0
+    @test mean(eros.soilloss) ≈ 0.008596682196335555f0
 end
 
 # run the second timestep
@@ -26,10 +26,10 @@ model = Wflow.update(model)
 @testset "second timestep sediment model (vertical)" begin
     eros = model.vertical
 
-    @test mean(eros.soilloss) ≈ 0.07441259472436522f0
-    @test mean(eros.erosclay) ≈ 0.002198312957830507f0
-    @test mean(eros.erossand) ≈ 0.02515898751801652f0
-    @test mean(eros.eroslagg) ≈ 0.021613617766119256f0
+    @test mean(eros.soilloss) ≈ 0.07601393657280235f0
+    @test mean(eros.erosclay) ≈ 0.0022388720384961766f0
+    @test mean(eros.erossand) ≈ 0.02572046244407882f0
+    @test mean(eros.eroslagg) ≈ 0.022126541806118796f0
     @test mean(eros.TCsed) == 0.0
     @test mean(eros.TCsilt) ≈ 1.0988158364353527f6
     @test mean(eros.TCsand) ≈ 1.0987090622888755f6
@@ -40,15 +40,15 @@ end
 @testset "second timestep sediment model (lateral)" begin
     lat = model.lateral
 
-    @test mean(lat.land.inlandsed) ≈ 0.0727443455767198f0
-    @test mean(lat.land.inlandclay) ≈ 0.002190037611542877f0
-    @test mean(lat.land.inlandsand) ≈ 0.024512058443074893f0
-    @test mean(lat.land.olclay) ≈ 0.006212441628146587f0
+    @test mean(lat.land.inlandsed) ≈ 0.07463801685030906f0
+    @test mean(lat.land.inlandclay) ≈ 0.0022367786781657497f0
+    @test mean(lat.land.inlandsand) ≈ 0.02519222037812127f0
+    @test mean(lat.land.olclay) ≈ 0.006443036462118322f0
 
-    @test mean(lat.river.SSconc) ≈ 0.8454728220852514f0
-    @test mean(lat.river.inlandclay) ≈ 0.019393435838709512f0
+    @test mean(lat.river.SSconc) ≈ 0.8259993252994058f0
+    @test mean(lat.river.inlandclay) ≈ 0.01980468760667709f0
     @test lat.river.h_riv[network.river.order[end]] ≈ 0.006103649735450745f0
-    @test lat.river.outclay[1] ≈ 3.3447705344961993f-6
+    @test lat.river.outclay[5649] ≈ 2.359031898208781f-9
 end
 
 Wflow.close_files(model)
diff --git a/test/runtests.jl b/test/runtests.jl
index 046b9b927..c32e1b579 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,6 @@
 ## load test dependencies and set paths to testing data
 using Dates
+using Downloads
 using LightGraphs
 using NCDatasets
 using StaticArrays
@@ -27,7 +28,7 @@ function testdata(version, source_filename, target_filename)
     target_path = joinpath(datadir, target_filename)
     base_url = "https://github.com/visr/wflow-artifacts/releases/download"
     url = string(base_url, '/', string('v', version), '/', source_filename)
-    isfile(target_path) || download(url, target_path)
+    isfile(target_path) || Downloads.download(url, target_path)
     return target_path
 end
 
@@ -37,16 +38,16 @@ staticmaps_lahn_path = testdata(v"0.2.1", "staticmaps-lahn.nc", "staticmaps-lahn
 forcing_moselle_path = testdata(v"0.2", "forcing-2000.nc", "forcing-moselle.nc")
 forcing_lahn_path = testdata(v"0.2", "forcing-lahn.nc", "forcing-lahn.nc")
 forcing_moselle_sed_path =
-    testdata(v"0.2", "forcing-moselle-sed.nc", "forcing-moselle-sed.nc")
+    testdata(v"0.2.3", "forcing-moselle-sed.nc", "forcing-moselle-sed.nc")
 staticmaps_moselle_sed_path =
-    testdata(v"0.2.2", "staticmaps-moselle-sed.nc", "staticmaps-moselle-sed.nc")
+    testdata(v"0.2.3", "staticmaps-moselle-sed.nc", "staticmaps-moselle-sed.nc")
 instates_moselle_sed_path =
     testdata(v"0.2", "instates-moselle-sed.nc", "instates-moselle-sed.nc")
 instates_moselle_path = testdata(v"0.2.1", "instates-moselle.nc", "instates-moselle.nc")
 forcing_sbm_gw_path =
     testdata(v"0.2.1", "forcing-sbm-groundwater.nc", "forcing-sbm-groundwater.nc")
 staticmaps_sbm_gw_path =
-    testdata(v"0.2.1", "staticmaps-sbm-groundwater.nc", "staticmaps-sbm-groundwater.nc")
+    testdata(v"0.2.2", "staticmaps-sbm-groundwater.nc", "staticmaps-sbm-groundwater.nc")
 lake_sh_1_path = testdata(v"0.2.1", "lake_sh_1.csv", "lake_sh_1.csv")
 lake_sh_2_path = testdata(v"0.2.1", "lake_sh_2.csv", "lake_sh_2.csv")
 lake_hq_2_path = testdata(v"0.2.1", "lake_hq_2.csv", "lake_hq_2.csv")
diff --git a/test/sbm_config.toml b/test/sbm_config.toml
index 76bd598c2..c043faeaa 100644
--- a/test/sbm_config.toml
+++ b/test/sbm_config.toml
@@ -165,7 +165,7 @@ parameter = "vertical.temperature"
 location = "temp_byindex"
 name = "temp_index"
 index.x = 100
-index.y = 50
+index.y = 264
 parameter = "vertical.temperature"
 
 [csv]
@@ -190,7 +190,7 @@ parameter = "vertical.temperature"
 [[csv.column]]
 header = "temp_byindex"
 index.x = 100
-index.y = 50
+index.y = 264
 parameter = "vertical.temperature"
 
 [[csv.column]]
diff --git a/test/sbm_gw.toml b/test/sbm_gw.toml
index ad9d10d84..00d9dc43a 100644
--- a/test/sbm_gw.toml
+++ b/test/sbm_gw.toml
@@ -156,7 +156,7 @@ parameter = "vertical.temperature"
 location = "temp_byindex"
 name = "temp_index"
 index.x = 100
-index.y = 50
+index.y = 264
 parameter = "vertical.temperature"
 
 [csv]
@@ -181,7 +181,7 @@ parameter = "vertical.temperature"
 [[csv.column]]
 header = "temp_byindex"
 index.x = 100
-index.y = 50
+index.y = 264
 parameter = "vertical.temperature"
 
 [[csv.column]]
diff --git a/test/sediment_config.toml b/test/sediment_config.toml
index 38835b4ea..8c2f10e3c 100644
--- a/test/sediment_config.toml
+++ b/test/sediment_config.toml
@@ -53,8 +53,8 @@ forcing = [
   "vertical.interception",
   "vertical.precipitation",
   "vertical.q_land",
-  "lateral.land.h_riv",
-  "lateral.land.q_riv",
+  "lateral.river.h_riv",
+  "lateral.river.q_riv",
 ]
 
 cyclic = ["vertical.leaf_area_index"]
@@ -86,11 +86,11 @@ resareas = "wflow_reservoirareas"
 lakeareas = "wflow_lakeareas"
 
 [input.lateral.land]
-h_riv = "h"
-q_riv = "q"
 slope = "Slope"
 
 [input.lateral.river]
+h_riv = "h"
+q_riv = "q"
 cbagnold = "c_Bagnold"
 d50 = "D50_River"
 d50engelund = "D50_River"