# Hydrology Toolbox Guide

Tools around the hydrology model of the Generic PEM (`do_hydrology=.true.`):

- build the hydrologic database (depression hierarchy, lake volume / area / elevation
  curves) from a global `.tif` DEM;
- display the hydrology state `starthydro.nc` written by the PEM.

## Contents

| File | Purpose |
|---|---|
| `database.py` | Builds the hydrologic database from a DEM: runs `dephier.exe`, stitches the map east-west, rebuilds the depression tree, then runs `database` |
| `database.F90` | Source of the Fortran `database` executable: lake volumes, areas and elevations of every depression, written in the final database |
| `display_starthydro.py` | Maps and figures of a `starthydro.nc` file against its database: water depth/elevation, outflow, streams, basins, drainage systems, comparison with observations |

## Build the Hydrologic Database

The workflow uses three programs:

| Step | Program | Work | Writes |
|---|---|---|---|
| 1 | `dephier.exe` (C++, Barnes et al. 2020) | Map of the leaf depressions and depression hierarchy graph | `dephier-label.tif`, `dephier-graph-detail.csv` |
| 2 | `database.py` | East-west stitching of the map, rebuild of the depression tree; then runs step 3 | `hydrologic_database.nc` |
| 3 | `database` (Fortran, compiled from `database.F90`) | Lake volumes, areas and elevations | `hydrologic_database_final.nc` |

The Python script `database.py` runs `dephier.exe` (step 1) and the Fortran executable
(step 3) itself, so both executables must exist **before** running the Python script. Their
paths are set by `DEPHIER_EXE` and `DATABASE_EXE` (see section 3 below).

### 1. Build dephier.exe

The source code is available here: <https://github.com/r-barnes/Barnes2019-DepressionHierarchy>.

Requirements: `cmake` ≥ 3.9, a C++17 compiler and GDAL (`gdal`, `libgdal.so`).

```bash
git clone --recurse-submodules -j8 https://github.com/r-barnes/Barnes2019-DepressionHierarchy.git
cd Barnes2019-DepressionHierarchy/src
sed -i 's/std::ofstream fgraph(out_name+"-graph.csv");/std::ofstream fgraph(out_name+"-graph-detail.csv");/g' main.cpp
sed -i 's/fgraph<<"parent,child,oceanlink/fgraph<<"child,parent,dep_vol,pit_cell,out_cell,out_elev,odep,geolink,oceanlink/g' main.cpp
sed -i '59d' main.cpp
sed -i '59ifgraph<<i<<","<<deps[i].parent<<","<<deps[i].dep_vol<<","<<deps[i].pit_cell<<","<<deps[i].out_cell<<","<<deps[i].out_elev<<","<<deps[i].odep<<","<<deps[i].geolink<<",";\' main.cpp
cd ..
rm -rf build && mkdir build && cd build
cmake -DCMAKE_BUILD_TYPE=Release -DUSE_GDAL=On ..
make -j 4
cp dephier.exe ../../dephier.exe
```

Quick test: `./dephier.exe` without argument must print
`Syntax: ... <Input> <Output Prefix> <Ocean Level>`.

### 2. Compile database.F90

The Python script calls `database`, which is compiled from `database.F90`.

Requirements (versions used on Spirit):

| Tool | Version | Note |
|---|---|---|
| gfortran | 9.4.0 | with OpenMP (`-fopenmp`, `libgomp`) |
| NetCDF-C | 4.7 | `libnetcdf` |
| NetCDF-Fortran | 4.5 | `libnetcdff` |

```bash
gfortran database.F90 -I/usr/include -L/usr/lib -lnetcdf -lnetcdff \
    -ffree-line-length-0 -O3 -march=native -fopenmp -fno-range-check -o database
```

### 3. Run database.py

#### Parameters

Set them in the `PARAMETERS` block at the top of `database.py`:

