Using the Virtual Ecosystem#

This page provides a brief demonstration of the Virtual Ecosystem model in operation. Once you have installed the Virtual Ecosystem, you should be able to replicate this example on your own computer using the commands below.

Example model data#

The demonstration requires an installation of the example data provided with the Virtual Ecosystem package. If you have previously attempted to run this example then the simulation will refuse to overwrite existing output files. You can either:

  • delete the existing example data folder and reinstall it,

  • create a fresh installation using a different location, or

  • create and use a new output directory with the existing example data folder.

It is worth re-reading the example data page to get an overview of the directory structure and the configuration and data files.

The ve_run command#

You’ve already used this command to install the example data but most of the options to the ve_run command are used to run the simulation. The --help option can be used to show the various arguments that can be used to set how a model runs:

ve_run --help

ve_run --help

ve_run --help

usage: ve_run [-h] [--version] [--install-example INSTALL_EXAMPLE]
              [-o OUTPATH] [-c CLI_CONFIG] [--validate-config-only]
              [-p CLI_PATHS] [--logfile LOGFILE] [-q]
              [cfg_paths ...]

Configure and run a Virtual Ecosystem simulation.

This program sets up and runs a Virtual Ecosystem simulation. The program expects
to be provided with paths to TOML formatted configuration files for the simulation.
The configuration is modular: a directory path can be used to add all TOML
configuration files in the directory, or individual file paths can be used to select
specific combinations of configuration files. These are combined and validated and
then used to initialise and run the model.

As an alternative to providing configuration paths, the `--install-example` option
allows users to provide a location where a simple example set of datasets and
configuration files provided with the Virtual Ecosystem package can be installed.
This option will create a `ve_example` directory in the location, and users can
examine the input files and run the simulation from that directory:

`ve_run /provided/install/path/ve_example`

The output directory for simulation results is typically set in the configuration
files, but can be overwritten using the `--outpath` option. A log file path can be
provided for logging output. If this is not provided then the log will be written to
the console, but the logging is typically verbose and it is usually better to
redirect the log to a file.

When logging is redirected to a file, a short progress report is written to stdout.
By default, the command reports: the start and end of the simulation and log
location; the completion of simulation stages; and a progress bar over the time
steps of the model. The `--quiet` command can be used to incrementally mute this
output: `-q` will remove the progress bar, `-qq` just prints the start and stop and
`-qqq` mutes the report entirely.

The `--config` option can be used to override configuration settings provided in the
file or to add additional settings. This is typically used to run a set of parallel
simulations that vary configuration settings of interest around a central
configuration setup, without the need to write a specific configuration file for
each permutation.

The `--data-path` option can be used to dynamically set the location of data paths
in the configuration. A file path in the config can be set as a path marker, which
must be a string starting with a "$", for example "$CLIMATE_DATA". This option can
then be used to substitute different files into that marker for different runs:
`--data-path CLIMATE_DATA=/path/to/file.nc`.

The `--validate-config-only` flag can be used to only run the configuration
validation part of the model setup and the exit before running any models.

The resolved complete configuration will then be written to a single consolidated
config file in the output path with a default name of
`ve_full_model_configuration.toml`. This can be disabled by setting the
`core.data_output_options.save_merged_config` option to false. Note that the merged
configuration automatically converts all file paths within the merged configurations
to absolute file paths - this ties the merged configuration to the file system where
the run is executed.

positional arguments:
  cfg_paths             Paths to config files

options:
  -h, --help            show this help message and exit
  --version             show program's version number and exit
  --install-example INSTALL_EXAMPLE
                        Install the Virtual Ecosystem example data to the
                        given location
  -o, --outpath OUTPATH
                        Path for output files
  -c, --config CLI_CONFIG
                        Override configuration settings
  --validate-config-only
                        Exit after validating configuration
  -p, --data-path CLI_PATHS
                        Set data paths used for input data
  --logfile LOGFILE     A file path to use for logging a Virtual Ecosystem
                        simulation
  -q, --quiet           Quieten the default progress reporting

Running the example model#

The code below runs a simulation using the example data. The command uses the command line options to set three things:

  1. It points to the files in config directory that should be used configure the model. N.B. this directory contains configuration for all possible model combinations so providing just the config directory as a path will result in an invalid configuration.

  2. It sets the output directory to be used by the simulation to out. You could create a new output directory (e.g. out_test_2) and change this to run a new simulation using the existing data.

  3. It redirects the model logging to a file in the output directory, rather than printing it all to screen.

