User Guide

This user guide provides a quick overview of how to use DYCOVE. Users can get started by reviewing this guide and running the Examples, but advanced use of DYCOVE requires an understanding of the concepts discussed in the Background.

General Steps

  1. Create a VegetationSpecies instance using a vegetation attribute .json file.

  2. Create a HydroSimulationBase instance by combining numerical model inputs/pointers with the VegetationSpecies object. This would be called as either ANUGA or DFM, depending on which numerical model is used.

  3. Run the simulation using run_simulation() by providing simulation time parameters.

  4. Inspect outputs using built-in plotting functions or external tools.

Imports

Import DYCOVE’s public API as follows:

from dycove import VegetationSpecies                      # main vegetation class
from dycove import MultipleVegetationSpecies as MultiVeg  # for combining multiple species
from dycove import ANUGA_hydro                            # ANUGA model interface
from dycove import DFM_hydro                              # Delft3D FM model interface
from dycove import plotting                               # model plotting module

Of course, users would only import the numerical model interface they plan to use.

Defining the Vegetation Species

Creating the vegetation species object is the first step in setting up a DYCOVE simulation. This object contains all of the vegetation attributes needed to run the simulation, read from a user-provided .json file (see the Background for details on the required inputs).

veg = VegetationSpecies("path/to/vegetation_attributes.json")

Creating a VegetationSpecies object automatically creates a VegetationAttributes instance, which contains all of the vegetation attributes read from the input file. Following the first colonization of the simulation, a VegCohort object is also created.

Creating the Model Interface (Delft3D FM)

Instantiate the model interface by calling DFM. This class inherits HydroSimulationBase, which actually runs the simulation, and instantiates the Delft3D FM “engine” DFMEngine.

model = DFM_hydro.DFM("path/to/system/files/", 'dimr_config.xml', 'model.mdu', vegetation=veg)

The first argument dfm_path points to the Delft3D FM system files (.dll library). The directory must contains dflowfm/ and dimr/ directories. On Windows machines, this directory might look like this:

dfm_path = 'C:\\Program Files (x86)\\Deltares\\Delft3D Flexible Mesh Suite HM (2021.03)\\plugins\\DeltaShell.Dimr\\kernels\\x64'

The second argument dimr_config is the path to the DIMR configuration file, which will be located in your working directory, after exporting the model as DIMR from the GUI (see tide channel example). The third argument mdu_file is the name of the Delft3D FM .mdu file.

Creating the Model Interface (ANUGA)

Instantiate the model interface by calling ANUGA. This class inherits HydroSimulationBase, which actually runs the simulation, and instantiates the ANUGA “engine” AnugaEngine.

model = ANUGA_hydro.ANUGA(anuga_domain, vegetation=veg)

The anuga_domain object is the domain instance created for ANUGA simulations, after all aspects like elevation, friction, forcings, and boundary conditions have been assigned to it.

Multiple Vegetation Species

Users have the option to model multiple species together, using MultipleVegetationSpecies. A separate .json file and VegetationSpecies object are needed for that species, which can be added to model like so (using the ANUGA engine as an example):

veg_alt = VegetationSpecies("path/to/vegetation_attributes_alt.json")
model = ANUGA_hydro.ANUGA(anuga_domain, vegetation=MultiVeg[veg, veg_alt])

Running the Simulation

Users can run a simulation using the following method of HydroSimulationBase, with only the total desired simulation time as an input:

model.run_simulation(3, sim_time_unit="eco-morphodynamic years")

Note that there are several optional arguments of run_simulation() that define the time scaling of DYCOVE models, which are described in the API and in greater detail in the Background.

Outputs

Numerical model outputs can be accessed and inspected as they typically would be without vegetation. For Delft3D FM, the _his.nc file for point output and _map.nc file for map output are found in the /dflowfm/output directory. For ANUGA, the output .sww file is located in the working directory.

DYCOVE vegetation outputs are found in /dflowfm/veg_output for Delft3D FM and in /veg_output for ANUGA. DYCOVE output consists of a series of .nc (netCDF) files, one file per vegetation cohort per ecological time step. One .nc file represents the attribute state of a single vegetation cohort at a point in time. It contains all of the attributes that are stored in a VegCohort object (fractions, stem heights, mortality fractions, etc.). It also contains metadata with the current ecological year and Ecological Time Step (ETS), species name, and cohort ID number.

Plotting

DYCOVE features a built-in plotting class ModelPlotter that can be used to plot time series of 2-D quantities from both the numerical model and DYCOVE:

from pathlib import Path

plotter = plotting.ModelPlotter(
    simdir = Path("."),
    # or: quantity = "Velocity", etc. See examples and ModelPlotter API for list of quantities
    quantity = "Stem Height",
    plot_times = {
        "plotHR_0": 0*24.,
        "plotHR_f": 21*24.,
        "mapHR_int": 1,
        "plotHR_int": 1,
    },
)

plotter.run()

Attributes in the .nc file that are 1-D arrays (fractions, mortality fractions) use the same indexing as the underlying mesh, and can be interpolated to 2-D grids using the 1-D X and Y arrays in the numerical model output files. Users can create an interpolation function directly using create_nn_interpFunc(). This function is created and called as shown in get_quantity_grids(). Loading of vegetation output files is performed in BaseMapLoader, while loading of numerical model outputs occurs in ANUGAMapLoader or DFMMapLoader.

ModelPlotter is customizable and has many optional arguments, which are detailed in the API docs.

Plotting Results with Multiple Species