| Parameter | Meaning |
|---|---|
| `TIF_FILE` | Input DEM (`.tif`, latitude -90/90, longitude -180/180) |
| `RESULTS_FOLDER` | Output folder; files are written in `RESULTS_FOLDER/Database_x<RESHAPE_VALUE>/` |
| `R` | Planet radius (km). Written as the global attribute `R` of `hydrologic_database.nc` and used by the Fortran code for the cell areas (lake areas and volumes). If the attribute is missing (older files), the Fortran code uses 3389.5 km |
| `RESHAPE_VALUE` | DEM degradation factor (minimum over each block): `1` = full resolution (128 px/°), `128` = 1 px/° |
| `ADD_VALUE` | Width of the bands added on the east and west sides to stitch the map: `width / ADD_VALUE` columns on each side |
| `DEPHIER_EXE` | Path of the `dephier.exe` executable (section 1) |
| `DATABASE_EXE` | Path of the Fortran `database` executable (section 2) |

#### Python environment of database.py

The conda environment needs the following packages:

| Package | Version |
|---|---|
| python | 3.8.10 |
| numpy | 1.24.3 |
| pandas | 2.0.3 |
| netCDF4 | 1.6.2 |
| matplotlib | 3.7.2 |
| psutil | 5.9.0 |
| rasterio | 1.3.9 (GDAL 3.6.4) |
| pysheds | 0.4 |

Optional: `imageio` (`build_gif_dephier`) and `treelib` (`build_dephier_tree`), only for
functions that are not called by default.

To create an equivalent environment:

```bash
conda create -n hdb python=3.8 numpy=1.24 pandas=2.0 netcdf4 matplotlib psutil
conda activate hdb
pip install rasterio==1.3.9 pysheds==0.4
```

#### Running on Spirit (Slurm)

Example script (`database.sh`) on the Spirit cluster:

```bash
#!/bin/bash
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --cpus-per-task=8
#SBATCH --partition=zen16
#SBATCH --mem=100G
#SBATCH --time=6:00:00
#SBATCH --output database.out

source your_conda_env_path/conda.sh
conda activate hdb

export OMP_NUM_THREADS=32
ulimit -s unlimited

srun --threads-per-core=1 --cpu-bind=threads -- python -u database.py
```

```bash
sbatch database.sh
```

### Output Files

In `RESULTS_FOLDER/Database_x<RESHAPE_VALUE>/`:

| File | Content |
|---|---|
| `elevation.tif` | DEM (degraded if `RESHAPE_VALUE > 1`) |
| `elevation_dephier.tif` | DEM with the added bands (dephier input) |
| `dephier-label.tif`, `dephier-graph-detail.csv` | Raw dephier outputs |
| `hydrologic_database.nc` | Intermediate database written by the Python script (Fortran input) |
| **`hydrologic_database_final.nc`** | **Final database** (`label` map, `elevation`, depression table, `lake_vol` / `lake_area` / `lake_elev` for 11 filling levels; global attribute `R` = radius used) |

### Use the Database in the PEM

The PEM reads the final database from `hydrodata_path` in `run_pem.def` (default
`hydrological_database.nc` in the simulation directory), so either point `hydrodata_path`
to `hydrologic_database_final.nc` or copy it under the default name:

```bash
cp RESULTS_FOLDER/Database_x<RESHAPE_VALUE>/hydrologic_database_final.nc <simulation_dir>/hydrological_database.nc
```

At initialization, the PEM checks that the planet area of the database matches the one of
its grid: build the database with the radius `R` of the planet that is simulated.

## Display the Hydrology State

`display_starthydro.py` draws the hydrology state saved by the PEM, together with the
database it was computed on (the same file as `hydrodata_path`).

### Inputs

- `--file`: a hydrology state written by the PEM: `starthydro.nc` in the simulation
  directory, `starts/starthydro<N>.nc` at the end of the PEM run `N`, or
  `starts/starthydro<N>_ts<i>.nc` at the iteration `i` of that run when the overflow is
  computed (`overflow_compute_rate` in `run_pem.def`);
- `--database-path`: the hydrologic database (`hydrologic_database_final.nc`);
- `--output-path`: the output folder, where `tmp/` (intermediate maps, reused with
  `--load`) and `figure/` are created;
- `--obs-path` (optional): the folder of the observation data, `toolbox/data` by default.
  It is only needed by the actions comparing with observations.

The overflow, so the outflow of the lakes, is only computed every `overflow_compute_rate`
PEM iterations: the actions built on the outflow (`water-outflow`, `stream`,
`matching-streams`, `global-streams`, `basin-type`) need one of the `_ts<i>` states.