When the detailed logging is redirected to a file, the command generates a short progress report to show the model running. This can be made shorter or completely muted by using the -q argument: repeat the argument to remove more details (e.g. -qq or -qqq).

Warning

If the path provided for the log points to a file that already exists the detailed logging is added to the end of the file rather than creating a new file. We would recommend creating a new logfile for each simulation as reusing files in this way can create confusion.

In the example code below, the ve_example folder has been previously installed under the directory /tmp/

ve_run /tmp/ve_example/config/data_config.toml \
    /tmp/ve_example/config/abiotic_simple_config.toml \
    /tmp/ve_example/config/animal_config.toml \
    /tmp/ve_example/config/hydrology_config.toml \
    /tmp/ve_example/config/litter_config.toml \
    /tmp/ve_example/config/plant_config.toml \
    /tmp/ve_example/config/soil_config.toml \
    --out /tmp/ve_example/out \
    --logfile /tmp/ve_example/out/logfile.log

ve_run C:\tmp\ve_example\config\data_config.toml ^
    C:\tmp\ve_example\config\abiotic_simple_config.toml ^
    C:\tmp\ve_example\config\animal_config.toml ^
    C:\tmp\ve_example\config\hydrology_config.toml ^
    C:\tmp\ve_example\config\litter_config.toml ^
    C:\tmp\ve_example\config\plant_config.toml ^
    C:\tmp\ve_example\config\soil_config.toml ^
    --out C:\tmp\ve_example\out ^
    --logfile C:\tmp\ve_example\logfile.log

ve_run C:\tmp\ve_example\config\data_config.toml `
    C:\tmp\ve_example\config\abiotic_simple_config.toml `
    C:\tmp\ve_example\config\animal_config.toml `
    C:\tmp\ve_example\config\hydrology_config.toml `
    C:\tmp\ve_example\config\litter_config.toml `
    C:\tmp\ve_example\config\plant_config.toml `
    C:\tmp\ve_example\config\soil_config.toml `
    --out C:\tmp\ve_example\out `
    --logfile C:\tmp\ve_example\logfile.log

* Starting Virtual Ecosystem simulation using v0.2.0.
* Logging to: ve_example/out/logfile.log
* Loading configuration
* Configuration validated
* Saved compiled configuration: ve_example/out/compiled_configuration.toml
* Built core model components
* Initial data loaded
* Models initialised: abiotic_simple, animal, hydrology, litter, plants, soil
* Initialisation data export complete.
* Starting simulation
100%|██████████████████████████████████████████| 24/24 [01:05<00:00,  2.73s/it]
* Simulation completed
Virtual Ecosystem run complete.

The log file is very long and shows the step by step process of running the model - it is primarily used for diagnosing problems with the model. You can view a sample of the contents in the dropdown below:

