Doc updates for nfdata

Author: samharrison7

docs/_toc.yml

@@ -14,6 +14,7 @@ parts:
     - file: users/input-data
       sections:
         - file: users/nanofase-data
+        - file: users/flow-dir-dem-conditioning
         - file: users/netcdf-namelist-input
         - file: users/parameter-reference
     - file: users/output

docs/intro.md

@@ -8,7 +8,7 @@ The NanoFASE model is a multimedia spatiotemporal environmental fate and exposur
 
 ## Quickstart guide
 
-The easiest way to install the model is to use the [Fortran Package Manager]([https://fpm.fortran-lang.org/en/install/index.html#install](https://fpm.fortran-lang.org/install/index.html#install)), which can be installed using Conda (`conda install -c conda-forge fpm`). The model also requires a recent version of GFortran, NetCDF Fortran and Git installed (see [](./getting-started/building-the-model.md)).
+The easiest way to install the model is to use the [Fortran Package Manager](https://fpm.fortran-lang.org/install/index.html#installing-fpm), which can be installed using Conda (`conda install -c conda-forge fpm`). The model also requires a recent version of GFortran, NetCDF Fortran and Git installed (see [](./getting-started/building-the-model.md)).
 
 Clone the code (make sure to `--recurse-submodules`) and use *fpm* to install:
 

docs/users/flow-dir-dem-conditioning.md

@@ -0,0 +1,49 @@
+# Flow direction and DEM conditioning
+
+The model routes water by using a flow direction raster, which uses standard eight-direction (D8) flow direction enconding. In `nfdata`, this flow direction can either be provided by the `flow_dir` parameter, or `nfdata` can attempt to deduce it from a digital elevation model (DEM) provided by the `dem` parameter.
+
+:::{admonition} `flow_dir` or `dem` as a basis for the grid
+:class: important
+Either the `flow_dir`, or if `flow_dir` isn't present, the `dem` variable, is used as the basis for the grid that defines the NetCDF file output by `nfdata`. That is, parameters such as the CRS, bounds, resolution and mask are obtained from this variable. `nfdata` does not do any reprojection of rasters from other grid systems (therefore, all provided rasters should be the same CRS and resolution), but it does clip rasters to the extent of `flow_dir`/`dem`. This clipping includes masking over nodata values in `flow_dir`/`dem`, and therefore it is one of these variable that decides which grid cells are modelled.
+:::
+
+## Directly providing `flow_dir`
+
+```{figure} ../img/flowdir.gif
+---
+align: right
+---
+```
+
+If the `flow_dir` variable is present in `nfdata`'s config, this will be used as the flow direction (regardless of whether `dem` is present). This variable should be a raster (openable by [Rasterio](https://rasterio.readthedocs.io/en/latest/index.html)) and have an integer data type with either D8 values or nodata values.
+
+This flow direction _must_ be resolved, otherwise running the model with it will fail. By resolved, we mean that all non-masked (nodata) cells must have a non-zero D8 value (1, 2, 4, 8, 16, 32, 64, 128), and that these values must not result in "pools" where cells flow into each other (e.g. neighbouring cells with the northerly cell having `flow_dir=4` and a southerly cell having `flow_dir=64`).
+
+The `flow_dir` group in the config file must contain a `path`, and optionally a `crs` value. The CRS is usually deduced from the provided raster, but in the case of rasters with ill-defined CRSs, the `crs` option gives the ability to override this. Note that this _does not_ reproject the raster to this CRS, it simply assigns the CRS to the raster.
+
+## Calculating `flow_dir` from a DEM
+
+Alternatively, if the `dem` variable is present and `flow_dir` is not, the model will _attempt_ to deduce the flow direction from this DEM. To achieve this, either a _conditioned_ DEM should be provided, or the config option `condition_dem` must be present to instruct `nfdata` to attempt to condition the DEM.
+
+By _conditioned_, we mean that no pits, depressions or flats should be present in the DEM, as these will prevent a fully-resolved flow direction from being calculated. If an unconditioned DEM is provided (and the `condition_dem` option is not specified), `nfdata` will raise an error. There are [numerous](https://mattbartos.com/pysheds/) [tools](https://richdem.readthedocs.io/) [available](https://www.whiteboxgeo.com/manual/wbt_book/available_tools/hydrological_analysis.html) that can help condition DEMs.
+
+If `condition_dem` is present, `nfdata` will attempt to condition the DEM using [Whitebox Tools' `breach_depressions_least_cost` routine](https://www.whiteboxgeo.com/manual/wbt_book/available_tools/hydrological_analysis.html#BreachDepressionsLeastCost). This is a fairly robust tool that is usually able to resolve most DEMs without significant alterations, but it might not work in all cases. `condition_dem` can either be `True`, or be a dict of further config options:
+
+```yaml
+condition_dem:
+  # If `save_dem_to_path` is set, the conditioned DEM will be saved to this path
+  # (even if the conditioning hasn't fully worked). If not, a temporary directory
+  # is used and the resulting DEM is deleted after the processing.
+  save_dem_to_path: ./dem.tif
+  # Additional config options specified here will be passed to `breach_depressions_least_cost`
+  # See https://www.whiteboxgeo.com/manual/wbt_book/available_tools/hydrological_analysis.html#BreachDepressionsLeastCost
+  # If a value for `callback` is provided, it is ignored as we implement our own callback to
+  # deal with errors. The only required (by Whitebox) parameter is `dist`, which specifies
+  # the maximum search distance for breach paths in cells. The default value used in nfdata is 10000.
+  dist: 10000
+```
+
+:::{important}
+Currently, `nfdata`'s flow direction algorithm assumes that nodata (masked) cells and cells immediately outside of the DEM's bounding box are at a very low elevation, therefore forcing the flow direction from neighbouring cells towards them. Effectively, this means that nodata cells and the region outside of the DEM are treated as waterbodies into which neighbouring cell drain. Therefore, we only recommend to use this flow direction calculation functionality if this is a valid assumption for your scenario, otherwise you should provide a pre-resolved flow direction raster.
+:::
+

docs/users/nanofase-data.md

@@ -116,19 +116,40 @@ The examples are annotated and should be self-explanatory. However, there are a
 The `create` and `edit` config files follow a similar layout. A variety of setup data is required and in the examples is placed at the top of the file. This includes file paths to input and output path, and model config info (e.g. timestep info):
 
 ```yaml
-# Setup
-nanomaterial: TiO2                                    # Name of the nanomaterial. Not used in model.
+# Name of the nanomaterial. Not used in the model.
+nanomaterial: TiO2
+# Where do you want the output NetCDF and constants files to be stored?
 output:
-  nc_file: ./data.nc                                  # Where do you want the output NetCDF file to be stored?
-  constants_file: ./constants.nml                     # Where do you want the output constants file to be stored?
-constants_file: ./data.example/thames_tio2_2015/constants.yaml    # Where is the input constants file?
-land_use_config: ./data.example/thames_tio2_2015/land_use.yaml    # Where is the input land use config file?
-root_dir: ./data.example/thames_tio2_2015/            # Root dir, can be used in path variables below as <root_dir>
-iso3: GBR                                             # iso3 code for country modelled
+  nc_file: ./data.nc
+  constants_file: ./constants.nml
+# Where is the input constants file?
+constants_file: ./data.example/thames_tio2_2015/constants.yaml
+# Where is the input land use config file (optional)?
+land_use_config: ./data.example/thames_tio2_2015/land_use.yaml
+# Root directory can be used in path variables below as <root_dir> (optional)
+root_dir: ./data.example/thames_tio2_2015/
+# ISO3 code using for filtering point source emissions to specific countries
+iso3: GBR
+# Number of timesteps, length of each timestep (in seconds) and start date
+# for the model run
 time:
-  n: 365                                              # Number of timesteps to run the model for
-  dt: 86400                                           # Length of each timestep in seconds
-  start_date: 2015-01-01                              # Start date for the model run
+  n: 365
+  dt: 86400
+  start_date: 2015-01-01
+# If the `condition_dem` config option is present, the compiler will attempt to remove
+# pits, depressions and flat areas from the DEM using WhiteboxTools'
+# `breach_depressions_least_cost`. It can either be True or a dict of config options
+condition_dem:
+  # If `save_dem_to_path` is set, the conditioned DEM will be saved to this path
+  # (even if the conditioning hasn't fully worked). If not, a temporary directory
+  # is used and the resulting DEM is deleted after the processing.
+  save_dem_to_path: ./dem.tif
+  # Additional config options specified here will be passed to `breach_depressions_least_cost`
+  # See https://www.whiteboxgeo.com/manual/wbt_book/available_tools/hydrological_analysis.html#BreachDepressionsLeastCost
+  # If a value for `callback` is provided, it is ignored as we implement our own callback to
+  # deal with errors. The only required (by Whitebox) parameter is `dist`, which specifies
+  # the maximum search distance for breach paths in cells. The default value used in nfdata is 10000.
+  dist: 10000
 ```
 
 ```{warning}
@@ -201,6 +222,11 @@ For reference, the NanoFASE land use categories are:
   - `other`
 ```
 
+(nanofase-data:dem-conditioning)
+#### DEM conditioning - `condition_dem`
+
+See [](./flow-dir-dem-conditioning).
+
 
 (nanofase-data:parameters)=
 ### Parameters
@@ -222,6 +248,10 @@ A few parameters require additional information:
 Individual model runs have constant emissions for the whole run. However, individual model runs can be chained together in a multi-year model run, each year having different areal emissions. See [](batch). This is how multi-year simulations with varying areal emissions are currently performed.
 ```
 
+#### Flow direction and DEM
+
+See [](flow-dir-dem-conditioning) for details on the `flow_dir` and `dem` variables.
+
 (nanofase-data:point-emissions)=
 #### Point source emissions and temporal profiles
 

docs/users/netcdf-namelist-input.md

@@ -25,6 +25,17 @@ If you are manually creating this NetCDF file, pay careful attention to units, w
 
 This selection of parameters pertains to the geographical region modelled, such as its climate, topography and soil properties.
 
+(input-data:flow_dir)=
+`int flow_dir(y, x)`
+: *Required, unitless* \
+Flow direction of water in the grid cell, relating to the topography. In `nfdata`, this can be [deduced from the DEM](flow-dir-dem-conditioning). Uses the standard eight-direction (D8) flow direction encoding.
+
+```{figure} ../img/flowdir.gif
+---
+align: right
+---
+```
+
 (input-data:dem)=
 `int dem(y, x)`
 : *Optional, if not present then default slope of 0.0005 m/m is assumed. Units: dm. Standard name: `height_above_mean_sea_level`.* \
@@ -245,17 +256,8 @@ Grid resolution - the size of each grid cell.
 : *Required, units: metres* \
 Coordinates of the bounding box for the model grid, in order (left, top, right, bottom).
 
-The following variables are ultimately deduced from the [digital elevation model](input-data:dem) (DEM) and are used to define the surface water network within the model, meaning the surface water network is automatically deduced and does not need to be provided as an input (see [](conceptual-structure:surface-water-network)).
+The following variables are ultimately deduced from the [digital elevation model](input-data:dem) (DEM) and/or [flow direction](input-data:flow_dir) and are used to define the surface water network within the model, meaning the surface water network is automatically deduced and does not need to be provided as an input (see [](conceptual-structure:surface-water-network)).
 
-```{figure} ../img/flowdir.gif
----
-align: right
----
-```
-
-`int flow_dir(y, x)`
-: *Required, unitless* \
-Flow direction of water in the grid cell, relating to the topography and deduced from the [DEM](input-data:dem). Uses the standard eight-direction (D8) flow direction encoding.
 
 `short outflow(y, x, d)`
 : *Required, unitless* \

docs/users/parameters/met-hydro-parameters.csv

@@ -1,5 +1,6 @@
 Name(dimensions),Specified in,Data type,Required or optional,Model units,Description,Potential source
-"`dem(y, x)`",{bdg-primary}`data.nc` {bdg-success}`config.yaml`          ,float,"Required in {bdg-success}`config.yaml`, optional in {bdg-primary}`data.nc`",dm,"Digital elevation model (height above mean sea level). In the model, this is only used to calculate slope, which defaults to 0.0005 m/m. In the NanoFASE data module, it is used to calculate numerous [secondary variables](netcdf-namelist-input:secondary-derived-variables) for input to the model. Standard name: `height_above_mean_sea_level`",[HydroSHEDS](https://www.hydrosheds.org/) or [Copernicus DEM](https://dataspace.copernicus.eu/explore-data/data-collections/copernicus-contributing-missions/collections-description/COP-DEM)
+"`flow_dir(y, x)`",{bdg-primary}`data.nc` {bdg-success}`config.yaml`          ,int,"Optional in {bdg-success}`config.yaml` if `dem` present",-,"Flow direction using D8 indices. `nfdata` will attempt to calculate this from `dem` if this is not provided in {bdg-success}`config.yaml`.,
+"`dem(y, x)`",{bdg-primary}`data.nc` {bdg-success}`config.yaml`          ,float,"Required in {bdg-success}`config.yaml`, optional in {bdg-primary}`data.nc`",dm,"Digital elevation model (height above mean sea level). In the model, this is only used to calculate slope, which defaults to 0.0005 m/m. In the NanoFASE data module, it can be used to calculate flow direction. Standard name: `height_above_mean_sea_level`",[HydroSHEDS](https://www.hydrosheds.org/) or [Copernicus DEM](https://dataspace.copernicus.eu/explore-data/data-collections/copernicus-contributing-missions/collections-description/COP-DEM)
 "`quickflow(t, y, x)`",{bdg-primary}`data.nc` {bdg-success}`config.yaml`          ,float,Required,m/dt,"Also called surface runoff, this is the component of runoff that is responsible for driving soil erosion.  Standard name: `surface_runoff_flux`",Usually obtained from a hydrological model
 "`runoff(t, y, x)`",{bdg-primary}`data.nc` {bdg-success}`config.yaml`          ,float,Required,m/dt,This is the total runoff (subsurface and surface) that is routed by the model to provide surface water flows.  Standard name: `runoff_flux`,Usually obtained from a hydrological model
 "`precip(t, y, x)`",{bdg-primary}`data.nc` {bdg-success}`config.yaml`          ,float,Required,m/dt,"Amount of precipitation, used to drive soil percolation. Standard name: `rainfall_flux`",[Copernicus E-OBS](https://cds.climate.copernicus.eu/cdsapp#!/dataset/10.24381/cds.151d3ec6?tab=overview) or [CEH-GEAR](https://doi.org/10.5285/dbf13dd5-90cd-457a-a986-f2f9dd97e93c) for the UK