Funmixer - Unmixing nested observed concentrations in river networks for source regions

This repository implements an efficient solution to the unmixing of nested concentrations in a (river) network using convex optimisation. The method is described in full in our article, Using Convex Optimization to Efficiently Apportion Tracer and Pollutant Sources From Point Concentration Observations, published in Water Resources Research (DOI: 10.1029/2023WR036159). A copy of the paper is included in this repository at docs/Barnes_Lipp_2024.pdf, and is the reference for the mathematics, assumptions and implementation details of the method.

Data input assumptions

The algorithm requires:

1. A GDAL readable raster of D8 flow directions. We use the ESRI/Arc D8 convention of representing directions with increasing powers of 2 (i.e., 1, 2, 4, 8 etc.) with sink pixels indicated by 0. We assume that every cell in the domain eventually flows into a sink node within the domain (or is itself a sink node). This assumption requires that every boundary pixel is set to be a sink.

2. A .csv file which contains the names, locations and geochemical observations at the sample sites. Sample names (e.g., 'SampleA') are expected in column 1 and the x and y-coordinates of the sample sites in columns 2 and 3. The x and y-coordinates of the sample sites need to be in the same reference system as the D8 raster (for the supplied example data this is EPSG:27700). It is assumed that the sample sites have already been manually aligned onto the drainage network. Subsequent columns contain the name of a given tracer (e.g., Mg) and their concentrations (arbitrary units).

Tracer values may also be reported as below the detection limit by writing them as a string of the form <X, where X is the detection limit (e.g. <0.5). Different elements, and different samples of the same element, may carry different detection limits. Such tables are read with get_bdl_element_obs and unmixed with BDLSampleNetworkUnmixer; see Below detection limit data below. Blank cells remain missing data, which is not the same thing as a censored value.

funmixer does include some basic data preprocessing functions that can be used to align the sample sites to the drainage network and fix the boundary conditions of the D8 raster. An example of use is given in the examples/ directory. Example, valid, datasets are contained in data/d8.tif and data/sample_data.csv.

The example dataset

The example data supplied in data/ are a real geochemical survey of north-east Scotland (the Cairngorms and surrounding catchments: the Spey, Dee, Don, Deveron and Tay). Both the D8 raster (data/d8.tif, 50 m resolution) and the sample site coordinates (data/sample_data.csv) are georeferenced to the OSGB36 / British National Grid (EPSG:27700) projection, so the x_coordinate and y_coordinate columns are eastings and northings in metres. The files carry their coordinate reference system, so they can be loaded directly into GIS software such as QGIS and will overlay correctly. data/noisy_sample_data.csv is the same table with the sample sites deliberately displaced from the drainage network, and is used by examples/preprocessing.py to demonstrate snap_to_drainage. data/BDL_sample_data.csv is the same table again with the lowest decile of each element replaced by censored <X entries, and is regenerated by python data/make_bdl_sample_data.py.

Some common data input problems can be solved by:

• Checking that there is not trailing white-space at the end of the sample site data table.
• Ensuring that the D8 flow-direction raster is in the same reference system as the sample site coordinates.
• Ensuring that the D8 raster is surrounded by a boundary of sink pixels.

Below detection limit (BDL) data

Geochemical surveys routinely report a concentration only as <X: the concentration is somewhere at or below the detection limit X. Passing X in as though it were a measurement asserts more than the datum says, and biases recovered sources upwards.

BDLSampleNetworkUnmixer behaves exactly like SampleNetworkUnmixer except that censored observations use a one-sided misfit:

observation	misfit
detected, c_obs	max(c_pred / c_obs, c_obs / c_pred)
censored, <X	max(1, c_pred / X)

The censored form is flat for any prediction at or below the detection limit, encoding the idea that any value below the limit is equally plausible, and grows in the usual relative way above it. When the geometric mean used internally to normalise the problem is computed, censored values enter at half their detection limit, following standard (if imperfect) practice for censored data. The misfit itself always compares against the full limit.

Regularization is strongly recommended for censored data. The one-sided misfit is still convex, so the problem remains a well-formed convex program and the solver will usually report optimal. What it loses is a unique minimiser: below a detection limit the objective is flat, so many source concentrations fit equally well and the solver may return an arbitrary point from that flat region. In the pathological limit where every observation in a network is censored, the data constrain the sources only from above. Regularization pulls each source towards the geometric mean of the observations and restores a well-posed answer, so constructing a BDLSampleNetworkUnmixer with use_regularization=False raises a warning.