Partial log output
[INFO] - main - ve_run(286) - Using Virtual Ecosystem v0.2.0.
[INFO] - config_builder - _collect_config_paths(427) - Config paths resolve to 7 files
[INFO] - config_builder - _load_config_toml(450) - Config TOML loaded from ve_example/config/data_config.toml
[INFO] - config_builder - _load_config_toml(450) - Config TOML loaded from ve_example/config/abiotic_simple_config.toml
[INFO] - config_builder - _load_config_toml(450) - Config TOML loaded from ve_example/config/animal_config.toml
[INFO] - config_builder - _load_config_toml(450) - Config TOML loaded from ve_example/config/hydrology_config.toml
[INFO] - config_builder - _load_config_toml(450) - Config TOML loaded from ve_example/config/litter_config.toml
[INFO] - config_builder - _load_config_toml(450) - Config TOML loaded from ve_example/config/plant_config.toml
[INFO] - config_builder - _load_config_toml(450) - Config TOML loaded from ve_example/config/soil_config.toml
[INFO] - config_builder - _compile_data(374) - Configuration data compiled.
[INFO] - registry - _register_module(163) - Registering module: virtual_ecosystem.core
[INFO] - registry - _register_module(176) - Configuration class registered for virtual_ecosystem.core
[INFO] - registry - _register_module(163) - Registering module: virtual_ecosystem.models.abiotic_simple
[INFO] - registry - _get_model(237) - Registering model class for virtual_ecosystem.models.abiotic_simple: AbioticSimpleModel
[INFO] - registry - _register_module(176) - Configuration class registered for virtual_ecosystem.models.abiotic_simple
[INFO] - registry - _register_module(163) - Registering module: virtual_ecosystem.models.animal
[INFO] - registry - _get_model(237) - Registering model class for virtual_ecosystem.models.animal: AnimalModel
[INFO] - registry - _register_module(176) - Configuration class registered for virtual_ecosystem.models.animal
[INFO] - registry - _register_module(163) - Registering module: virtual_ecosystem.models.hydrology
[INFO] - registry - _get_model(237) - Registering model class for virtual_ecosystem.models.hydrology: HydrologyModel
--- many lines omitted ---
[INFO] - data - __setitem__(237) - Replacing data array for 'plant_pft_propagules'
[INFO] - plants_model - update_canopy_layers(864) - Updated canopy data on 1
[INFO] - data - __setitem__(237) - Replacing data array for 'shortwave_absorption'
[INFO] - data - __setitem__(237) - Replacing data array for 'plant_ammonium_uptake'
[INFO] - data - __setitem__(237) - Replacing data array for 'plant_nitrate_uptake'
[INFO] - data - __setitem__(237) - Replacing data array for 'plant_phosphorus_uptake'
[INFO] - data - __setitem__(237) - Replacing data array for 'stem_lignin'
[INFO] - data - __setitem__(237) - Replacing data array for 'senesced_leaf_lignin'
[INFO] - data - __setitem__(237) - Replacing data array for 'plant_reproductive_tissue_lignin'
[INFO] - data - __setitem__(237) - Replacing data array for 'root_lignin'
[INFO] - data - __setitem__(237) - Replacing data array for 'subcanopy_vegetation_cnp'
[INFO] - data - __setitem__(237) - Replacing data array for 'subcanopy_seedbank_cnp'
[INFO] - data - __setitem__(237) - Replacing data array for 'subcanopy_vegetation_litter_cnp'
[INFO] - data - __setitem__(237) - Replacing data array for 'subcanopy_seedbank_litter_cnp'
[INFO] - data - __setitem__(237) - Replacing data array for 'subcanopy_vegetation_litter_lignin'
[INFO] - data - __setitem__(237) - Replacing data array for 'subcanopy_seedbank_litter_lignin'
[INFO] - data - __setitem__(237) - Replacing data array for 'subcanopy_ammonium_uptake'
[INFO] - data - __setitem__(237) - Replacing data array for 'subcanopy_nitrate_uptake'
[INFO] - data - __setitem__(237) - Replacing data array for 'subcanopy_phosphorus_uptake'
[INFO] - base_model - update(439) - Updating animal model

Looking at the results#

The Virtual Ecosystem writes data out to a Zarr data store. This is a relatively new open-source format that is designed for use with large multi-dimensional data. It is used in the Virtual Ecosystem because the format supports appending data at each time step, making it easy to build up a single integrated dataset during the model run. This also means that if a model crashes, then the output variables will be structured in a single datastore as time series up until the crash point.

The default output filename is model_data.zarr but this can be changed in the output configuration. Although the file suffix makes that look like a single file, the Zarr format is actually a directory structure, containing variables or groups of variables. The Virtual Ecosystem outputs data into one of three groups, depending on when the values are calculated.

  • The inputs group: variables that loaded at the start of the simulation from the model data configuration. Obviously, these data are in your input files, but it can be convenient to have them packaged within a single data source alongside the model outputs.

  • The init group: variables that are calculated during the initialisation process of science models are written to this group. These values capture the state of the model before the model starts to iterate through time. You might not plot them as part of a time series but they can be very useful for understanding how your input data is used to set the model running.

  • The outputs group: variables that are calculated by the science models at each time step. The data are added to the data store at the end of each time step.

For more details on the variables used in the Virtual Ecosystem and to see which variables are part of the inputs, initialisation and outputs, see the variables table.

Note

The configuration for the example model exports all of the variables used in the model. If you only require some of the variables for your analayses, you can alter the output configuration to export a subset of the variables.

The sections below go into more detail on each data group, the code below uses the xarray and matplotlib Python packages to load and visualise output data. You may need to install these to replicate these outputs on your own computer.

import matplotlib.pyplot as plt
import numpy as np
import xarray

Input data#

The inputs group contains variables that users are required to provide to run the model.

  • For some variables - such as elevation and soil pH - these are constant values for each grid cell that are used throughout the simulation.

  • Other variables, - such as precipitation and temperature - provide a time series that is provided which is used to that drive (or force) the behaviour of the model through time.

