.. _user-guide: 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 :ref:`Examples `, but advanced use of DYCOVE requires an understanding of the concepts discussed in the :ref:`Background `. General Steps ------------- 1. Create a :class:`~dycove.sim.vegetation.VegetationSpecies` instance using a vegetation attribute `.json` file. 2. Create a :class:`~dycove.sim.base.HydroSimulationBase` instance by combining numerical model inputs/pointers with the :class:`~dycove.sim.vegetation.VegetationSpecies` object. This would be called as either :class:`~dycove.sim.engines.ANUGA_hydro.ANUGA` or :class:`~dycove.sim.engines.DFM_hydro.DFM`, depending on which numerical model is used. 3. Run the simulation using :meth:`~dycove.sim.base.HydroSimulationBase.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: .. code-block:: python 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 :ref:`Background ` for details on the required inputs). .. code-block:: python veg = VegetationSpecies("path/to/vegetation_attributes.json") Creating a :class:`~dycove.sim.vegetation.VegetationSpecies` object automatically creates a :class:`~dycove.sim.vegetation_data.VegetationAttributes` instance, which contains all of the vegetation attributes read from the input file. Following the first colonization of the simulation, a :class:`~dycove.sim.vegetation_data.VegCohort` object is also created. Creating the Model Interface (Delft3D FM) ----------------------------------------- Instantiate the model interface by calling :class:`~dycove.sim.engines.DFM_hydro.DFM`. This class inherits :class:`~dycove.sim.base.HydroSimulationBase`, which actually runs the simulation, and instantiates the Delft3D FM "engine" :class:`~dycove.sim.engines.DFM_hydro.DFMEngine`. .. code-block:: python 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: .. code-block:: python 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 :ref:`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 :class:`~dycove.sim.engines.ANUGA_hydro.ANUGA`. This class inherits :class:`~dycove.sim.base.HydroSimulationBase`, which actually runs the simulation, and instantiates the ANUGA "engine" :class:`~dycove.sim.engines.ANUGA_hydro.AnugaEngine`. .. code-block:: python 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 :class:`~dycove.sim.vegetation.MultipleVegetationSpecies`. A separate `.json` file and :class:`~dycove.sim.vegetation.VegetationSpecies` object are needed for that species, which can be added to model like so (using the ANUGA engine as an example): .. code-block:: python 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 :class:`~dycove.sim.base.HydroSimulationBase`, with only the total desired simulation time as an input: .. code-block:: python model.run_simulation(3, sim_time_unit="eco-morphodynamic years") Note that there are several optional arguments of :meth:`~dycove.sim.base.HydroSimulationBase.run_simulation` that define the time scaling of DYCOVE models, which are described in the API and in greater detail in the :ref:`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 :class:`~dycove.sim.vegetation_data.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 :class:`~dycove.utils.plotting.ModelPlotter` that can be used to plot time series of 2-D quantities from both the numerical model and DYCOVE: .. code-block:: python 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 :func:`~dycove.utils.plotting.create_nn_interpFunc`. This function is created and called as shown in :meth:`~dycove.utils.plotting.ModelPlotter.get_quantity_grids`. Loading of vegetation output files is performed in :class:`~dycove.utils.model_loader.BaseMapLoader`, while loading of numerical model outputs occurs in :class:`~dycove.utils.model_loader.ANUGAMapLoader` or :class:`~dycove.utils.model_loader.DFMMapLoader`. :class:`~dycove.utils.plotting.ModelPlotter` is customizable and has many optional arguments, which are detailed in the API docs. .. _plot-multi-species: Plotting Results with Multiple Species -------------------------------------- Plotting results of simulations with multiple species using :class:`~dycove.utils.plotting.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, :class:`~dycove.utils.plotting.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``): .. code-block:: python 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 :func:`~dycove.utils.plotting.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: .. code-block:: python 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 :class:`~dycove.sim.vegetation_data.VegCohort` class. An example of how to load and plot quantities (outside of :class:`~dycove.utils.plotting.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 :class:`~dycove.utils.plotting.ModelPlotter`. Either way, we need to read in the underlying hydrodynamic model file for interpolating to a regular grid: .. code-block:: python 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']``.