API Reference

waterflows(dem, cellarea=fill!(similar(dem),1);
            drain_pits=true,
            bnd_as_sink=true,
            nan_as_sink=true,
            extra_sinks=CartesianIndex{2}[],
            extra_barriers=CartesianIndex{2}[],
            flowdir_fn=d8dir_feature,
            feedback_fn=nothing)

Water flow routing according to the D8 algorithm. Local minima are filled, by default, using a breach-type algorithm, this means that the input DEM does not need to be pre-filled.

args:

• dem – the DEM (or hydro-potential); array
• cellarea=fill!(similar(dem),1) – the source per cell, defaults to 1.
  • if cellarea is negative in places, flux may go to zero but not below.
  • in areas where no routing takes place, typically NaNs in the dem, cellarea is ignored. This may affect mass-conservation.
  • if using physical units then use a volumetric flux per cell, e.g. m3/s.
  • Alternatively, cellarea can be a tuple of arrays. Then they are treated/routed separately, for instance (water, tracer). All quantities need to be extensive (i.e. additive, e.g. use internal energy and not temperature)

kwargs:

• drain_pits – whether to route through pits (true)
• bnd_as_sink (true) – whether the domain boundary should be sinks, i.e. adjacent cells can drain into them, or whether to ignore them.
• nan_as_sink (true) – whether NaN cells in the DEM should make adjacent cells a sink. Note that on the NaN-cell itself no routing occurs (i.e. a BARRIER-cell).
• extra_sinks=CartesianIndex{2}[] – additional cells that act as sinks, i.e. outlets where flow leaves the active domain.
• extra_barriers=CartesianIndex{2}[] – additional cells that act as barriers, i.e. cells that are excluded from routing and do not conduct flow.
• flowdir_fn=d8dir_feature – the routing function. Defaults to the built-in d8dir_feature function but could be customized
• feedback_fn – function which is applied to area-value(s) at each cell once all water of the cell has been accumulated but before the water is routed further downstream. Signature (uparea, ij, dir) ->new_uparea--> for example,(uparea, ij, dir) -> max(uparea, 0)` would ensure that all upareas are non-negative.

Returns a NamedTuple with fields:

• area – upslope area (or a tuple of upslope areas if cellarea is a tuple too)
• slen – length of stream to the farthest source (number of cells traversed)
• dir – flow direction at each location
• nout – whether the point has outlflow. I.e. nout[I]==0 –> I is a pit
• nin – number of inflow cells
• sinks – location of sinks as Vector{CartesianIndex{2}}
• pits – location of pits as Vector{CartesianIndex{2}}
• c – catchment map (color numbers ∈ 1:length(sinks) are for sinks, others for pits)
• bnds – boundaries between catchments. The boundary to the exterior/NaNs is not in here.
• flowdir_extra_output – extra output of the flowdir_fn, which is nothing for the default

fill_dem(dem, sinks, dir; small=0)

Fill the pits of a DEM (apply this after applying "drainpits", which is done by default in waterflows). Returns the filled DEM.

Notes:

• this is not needed as pre-processing step to use the flow-routing function waterflows.
• routing on a filled DEM will not produce exactly the same flow pattern: on the shores of lakes streams which entered the lake can now go the other way. I suspect on most DEMs the differences will be very minimal.
• This uses a tree traversal to fill the DEM. It does it depth-first (as it is easier) which may lead to a stack overflow on a large DEM.

Direction number/constant: indicating no flow, this is a "pit", i.e. a local minimum or a cell in a completely flat area. Note: PITs are sinks if drain_pits==false.

Direction number/constant: indicating a cell where water disappears, typically located at the domain boundary (when setting bndassink=true) adjacent to NaN-cells of the DEM.

Direction number/constant: indicating no flow into or out of this cell. All DEM cells which are NaNs map to this. If NaNs are not acting as sinks, then water routes around these cells.

catchment(dir, ij)

Calculates the catchment of one or several grid-point(s) ij. If desired, the catchment boundary can be calculated with make_boundaries([c], [1]).

Input

• dir – direction field
• ij – index of the point (2-Tuple, or CartesianIndex) or a Vector{CartesianIndex} for several

Returns

• catchment – BitArray

Tip: only being off by one grid-point can make the difference between a tiny and a huge catchment!

See also: catchments

catchments(dir, sinks::Union{Vector{Vector{CartesianIndex{2}}}, Vector{<:CartesianIndices{2}}}, dem=nothing;
                check_sinks_overlap=true)

Make a map of catchments from different (non-overlapping) sinks.

See also: catchment

catchment_flux(cellarea, c, color) = sum(cellarea[c.==color])
catchment_flux(cellarea, c::Union{BitArray, Matrix{Bool}})

The total flux, i.e. input, in one catchment.

prune_catchments(catchments, minsize; val=0)

Sets all catchments which are smaller than minsize to val (=0). Catchments with number <1 are ignored.

Note that, as the catchments are re-numbered, the number will not correspond to the sinks and pits anymore.

d8dir_feature(dem, bnd_as_sink, nan_as_sink, extra_sinks=CartesianIndex{2}[], extra_barriers=CartesianIndex{2}[])

D8 directions of a DEM and drainage features (nin & nout).

Elevations with NaN map to dir==BARRIER, cells around them will be set to SINK if nan_as_sink==true or receive no special treatment otherwise.

The argument bnd_as_sink determines whether cells at the domain boundary act as sinks.

The extra_sinks and extra_barriers arguments can be used to manually mark additional cells as sinks or barriers, respectively. Cells in extra_sinks act as outlets where flow leaves the domain. Cells in extra_barriers are excluded from routing and do not conduct flow.

Return

• dir - direction, encoded as dirnums
• nout - number of outflow cells of a cell (0 or 1)
• nin - number of inflow cells of a cell (0-8)
• sinks - location of sinks as a Vector{CartesianIndex{2}} (sorted) (here dir==SINK)
• pits - location of pits as a Vector{CartesianIndex{2}} (sorted) (here dir==PIT)
• dem - DEM, unchanged
• flowdir_extra_output – nothing (not used by this function, but could be by custom ones)

flowrouting_catchments(dir, pits, cellarea, feedback_fn)

Recursively calculate flow-routing and catchments from

• dir - direction field
• pits - pit coordinates
• cellarea - water input

Returns:

• area – upslope contributing area
• slen – stream length, i.e. the number of cells traversed along the longest upstream flow path
• c – catchment map (Matrix{Int}); c==0 corresponds to NaN/BARRIER regions where no water flows into.

Note: this function may cause a stackoverflow on very big catchments.

drainpits!(dir, nin, nout, sinks, pits, ctch, bnds, dem)

Update in place the direction field such that it drains pits. This is done by reversing the flow connecting the lowest point (which can drain) on the catchment boundary to the pit for each such catchment.

Update in place dir, nin, nout, pits (sorted), c; returns an empty bnds

make_boundaries(catchments, pit_colors, bnds=Origin(firstindex(pit_colors))([CartesianIndex{2}[] for i in 1:length(pit_colors)]) )

Make vectors of boundary cells for catchments of pit_colors (as the name suggests, typically just the pit-catchment colors.) Assumes that pit_colors is sorted.

Note that cells along the edge of the domain[1] as well as cells only bordering BARRIER-cells are not included (because the algorithms do not need to traverse them).

Return:

• bnds – Vector of Vector{CartesianIndex{2}} containing the cells which are on the boundary of said catchment.

[1] Note that when bnd_as_sink==true then no cells along the boundary will belong to a pit-catchment.

make_flowfeatures(dir)

Makes nout, nin, sinks, pits by iterating once over the dir array

dir2ind(dir, map_special_to_PIT=false)

Translate a D8 direction number into a CartesianIndex (i.e. a flow vector), also maps SINK to CartesianIndex(0,0).

If map_special_to_PIT==true, then any dir>9 is mapped to CartesianIndex(0,0).

Translate a D8 direction number into a 2D vector.

Translate a CartesianIndex, giving the offset, into a direction number

Adjustment factor for gradient in diagonal vs non-diagonal

Tests whether a cell J with flowdir dirJ flows into cell I.

on_outer_boundary(ar, I::CartesianIndex)

Check whether on outer boundary of ar

Return CartesianIndices corresponding to the 8 neighbors and the point itself

Show an array with its axes oriented such that they correspond to x and y direction.

Direction numbers. E.g. dirnums[1,1] will return the number corresponding to the direction top-left.

Note, I use the conversion that the x-axis corresponds to the row of the matrix and the y-axis the columns. To print them in this "normal" coordinate system use showme

Translation from dirnums to CartesianIndex

Update dir, nin, and nout such that flow at P1 is now from P1 to P2.

It can potentially modify dir, nin, and nout at three locations:

• P1: dir, nout, nin
• P2: nin
  • if allow_P2_pit==true, then P2's dir and nout can also be modified.
• P3 (previous receiver cell of P1): nin

Note that:

• if P1 and P2 lie in the same catchment, then P3 is also in that catchment.
• if flow was from P2 to P1, then P2 has to become a pit (PIT) to keep dir consistent.

boxcar(A::AbstractArray, window::AbstractArray, weights)

Moving average filter with spatially dependent filter window and weighting of cells. Filters over ±window.

catchment_sinks(ctch_polygons::Vector{Vector{T}} where T<:Point2, routing_mask, x, y, sinks=get_boundary_cells(routing_mask, true);
                check=true, threshold_dist=Inf,
                D9_iters=4, D25_iter=2)
catchment_sinks(ctch_polygons::AbstractVector, routing_mask, x, y, sinks=get_boundary_cells(routing_mask, true);
                kws...)

Calculate the catchment sink cells for basins given by polygons. The output is intended to be fed to waterflowssubglacial via its `ctchsinks` argument. This is intended to mark the cells of grounding line of a catchment as sinks.

Correctness check can be done with the check_catchment_sinks function, which is done by default.

Input:

• ctch_polygons::Vector{AbstractVector{<:Point2}} – a Vector of catchment polygons, presumably sourced from the internet. Need to be non-overlapping
• routing_mask - cells where water is routed set to true; assumed that boundary are all sinks
• x,y – coordinate vectors
• sinks – sink cells, by default calculated as all cells on the border of the routing_mask

KW-args

• check=true – do some checks (takes a bit of time)
• threshold_dist=Inf – probably set to something lower

Output

• ctch_sinks – a vector

Algo:

• get all boundary cells (make three vectors: cartesian-index, closest catchment, distance)
• calc signed_distance for all pixels for all catchments
• assign groundingline cells to catchment to which it is closest
• do some cleanup...

check_catchment_sinks(ctch_sinks_map, ctch_polygons, sinks)

d8dir_pressmelt((phi, phim), bnd_as_sink, nan_as_sink, gamma, avoid_sc)

D8 directions of a DEM and drainage features, taking pressure-melting point effects into account (i.e. flow deflection and (almost) no flow when the supercooling threshold is met).

Elevations with NaN map to dir==BARRIER, cells around them will be set to SINK if nanassink==true or receive no special treatment otherwise.

The argument bnd_as_sink determines whether cells at the domain boundary act as sinks.

Places where supercooling occurs can be treated in two ways: water routes around them (avoidsc==true, this sets supercooled cells to BARRIER) or water still flows through then (avoidsc==false).

Args

• gamma – the Röthlisberger constant. Set to 0 to get no supercooling, typical value is GAMMA.
• avoid_sc – if set, then they are set as barrier cells

Return

• dir - direction, encoded as dirnums
• nout - number of outflow cells of a cell (0 or 1)
• nin - number of inflow cells of a cell (0-8)
• sinks - location of sinks as a Vector{CartesianIndex{2}} (sorted) (here dir==SINK)
• pits - location of pits as a Vector{CartesianIndex{2}} (sorted) (here dir==PIT)
• dem - DEM, unchanged
• flowdirextraoutput – location of the supercooled cells (more precisely: the outflow edge of such a cell is supercooled)

get_boundary_cells(mask, value)

Get all the cells of mask of value value which border on cells of different value.

getkappa_D8(R, lambda::Int, gamma)

Return kappa (angle between -∇ϕ and Q) which maximizes melt for a given ratio R = ϕₘ/ϕ. A crude approximation which is exact for D8.

Notes:

• kappa ∈ {-π/2, -π/4, 0, π/4, π/2}
• Only implemented for gamma ∈ {0, -0.31}.
• Takes lambda, the angle between surface and bed slope, in units of π/4, and returns kappa in those units too.
• It important is to check that water is not flowing phi-uphill after the deflection. If it is the case the deflection needs to be reduced/removed in the calling function.

Return CartesianIndices corresponding to the 24 neighbors and the point itself.

lake_depth(phi_filled, phi_orig; fixed_surface=true, rhow=RHOW, rhoi=RHOI)

The lake depth can be calculated assuming

1. adjusting the ice surface (fixed_surface=false)
2. adjusting the ice thickness only but keeping the ice surface fixed (fixed_surface=true)

The formulas are, respectively:

(phifilled - phioriginal)

(phifilled - phioriginal) * (rhow / (rhow - rhoi) )

make_feedback_fn(phi, phim, gamma, dx)

Make the feedback_fn as used by WWF with the appropriate closures and taking care of https://github.com/JuliaLang/julia/issues/15276 (if needed)

It returns:

• total discharge
• extra discharge due to dissipation and pressure melting
• dissipation melt rate (m/s)
• pressure melt rate (m/s)

mask_contiguous(mask, IJ, out=similar(mask, Bool))

Create a new mask (or update the optional 3rd argument) which

• marks all points connected to IJ and
• where each masked point has the same value as mask[IJ]

melt_rates(Q::Number, phi::AbstractMatrix, phim::AbstractMatrix, gamma, ij::CartesianIndex, dir::AbstractMatrix, dx::Number, rhow::Number)

Calculate melt rates due to dissipation of potential energy and due to pressure melting point effects given by

(-Q ∇ϕ - γ Q ∇pw)/(GRAV * L),

for ϕ and pw given in [m H2O].

This is calculated for the water flow from cell ij to its downstream cell.

This is used as part of the feedback_fn in WWF.waterflows.

Returns

• dissipation melt rate (first term) [m/s]
• pressure melt rate (first term) [m/s]

Note: the dissipation melt rate is negative in places where water is routed out of depressions as the water flows up the hydraulic potential (this is due to the breach type algorithm in drainpits!). However, freezing will balance the extra melt occurring when the water descends into the filled local minimum (pit).

phi_fn(surface, thick, floatfrac=1; rhow=RHOW, rhoi=RHOI)

Shreve hydraulic potential [m H2O].

signed_distance(p::Point2{T}, poly::AbstractVector{Point2{T}}) signed_distance(p::Point2, poly::AbstractVector{<:Point2})where {T}

Returns the distance of a point p to the polygon poly.

smooth_surface(x, y, surfdem, beddem, icethicknesses, mask= surfdem.>=beddem; minwindow=0)

Smooth the surface DEM of ±icethicknesses. Only smoothes where there is ice and only using ice-covered cells.

waterflows_subglacial(surfdem, beddem, dx, 
                           floatfrac=1,
                           source=ones(size(surfdem)),
                           mask=trues(size(surfdem));
                           gamma=-0.31,
                           ctch_sinks=[],
                           rhow = 1000.0,
                           rhoi = 910.0,
                           drain_pits=true,
                           bnd_as_sink=true,
                           nan_as_sink=true)

Does the water flow routing according the D8 algorithm for a subglacial setting using the Shreve-potential for routing. Utilizes WhereTheWaterFlows.waterflows.

args:

• surfdem, beddem – the surface and bed DEM
• floatfrac – flotation fraction
• dx – grid size (must be equal in x and y-direction)
• source=ones(size(dem)) – the source per cell, defaults to 1. If using physical units then use a source in volume per unit area, e.g. m/s
• mask – routing mask: where to do routing

kwargs:

• drain_pits – whether to route through pits (true)
• bndassink (true) – whether the domain boundary be sinks, i.e. i.e. adjacent cells can drain into them or whether to ignore them.
• nanassink (true) – whether NaNs are sinks
• gamma – Röthlisberger constant (-0.31, note it's negative)
• avoid_sc (false) – If set to true, then supercooled cells become BARRIERs. Note that this then also means that source in those cells will just vanish, i.e. mass is lost. Probably the default ==false is more physical as waterflow should continue just not in R-channel. Thus recommended to set to false.
• rhow, rhoi – mean density of water and glacier ice (1000, 910)
• ctch_sinks – List of sink-areas for which catchments will be calculated. Each sink can be given as a vector of CartesianIndex or as CartesianIndices. Note, sink-areas and catchments can be overlapping, if desired.

Returns one nested NamedTuple with keys:

• routing: fields area, slen, dir, nout, nin, sinks, pits, c, bnds, phi
  • routing.area: named tuple with total, extra, dissipation_melt_rate, pressure_melt_rate
• pressmelt: fields sc_locs, kappas, dir_og
• lakes: fields depth_fixed_surface, depth_free_surface
• sink_catchments: fields masks, fluxes
  • sink_catchments.fluxes: named tuple with total, dissipation, pressmelt

Type to hold the uncertainty information for a field.

• absuc=0 – absolute uncertainty as std (array or scalar)
• reluc=0 – relative (to local field value) uncertainty as std (array or scalar)
• correlation_length=1 – ditto (scalar)
• abs_bounds=(-Inf,Inf) – values of the GRF are constrained to lie within these bounds (simple crop, after being scaled with absuc and reluc)

depad(out::Array, nx::Int, ny::Int, padding::Tuple) -> Array

Removes padding from an array.

distance_squared(nx, ny) -> Array

Computes the square of the Euclidean distance from the center in a 2D grid.

eq18(nx::Int, ny::Int, len::Number, kernel::Array) -> Array

Draws a sample of the Gaussian random field using the precomputed kernel.

exponential_kernel(nx::Int, ny::Int, len::Number) -> Array

Generates an exponential kernel based on the input dimensions and length scale.

find_next_2357(n::Integer) -> Integer

Finds the next integer following n that has only 2, 3, 5, and 7 as its prime factors. This is useful for optimization in FFT operations, which are more efficient on sizes factorable into these primes.

gaussian_kernel(nx::Int, ny::Int, len::Number) -> Array

Generates a Gaussian kernel based on the input dimensions and correlation length.

make_field_realization(field, sampler, uncert::Uncertainty)

Make a realization of a field.

make_fns_subaerial(dx,
                    dem, dem_uc,
                    source, source_uc,
                    ctch_sinks;
                    drain_pits=true,
                    bnd_as_sink=true,
                    nan_as_sink=true)

Build model, sample, and reduce! functions for stochastic subaerial routing.

Args

• dx – grid spacing
• dem, source – baseline elevation and source fields
• dem_uc, source_uc – corresponding Uncertainty objects
• ctch_sinks – sink groups used for catchment masks/flux aggregation

Kwargs

• drain_pits, bnd_as_sink, nan_as_sink – forwarded to WhereTheWaterFlows.waterflows

Return

• model(dem, source) -> returns (; dem, dx, source), waterflows(...)
• sample() -> returns one sampled (dem_, source_)
• reduce! -> three-method reduction function (reduce!(), reduce!(aggr, out), reduce!(aggr)).

make_fns_subglacial(dx,
                      surfdem, surfdem_uc,
                      beddem, beddem_uc,
                      floatfrac, floatfrac_uc,
                      source, source_uc,
                      ctch_sinks;
                      mask=mask::AbstractMatrix=fill!(similar(surfdem, Bool), true),
                      gamma=0.0,
                      min_lake_depth=10.0, # default min lake depth under which value are not aggregated
                      rhow=RHOW, rhoi=RHOI)

Build model, sample, and reduce! functions for stochastic subglacial routing which are then used in map_mc.

Args

• dx – grid spacing
• surfdem, beddem, floatfrac, source – baseline fields used for sampling and routing
• floatfrac must be an array with the same size as surfdem. A scalar is not accepted here; use fill(value, size(surfdem)...) for a spatially uniform flotation field.
• surfdem_uc, beddem_uc, floatfrac_uc, source_uc – corresponding Uncertainty objects
• ctch_sinks – sink groups used for catchment masks/flux aggregation

Kwargs

• mask – routing mask passed to Subglacially.waterflows_subglacial
• gamma – pressure-melting deflection parameter; defaults to 0.0 when omitted. Otherwise set to WWF.GAMMA for a standard value.
• min_lake_depth – threshold used for lake-occurrence masks and lake-volume vectors
• rhow, rhoi – water and ice density overrides

Return:

• model(surf, bed, floatfrac, source) -> runs one model realisation returning (;surf, bed, dx, floatfrac, source), waterflows_subglacial(...)
• sample() -> surf, bed, float, source
• reduce! -> three-method reduction function (reduce!(), reduce!(aggr, out), reduce!(aggr)).

make_grf_sampler(nx::Int, ny::Int, kernel_fn::Function, len::Number) -> Function

Creates a sampler function for generating Gaussian random fields (GRF) using pre-specified parameters.

KW-arguments:

• fftw_plan: if not nothing then FFTW is planned for improved performance. Use :estimate or :measure (recommended) to use a plan.

make_kernel(kernel_fn::Function, nx::Int, ny::Int, len::Number) -> Array

Precomputes and returns the Fourier transform of a kernel function over specified dimensions and correlation length.

map_mc(model, sample, reduce!, n)

A general Monte Carlo map function

Args

• model – the model function, runs with output = model(sample()...)
• sample – input = sample() returns a new sample to use as input to model
• reduce! – aggregates the model output with reduce!(aggr, output). Also initializes the aggregate storage with aggr = reduce!() and finalizes it with reduce!(aggr).
• n – number of samples to take (must be <= 2048; Float16 aggregation for catchment frequencies loses precision for larger counts)

Thread safety: only reduce! is allowed to not be thread-safe.

pad(nx, ny, len) -> Tuple

Calculates the padding needed on both sides to facilitate efficient and non-periodic FFT processing.

Arguments

• nx, ny::Int`: Original x and y dimension.
• len::Number: correlation length (in number of grid cells)

Returns

• Tuple: A tuple containing the new x dimension, new y dimension, and the padding to be applied on each side to reach these new dimensions.

plt_area(x, y, area; prefn=log10, sinks=[], threshold=Inf, colorbar=true,
         colorbar_label="log₁₀(upstream_area)", colorbar_kwargs=(;))

Plot uparea, or another variable.

Kwargs:

• pre-proc with prefun, typically and by default this is log10
• if sinks (or pits) are passed, plot as points
• threshold the area: do not plot pixels with area below threshold
• add a colorbar with colorbar=true
• set colorbar label with colorbar_label
• pass additional kwargs to Colorbar via colorbar_kwargs

plt_bnds(x, y, bnds)

Plot boundary points.

plt_catchments(x, y, c; minsize=0)

Plot catchments. Catchments below minsize size are not plotted. With the default minsize=0 all catchments are plotted.

Note, minsize>0 can be quite slow to compute.

plt_dir( x, y, dir)

Plot dir as flow field

plt_it(x, y, waterflows_output, dem)

Plot DEM, uparea, flow-dir

plt_sinks(x, y, sink_pits)

Plot sinks or another vector of CartesianIndices