The code below plots both constant and forcing variables.

# Load the generated data files
inputs = xarray.load_dataset("ve_example/out/model_data.zarr", group="inputs")
inputs

<xarray.Dataset> Size: 175kB
Dimensions:                            (cell_id: 81, time_index: 24,
                                        element: 3, pft: 2)
Coordinates:
  * cell_id                            (cell_id) int64 648B 0 1 2 3 ... 78 79 80
    x                                  (cell_id) int32 324B 0 90 180 ... 630 720
    y                                  (cell_id) int32 324B 720 720 720 ... 0 0
  * time_index                         (time_index) int64 192B 0 1 2 ... 22 23
  * element                            (element) <U1 12B 'C' 'N' 'P'
  * pft                                (pft) <U9 72B 'broadleaf' 'shrub'
Data variables: (12/41)
    air_temperature_ref                (cell_id, time_index) float64 16kB 22....
    atmospheric_co2_ref                (cell_id, time_index) float64 16kB 400...
    atmospheric_pressure_ref           (cell_id, time_index) float64 16kB 100...
    clay_fraction                      (cell_id) float64 648B 0.27 ... 0.27
    downward_longwave_radiation        (cell_id, time_index) float64 16kB 406...
    downward_shortwave_radiation       (cell_id, time_index) int64 16kB 2040 ...
    ...                                 ...
    soil_p_pool_labile                 (cell_id) float64 648B 2.5e-05 ... 2.5...
    soil_p_pool_primary                (cell_id) float64 648B 0.001 ... 0.001
    soil_p_pool_secondary              (cell_id) float64 648B 0.005 ... 0.005
    subcanopy_seedbank_biomass         (cell_id) float64 648B 0.07 0.07 ... 0.07
    subcanopy_vegetation_biomass       (cell_id) float64 648B 0.07 0.07 ... 0.07
    wind_speed_ref                     (cell_id, time_index) float64 16kB 0.1...
extent = [
    float(inputs["x"].min()),
    float(inputs["x"].max()),
    float(inputs["y"].min()),
    float(inputs["y"].max()),
]

# Make two side by side plots
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10, 5))

# Elevation
im1 = ax1.imshow(inputs["elevation"].to_numpy().reshape((9, 9)), extent=extent)
ax1.set_title("Elevation (m)")
fig.colorbar(im1, ax=ax1, shrink=0.7)

# Initial soil carbon
im2 = ax2.imshow(inputs["pH"].to_numpy().reshape((9, 9)), extent=extent)
ax2.set_title("Soil pH (-)")
fig.colorbar(im2, ax=ax2, shrink=0.7)

plt.tight_layout()
../_images/f81ff301627e63a336b158462030f2263037cecdf99d6dd1036bcde65954d3b7.png

For some variables, it may be useful to visualise spatial structure in 3 dimensions. The obvious candidate is elevation.

# Extract the elevation data for a 3D plot
top = inputs["elevation"].to_numpy()
x = inputs["x"].to_numpy()
y = inputs["y"].to_numpy()
bottom = np.zeros_like(top)
width = depth = 90

# Make a 3D barplot of the elevation
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(projection="3d")
colors = plt.cm.turbo(top.flatten() / float(top.max()))

poly = ax.bar3d(x, y, bottom, width, depth, top, shade=True, color=colors)
ax.set_title("Elevation (m)")

cell_bounds = range(0, 811, 90)
ax.set_xticks(cell_bounds)
_ = ax.set_yticks(cell_bounds)
../_images/b65b4541ba41f2af6f7af1d77d034e596641fb671efc32f1269b06f830539fc1.png

For other variables, such as air temperature and precipitation, the initial data also provides time series data at reference height that are used to force the simulation across the configured time period.

# Make two side by side plots
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 5))

# Air temperature
ax1.plot(inputs["time_index"], inputs["air_temperature_ref"].T)
ax1.set_title("Air temperature forcing across grid cells")
ax1.set_ylabel("Air temperature (°C)")
ax1.set_xlabel("Time step (months)")

# Precipitation
ax2.plot(inputs["time_index"], inputs["precipitation"].T)
ax2.set_title("Precipitation forcing across grid cells")
ax2.set_ylabel("Total monthly precipitation (mm)")
_ = ax2.set_xlabel("Time step (months)")
../_images/54fd664eba70631c9ea88aba4fb1bd31c7155ee3ebfb30fef739b0da6de8d0ca.png