### Python environment of display_starthydro.py

Required: `numpy`, `pandas`, `netCDF4`, `matplotlib`, `cartopy`, `rasterio`, `geopandas`,
`shapely`, `pyproj`, `pyshp` (`shapefile`) and `whitebox`.

Optional: `pysheds` and `affine` for the stream actions (`stream`, `matching-streams`,
`global-streams`), `vedo` for the 3D view of the water depth.

### Usage

```bash
python display_starthydro.py --file starts/starthydro1.nc \
    --database-path hydrological_database.nc --output-path hydro_figs \
    --action water-depth --save
```

Every run needs `--file`, `--database-path`, `--output-path` and `--action`. Run with `-h`
for the full list of options and their defaults.

| Action | Draws |
|---|---|
| `hydro-function`, `hydro-functions` | Lake volume and area against the water elevation, for the first depression or for several of them |
| `database` | One variable of the database (`--database-variable`, `--database-type`) |
| `watershed-distrib` | Mean watershed area of the elementary depressions per degree of latitude and of longitude of their outlet |
| `aw` | Watershed areas of one zoom region |
| `water-volume` | Flooded cells, where the water level of an active depression lies above the ground |
| `water-depth`, `water-elevation`, `watershed-area` | Maps of the water depth, water surface elevation and watershed area |
| `water-outflow` | Map of the outflow of the lakes (`--lake-only`, `--noflow`, `--zoom-only`) |
| `basin-type` | Open and closed (endorheic) lakes |
| `drainage-system`, `drainage-size` | Drainage systems and their size |
| `global-streams` | Global stream network (`--coarse`, `--min-area`, `--color-by`, `--label-systems`) |
| `stream` | Streams of one zoom region |
| `matching-streams` | Simulated streams of one zoom region against the valley networks of Hynek et al. (2010) |
| `raster` | One raster (`--raster`: elevation, water depth, water elevation, watershed area or label) |
| `observations`, `observed-data` | Load the observation layers (`--observations`) and draw them on the topography |
| `compare-sim-obs` | Score the state against the open-/closed-basin lakes of Goudge et al. (2016) and the deltas of Di Achille and Hynek (2010) |
| `check-depressions` | Number of active depressions owning each cell; 2 or more flags nested active depressions, which should never happen |
| `clear-tmp` | Remove the intermediate files of `tmp/` |

Common options: `--zoom` (`nirgal`, `jezero`, `gale`, `subregion` or `all`), `--save` and
`--save-name`, `--load` to reuse the maps of `tmp/`, `--shadow` (topographic shading),
`--axis`, `--vmin`/`--vmax`, `--cmap`, `--tex` (LaTeX fonts).

With `--save`, the figures are written without opening any window, which also works on a
node without a display.

### Observation data

Expected under `--obs-path`, at these places or anywhere below it by name:

| Layer | File |
|---|---|
| `basin` | `goudge/Goudge_et_al_2016_Geology_CBL_OBL_Catalogue.csv` |
| `delta` | `deltas.txt` |
| `river` | `hynek/Bigmaprejoin.shp` |
| `shoreline` | `shoreline/Sholes2021_Data/Elevation Data/Arabia/1_Parker1993_Arabia_Topo.csv`, `.../Arabia/7_Perron2007_Arabia_Topo.csv`, `.../Deuteronilus/0_Ivanov2017_Deuteronilus_Original_Topo.csv` |

These data are not versioned.

## Known Limitations

- `database.py`: the parameters are set in the script itself, not on the command line
- `database.py`: `dephier.exe` must be built from its sources with the patches of section 1
- `display_starthydro.py`: the zoom regions and the observation data are those of Mars
- `display_starthydro.py`: at full resolution (128 px/°) one whole-DEM map takes 4.25 GB; the maps are written to `tmp/` when the memory limit of the process is reached

## Related Docs

- `../README.md`: toolbox index
- `../../USER_GUIDE.md`: build, run and `--init-hydro-pcm`
- `../../deftank/generic/run_pem.def`: hydrology parameters (`do_hydrology`, `hydrodata_path`, `gel`, `overflow_compute_rate`)