import pandas as pd
import funmixer

sample_network, labels = funmixer.get_sample_graph(
    flowdirs_filename="data/d8.tif", sample_data_filename="data/BDL_sample_data.csv"
)
# Read as strings so that censored entries such as "<3821.8" survive
obs_data = pd.read_csv("data/BDL_sample_data.csv", dtype=str, keep_default_na=False)
element_data = funmixer.get_bdl_element_obs("Mg", obs_data)

problem = funmixer.BDLSampleNetworkUnmixer(
    sample_network=sample_network, observation_data=element_data
)
solution = problem.solve(element_data, regularization_strength=1e-3)

Each observation is a BDLObservation(value, is_bdl), where value holds the concentration or, for a censored entry, the detection limit. Which sites are censored is fixed when the problem is built, so observation values can still be changed between solves without rebuilding, but the censoring pattern cannot.

Installation

Funmixer has some python package dependencies so you will need to install these before you can run Funmixer. You can do this by either creating a conda environment or by creating an apptainer container

Making the Funmixer conda environment

The following assumes a UNIX operating using the conda package manager. conda is preferred as it (in general) allows for easier installation of the gdal dependency than, for instance, pip. Whilst this package was developed on a UNIX system, the following commands (or similar) should be possible on a Windows OS using an Anaconda prompt.

First, clone the repository into a local directory:

git clone https://github.com/AlexLipp/funmixer/ [LOCAL_DIRECTORY]

A conda environment file (requirements.yaml) is provided containing the python dependencies. A conda environment entitled funmixer can be generated from it by running:

conda env create -f requirements.yaml

The environment can then be activated using

conda activate funmixer

Next, install the funmixer python package using:

pip install -e .

This command installs the funmixer python package that can then be imported as normal (e.g., import funmixer).

Problem solving

If you encounter any problems with installation you can contact us or raise an issue on this repository. Based on user feedback, some common problems and solutions are given below:

• If you're having problems related to permissions, try using sudo before the pip command (e.g., sudo pip install -e .).

Using apptainer

Apptainer is a program for making containers--self contained environments for computing. We supply a file, funmixer_container.def, that builds a container for funmixer, and the file readme_apptainer.md which contains instruction on how to use it. The advantage of apptainer is you can build the container and then copy that file (a single file with the extension .sif) to any computer with apptainer installed (which includes all modern linux systems by default) and run Funmixer immediately.

Testing

To check if installation has happened correctly you can run the synthetic test script:

python3 tests/synthetic_test.py

This script generates a synthetic dataset and recovers the original input. The results are then plotted.

Unit-tests

Formal unit-tests can be run using:

pytest tests/random_networks_test.py

These tests randomly generate sample networks (full R-ary trees and balanced trees) up to 100 nodes in size, with random source concentrations and sub-basin areas drawn from distributions spanning two orders of magnitude. The tests pass if all the inputted upstream source chemistry is recovered to a relative accuracy of 1%.

Below-detection-limit unmixing is tested separately:

pytest tests/bdl_unmixer_test.py

Note that exact recovery is only asserted for uncensored data. Censoring genuinely destroys information - many source concentrations fit a censored observation equally well - so censored behaviour is instead tested through deterministic properties: that censoring at a limit above the true value cannot increase the misfit, that a limit below the true value pulls predictions down to it, and that the problem stays convex and DPP-compliant.

Runtime Benchmark

A timing benchmark can be run using:

python tests/runtime_benchmark.py run

This script benchmarks the runtime of the algorithm for the GUROBI, ECOS and SCS solvers for branching networks up to 500 nodes. This takes ~ 30 minutes to run on standard laptop hardware. The results are cached to file and be visualised using:

python tests/runtime_benchmark.py plot

Usage

Some documented example scripts are given in the directory examples/, and are run from the root directory of the repository, e.g.,

python examples/unmix_mwe.py

examples/unmix_BDL.py is the equivalent example for data containing below-detection-limit values.

Cite

If you use this please cite the paper, which is published at Water Resources Research (a copy is included at docs/Barnes_Lipp_2024.pdf):

Barnes, R., & Lipp, A. G. (2024). Using convex optimization to efficiently apportion tracer and pollutant sources from point concentration observations. Water Resources Research, 60, e2023WR036159. https://doi.org/10.1029/2023WR036159

A .cff citation file is also provided in the repository.