Model initialisation#

The init group contains the state of model variables that are calculated when the science models initialise. The science models are initialised in a specific order so that each model has access to required initialisation data either in the inputs or by the initialisation of other models.

init = xarray.load_dataset("ve_example/out/model_data.zarr", group="init")
init

<xarray.Dataset> Size: 374kB
Dimensions:                               (cell_id: 81, layers: 14, pft: 2,
                                           element: 3, groundwater_layers: 2,
                                           community_id: 81,
                                           functional_group_id: 20,
                                           time_index: 24)
Coordinates:
  * cell_id                               (cell_id) int64 648B 0 1 2 ... 79 80
    x                                     (cell_id) int32 324B 0 90 ... 630 720
    y                                     (cell_id) int32 324B 720 720 ... 0 0
  * layers                                (layers) int64 112B 0 1 2 ... 11 12 13
    layer_roles                           (layers) <U7 392B 'above' ... 'subs...
  * pft                                   (pft) <U9 72B 'broadleaf' 'shrub'
  * element                               (element) <U1 12B 'C' 'N' 'P'
  * community_id                          (community_id) int64 648B 0 1 ... 80
  * functional_group_id                   (functional_group_id) <U30 2kB 'car...
  * time_index                            (time_index) int64 192B 0 1 ... 22 23
Dimensions without coordinates: groundwater_layers
Data variables: (12/69)
    aerodynamic_resistance_canopy         (cell_id) float64 648B 12.1 ... 12.1
    aerodynamic_resistance_soil           (layers, cell_id) float64 9kB 12.5 ...
    air_temperature                       (layers, cell_id) float64 9kB 22.92...
    arbuscular_mycorrhizal_n_supply       (cell_id) float64 648B 0.002431 ......
    arbuscular_mycorrhizal_p_supply       (cell_id) float64 648B 3.213e-05 .....
    atmospheric_co2                       (layers, cell_id) float64 9kB 400.0...
    ...                                    ...
    total_animal_respiration              (cell_id) float64 648B 0.0 0.0 ... 0.0
    vapour_pressure                       (layers, cell_id) float64 9kB 2.351...
    vapour_pressure_deficit               (layers, cell_id) float64 9kB 0.439...
    vapour_pressure_deficit_ref           (cell_id, time_index) float64 16kB ...
    vapour_pressure_ref                   (cell_id, time_index) float64 16kB ...
    wind_speed                            (layers, cell_id) float64 9kB 0.188...

A good example is the layer heights within the simulation - this is calculated from the configuration of ground layers and the heights of trees in the initial plant community data and provides the vertical structure for the calculation of initial abiotic conditions and other model components

# Loop over vertical layers
for layer in init["layers"]:
    # Get the data for the layer
    layer_data = init["layer_heights"].sel(layers=layer)

    if not np.all(np.isnan(layer_data)):
        # Plot layers that contain some canopy or other data.
        plt.plot(
            init["cell_id"],
            init["layer_heights"].sel(layers=layer),
            label=init["layer_roles"].sel(layers=layer).data,
        )
    plt.xlabel("Cell ID")
    plt.ylabel("Layer height (metres)")
    plt.legend()
../_images/70259e6358309143a05cfafe4d8384204d5822d6e134b664230924dce7ea079e.png

Model outputs#

The main output of interest is the outputs group, which contains variables describing the change in the models through the simulation process.

outputs = xarray.load_dataset("ve_example/out/model_data.zarr", group="outputs")
outputs

<xarray.Dataset> Size: 8MB
Dimensions:                                   (time_index: 24, cell_id: 81,
                                               layers: 14, element: 3, pft: 2,
                                               groundwater_layers: 2)
Coordinates:
  * time_index                                (time_index) int64 192B 0 1 ... 23
  * cell_id                                   (cell_id) int64 648B 0 1 ... 79 80
    x                                         (cell_id) int32 324B 0 90 ... 720
    y                                         (cell_id) int32 324B 720 720 ... 0
  * layers                                    (layers) int64 112B 0 1 ... 12 13
    layer_roles                               (layers) <U7 392B 'above' ... '...
  * element                                   (element) <U1 12B 'C' 'N' 'P'
  * pft                                       (pft) <U9 72B 'broadleaf' 'shrub'