Plotting results of simulations with multiple species using ModelPlotter is no different from plotting results of a single-species simulation. Plotting stem height, for example, will average all stem heights in each grid cell among the various species and life stages present within that cell. When plotting vegetation fractions or mortality, however, ModelPlotter will create individual plots for each species. Example: A simulation of one ecological year and two species consists of, say, one colonization event for each species in the second ETS of the year. Plots of stem heights will consist of one plot per ETS, starting at ETS 2. Plots of fractions will consist of two plots per ETS, one for each species, starting at ETS 2.

Accessing and Plotting Outputs Without ModelPlotter

The following explanations and examples are for users who prefer to develop their own plotting codes. ANUGA model quantities can be read from the output .sww file using the xarray library (also netCDF4):

import xarray as xr

map_vars = xr.load_dataset("path/to/anuga_model.sww")

# Vertex quantities
x = map_vars["x"]
y = map_vars["y"]
z = map_vars["elevation"]

# Time-varying vertex quantities
stage = map_vars["stage"]
xmom = map_vars["xmomentum"]
depth = stage - z

# Some quantities are available directly as centroids
z_c = map_vars["elevation_c"]
stage_c = map_vars["stage_c"]
xmom_c = map_vars["xmomentum_c"]

These are 1-D arrays, that can be interpolated using create_nn_interpFunc() or other tools. Note that stage and other time-varying quantities have a time axis as the first axis (same with DFM). Centroid data is also included in the .sww file, but x and y values at centroids would need to be computed as they are in dycove.utils.model_loader.ANUGAMapLoader._load_outputs (or anuga.utilities.plot_utils._get_centroid_values). Note that DYCOVE vegetation arrays are based on centroids.

Similarly, DFM model centroid quantities can be read from the output _map.nc file:

import xarray as xr

map_vars = xr.load_dataset("path/to/dflowfm/output/FlowFM_map.nc")
x_c = map_vars["mesh2d_face_x"]
y_c = map_vars["mesh2d_face_y"]
z_c = map_vars["mesh2d_flowelem_bl"]

# Time-varying quantities
bedlevel = map_vars["mesh2d_mor_bl"]
stage = map_vars["mesh2d_s1"]
velocity = map_vars["mesh2d_ucmag"]
depth = stage - bedlevel

Vegetation quantities can be accessed via the files in the veg_output directory. Quantities that are saved in these output files are listed as attributes of the VegCohort class. An example of how to load and plot quantities (outside of ModelPlotter) is provided below. This example is based on ANUGA, but DFM users can modify it with just a few changes based on reading files as instructed above.

Vegetation fractions and cohorts are tracked individually in DYCOVE, and we can choose to plot quantities individually for each species or take a weighted average in each grid cell, with weights determined by the fractions present in the cell. In this example, we plot a weighted average stem height for each cell. The option to plot species results separately can be implemented via the plot_separate_species flag in ModelPlotter. Either way, we need to read in the underlying hydrodynamic model file for interpolating to a regular grid:

from pathlib import Path
import matplotlib.pyplot as plt
import json
import xarray as xr
import numpy as np
from dycove.utils.array_math import cell_averaging
from dycove.utils.plotting import create_nn_interpFunc
from dycove.utils.model_loader import get_anuga_centroid_coords

# Declare model directory paths
model_dir = Path("path/to/anuga_model_dir")
eco_dir = model_dir / "veg_output"

# Load ANUGA (or DFM) output, read X and Y coordinate arrays
map_vars = xr.load_dataset(model_dir / "rectang_beach.sww")

# Convert ANUGA vertex coordinates to centroids (may take a little time)
# DFM USERS: x_c and y_c can be pulled from the output file directly (see above)
x_c, y_c = get_anuga_centroid_coords(map_vars)

# Create interpolation function for 10-m grid. Can also provide as an argument a polygon .csv file to mask outside domain.
interp_func = create_nn_interpFunc(x_c, y_c, 10)

# Interpolate model bathymetry to use as a base map (using centroids)
z_grid = interp_func(map_vars["elevation_c"])

# Read output cohort index file, categorizing each cohort output file by ecological year and ETS
with open(eco_dir / "_cohort_files_ets_index.json", "r") as f:
    cohort_index = json.load(f)

# Loop through all saved cohort files by year and ETS, load and plot data
eco_years = sorted(cohort_index, key=int)  # output eco years in order
for year in eco_years:
    ets_list = sorted(cohort_index[year], key=int)  # ETS in order for current eco year
    for ets in ets_list:
        fractions, stem_heights = [], []
        for file in cohort_index[year][ets]:
            c = xr.load_dataset(eco_dir / (file + ".nc"))

            # Append to list all data from this ETS
            fractions.append(c["fraction"])  # each c["fraction"] is an array
            stem_heights.append(c.attrs["height"])  # each c.attrs["height"] is a scalar

        # Do weighted average based on vegetation fractions in each cell
        stemht_avg = cell_averaging(fractions, stem_heights)

        # Interpolate to grid using same interp_func as for the model elevation values
        stemht_grid = interp_func(stemht_avg)

        # Mask out non-vegetated cells
        stemht_grid = np.ma.masked_where(stemht_grid < 0.05, stemht_grid)

        # Do plotting (colorbars, etc can be added as well)
        fig, ax = plt.subplots()
        im_base = ax.imshow(z_grid, cmap="Greys_r")
        im_veg = ax.imshow(stemht_grid, cmap="Greens", vmin=0, vmax=1)
        ax.set_title(f"Stem height [m] -- eco year {year}, ETS {ets}")
        plt.show()
        plt.close()

For each cohort output file, scalar quantities like 'eco_year', 'height', 'diameter', and 'density' are stored as metadata under attrs. These can be accessed with c.attrs['eco_year'], etc. Array quantities can be accessed directly from the object, like c['fraction'].