Dimensions without coordinates: groundwater_layers
Data variables: (12/126)
    aerodynamic_resistance_soil               (time_index, cell_id) float64 16kB ...
    air_temperature                           (time_index, layers, cell_id) float64 218kB ...
    animal_arbuscular_mycorrhiza_consumption  (time_index, cell_id) float64 16kB ...
    animal_bacteria_consumption               (time_index, cell_id) float64 16kB ...
    animal_ectomycorrhiza_consumption         (time_index, cell_id) float64 16kB ...
    animal_pom_consumption_cnp                (time_index, cell_id, element) float64 47kB ...
    ...                                        ...
    total_runoff                              (time_index, cell_id) float64 16kB ...
    transpiration                             (time_index, layers, cell_id) float64 218kB ...
    vapour_pressure                           (time_index, layers, cell_id) float64 218kB ...
    vapour_pressure_deficit                   (time_index, layers, cell_id) float64 218kB ...
    vertical_flow                             (time_index, layers, cell_id) float64 218kB ...
    wind_speed                                (time_index, layers, cell_id) float64 218kB ...

This group contains the values of variables as they are updated during the model iteration. All the data in this group has a time dimension, so we can visualise it across time. This can happen in a number of formats such as spatial grids, individual time series within grid cells and as the three dimensional structure of the vertical layers within the simulation. We describe these plotting processes below.

Warning

At present when the Virtual Ecosystem reads in netCDF data it flips it (i.e. the spatial x dimension becomes the y dimension and vice versa). This doesn’t affect the simlation, but does need to be corrected for when analysing the spatial outputs, which will be flipped.

Spatial data#

Using the soil carbon held as mineral-associated organic matter as an example:

# Extract the carbon data
soil_maom_carbon = outputs["soil_cnp_pool_maom"].sel(element="C")

# Make three side by side plots
fig, axes = plt.subplots(ncols=3, figsize=(10, 5))

# Plot start and end MAOM
val_min = soil_maom_carbon.min()
val_max = soil_maom_carbon.max()

# Plot 3 time slices
for idx, ax in zip([0, 10, 23], axes):
    im = ax.imshow(
        soil_maom_carbon[idx, :].to_numpy().reshape((9, 9)),
        extent=extent,
        vmax=val_max,
        vmin=val_min,
    )
    ax.set_title(f"Time step: {idx}")

fig.colorbar(im, ax=axes, orientation="vertical", shrink=0.5)
_ = plt.suptitle("Soil carbon: mineral-associated organic matter", y=0.78, x=0.45)
../_images/ef649aabf32cd0a8162d8c044407998f516e25e4bbf4938f44f5c1941f60623c.png

Temporal data#

The plot below shows the mineral-associated organic matter data as a time series showing the values in each cell across time.

plt.plot(outputs["time_index"], soil_maom_carbon)
plt.xlabel("Time step")
_ = plt.ylabel("Soil carbon as MAOM")
../_images/ddf967b19b9de460b97569cd7d6d7370cec37690465f018bc460639a09ef2fd8.png

Vertical structure#

The layer heights shown in the initial data example above provide the vertical structure used throughout the simulation. The plot below shows the microclimatic estimates of air temperature through the vertical structure of the canopy.

# Extract the x and y location of the grid cell centres and layer heights
# for all observations at a given time step.
time_index = 0

x_3d = (
    outputs["x"]
    .broadcast_like(outputs["layer_heights"][time_index])
    .to_numpy()
    .flatten()
    + 45
)
y_3d = (
    outputs["y"]
    .broadcast_like(outputs["layer_heights"][time_index])
    .to_numpy()
    .flatten()
    + 45
)
z_3d = outputs["layer_heights"][time_index].to_numpy().flatten()

# Extract the air temperature for those points to colour the 3D data.
temp_vals = outputs["air_temperature"][time_index].to_numpy().flatten()

# Generate a 3 dimensional plot of layer heights showing temperature.

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(projection="3d")

cmap = plt.get_cmap("turbo")
paths = ax.scatter(x_3d, y_3d, z_3d, c=temp_vals, cmap=cmap)
fig.colorbar(
    paths,
    ax=ax,
    orientation="vertical",
    shrink=0.6,
    label="Air temperature (°C)",
    pad=0.1,
)

ax.set_xlabel("Easting (m)")
ax.set_ylabel("Northing (m)")
ax.set_zlabel("Layer height (m)")

ax.set_xticks(cell_bounds)
_ = ax.set_yticks(cell_bounds)
../_images/2cb66b5f2a37821b045f3b6cd30a0b7346982bb58b7cfd1768ee810913782e49.png