This tutorial demonstrates implementation of the Cellular Automaton Tree-GRass-Shrub Simulator (CATGRaSS) [Zhou et al., 2013] on a digital elevation model (DEM). This model is built using components from the Landlab component library. CATGRaSS is a spatially explicit model of plant coexistence. It simulates local ecohydrologic dynamics (soil moisture, transpiration, biomass) and spatial evolution of tree, grass, and shrub Plant Functional Types (PFT) driven by rainfall and solar radiation.
Each cell in the model grid can hold a single PFT or remain empty. Tree and shrub plants disperse seeds to their neighbors. Grass seeds are assumed to be available at each cell. Establishment of plants in empty cells is determined probabilistically based on water stress for each PFT. Plants with lower water stress have higher probability of establishment. Plant mortality is simulated probabilistically as a result of aging and drought stress. Fires and grazing will be added to this model soon.
This model (driver) contains:
A local vegetation dynamics model that simulates storm and inter-storm water balance and ecohydrologic fluxes (ET, runoff), and plant biomass dynamics by coupling the following components:
A spatially explicit probabilistic cellular automaton component that simulates plant competition by tracking establishment and mortality of plants based on soil moisture stress:
To run this Jupyter notebook, please make sure that the following files are in the same folder:
[Ref: Zhou, X, E. Istanbulluoglu, and E.R. Vivoni. "Modeling the ecohydrological role of aspect-controlled radiation on tree-grass-shrub coexistence in a semiarid climate." Water Resources Research 49.5 (2013): 2872-2895]
In this tutorial, we are going to work with a landscape in central New Mexico, USA, where aspect controls the organization of PFTs. The climate in this area is semi-arid with Mean Annual Precipitation (MAP) of 254 mm [Zhou et. al 2013]. We will do the following:
Import the required libraries:
In [ ]:
from __future__ import print_function
%matplotlib inline
In [ ]:
import time
import numpy as np
from landlab.io import read_esri_ascii
from landlab import RasterModelGrid as rmg
from landlab import load_params
from Ecohyd_functions_DEM import (Initialize_, Empty_arrays, Create_PET_lookup,
Save_, Plot_)
Note: 'Ecohyd_functions_DEM.py' is a utility script that contains 'functions', which instantiates components and manages inputs and outputs, and help keep this driver concise. Contents of 'Ecohyd_functions_DEM.py' can be a part of this driver (current file), however left out to keep driver concise.
We will use two grids in this driver. One grid will represent the actual landscape or domain (e.g., created from a DEM). Another grid, with one cell for each of the plant functional types (PFTs), will be used to create Potential Evapotranspiration (PET) lookup tables.
grid: This grid represents the actual landscape. Each cell can be occupied by a single PFT such as tree, shrub, grass, or can be empty (bare). In this example we assume that the elevation field and the vegetation field has the same resolution.
grid1: This grid will be used to compute plant-specific PET at a point. Spatially distributed PET Lookup arrays (for all days of the year) will be created for each PFT based on these point values.
Note: In this tutorial, the physical ecohydrological components and cellular automata plant competition will be run on grids with same resolution. To develop differential spatial resolutions for the two models, see the tutorial 'cellular_automaton_vegetation_flat.ipynb'.
In [ ]:
(grid, elevation) = read_esri_ascii('DEM_10m.asc') # Read the DEM
grid1 = rmg((5, 4), xy_spacing=(5., 5.)) # Representative grid
Include the input file that contains all input parameters needed for all components. This file can either be a Python dictionary or a text file that can be converted into a Python dictionary. If a text file is provided, it will be converted to a Python dictionary. Here we use an existing text file prepared for this exercise.
In [ ]:
InputFile = 'Inputs_Vegetation_CA_DEM.txt'
data = load_params(InputFile) # Creates dictionary that holds the inputs
Instantiate Landlab components to simulate corresponding attributes. In this example, we shall demonstrate the use of seasonal rainfall and PFT-specific potential evapotranspiration. The instantiated objects are:
Note: Almost every component in Landlab is coded as a 'class' (to harness the advantages of object oriented programming). An 'object' is the instantiation of the 'class' (for more information, please refer any object oriented programming book). A 'field' refers to a Landlab field (please refer to the Landlab documentation to learn more about Landlab fields).
Now let's instantiate all Landlab components that we are going to use for this tutorial:
In [ ]:
PD_D, PD_W, Rad, Rad_PET, PET_Tree, PET_Shrub, PET_Grass, SM, VEG, vegca = (
Initialize_(data, grid, grid1, elevation))
Lets look at the initial organization of PFTs
In [ ]:
import matplotlib.pyplot as plt
from landlab.plot import imshow_grid
import matplotlib as mpl
cmap = mpl.colors.ListedColormap(
['green', 'red', 'black', 'white', 'red', 'black'])
bounds = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5]
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
description = 'green: grass; red: shrub; black: tree; white: bare'
plt.figure(101)
imshow_grid(grid,
'vegetation__plant_functional_type',
values_at='cell',
cmap=cmap,
grid_units=('m', 'm'),
norm=norm,
limits=[0, 5],
allow_colorbar=False)
plt.figtext(0.2, 0.0, description, weight='bold', fontsize=10)
Specify an approximate number of years for the model to run.
IMPORTANT: This code in numerically extensive. It might take an hour or more to run this simulation for 300 years. It is suggested to run the simulation for 50 years which might take less than 7 minutes to execute.
In [ ]:
n_years = 50 # Approx number of years for model to run
# Calculate approximate number of storms per year
fraction_wet = (data['doy__end_of_monsoon'] -
data['doy__start_of_monsoon']) / 365.
fraction_dry = 1 - fraction_wet
no_of_storms_wet = (8760 * (fraction_wet) /
(data['mean_interstorm_wet'] + data['mean_storm_wet']))
no_of_storms_dry = (8760 * (fraction_dry) /
(data['mean_interstorm_dry'] + data['mean_storm_dry']))
n = int(n_years * (no_of_storms_wet + no_of_storms_dry))
Create empty arrays to store spatio-temporal data over multiple iterations. The captured data can be used for plotting model outputs.
In [ ]:
P, Tb, Tr, Time, VegType, PET_, Rad_Factor, EP30, PET_threshold = (
Empty_arrays(n, n_years, grid, grid1))
To reduce computational overhead, we shall create a lookup array for plant-specific PET values for each day of the year, and slope and aspect grid.
In [ ]:
Create_PET_lookup(Rad, PET_Tree, PET_Shrub, PET_Grass, PET_, Rad_Factor, EP30,
Rad_PET, grid)
Specify current_time (in years). current_time is the current time in the simulation.
In [ ]:
# # Represent current time in years
current_time = 0 # Start from first day of Jan
# Keep track of run time for simulation—optional
Start_time = time.clock() # Recording time taken for simulation
# declaring few variables that will be used in storm loop
time_check = 0. # Buffer to store current_time at previous storm
yrs = 0 # Keep track of number of years passed
WS = 0. # Buffer for Water Stress
Tg = 365 # Growing season in days
The loop below couples the components introduced above in a for loop until all "n" number of storms are generated. Time is advanced by the soil moisture object based on storm and interstorm durations that are estimated by the strom generator object. The ecohydrologic model is run on each storm whereas cellular automaton vegetation component is run once every year.
Note: This loop might take around 10 minutes (depending on your computer) to run for a 50 year simulation. Ignore any warnings you might see.
In [ ]:
# # Run storm Loop
for i in range(0, n):
# # Update objects
# Calculate Day of Year (DOY)
Julian = np.int(np.floor((current_time - np.floor(current_time)) * 365.))
# Generate seasonal storms
# for Dry season
if Julian < data['doy__start_of_monsoon'] or Julian > data[
'doy__end_of_monsoon']:
PD_D.update()
P[i] = PD_D.get_storm_depth()
Tr[i] = PD_D.get_precipitation_event_duration()
Tb[i] = PD_D.get_interstorm_event_duration()
# Wet Season—Jul to Sep—NA Monsoon
else:
PD_W.update()
P[i] = PD_W.get_storm_depth()
Tr[i] = PD_W.get_precipitation_event_duration()
Tb[i] = PD_W.get_interstorm_event_duration()
# Spatially distribute PET and its 30-day-mean (analogous to degree day)
grid['cell']['surface__potential_evapotranspiration_rate'] = (
(np.choose(grid['cell']['vegetation__plant_functional_type'],
PET_[Julian])) * Rad_Factor[Julian])
grid['cell']['surface__potential_evapotranspiration_30day_mean'] = (
(np.choose(grid['cell']['vegetation__plant_functional_type'],
EP30[Julian])) * Rad_Factor[Julian])
# Assign spatial rainfall data
grid['cell']['rainfall__daily_depth'] = P[i] * np.ones(
grid.number_of_cells)
# Update soil moisture component
current_time = SM.update(current_time, Tr=Tr[i], Tb=Tb[i])
# Decide whether its growing season or not
if Julian != 364:
if EP30[Julian + 1, 0] > EP30[Julian, 0]:
PET_threshold = 1
# 1 corresponds to ETThresholdup (begin growing season)
else:
PET_threshold = 0
# 0 corresponds to ETThresholddown (end growing season)
# Update vegetation component
VEG.update(PETthreshold_switch=PET_threshold, Tb=Tb[i], Tr=Tr[i])
# Update yearly cumulative water stress data
WS += (grid['cell']['vegetation__water_stress']) * Tb[i] / 24.
# Record time (optional)
Time[i] = current_time
# Cellular Automata
if (current_time - time_check) >= 1.:
if yrs % 5 == 0:
print('Elapsed time = {time} years'.format(time=yrs))
VegType[yrs] = grid['cell']['vegetation__plant_functional_type']
grid['cell']['vegetation__cumulative_water_stress'] = WS / Tg
vegca.update()
SM.initialize()
VEG.initialize()
time_check = current_time
WS = 0
yrs += 1
VegType[yrs] = grid['cell']['vegetation__plant_functional_type']
Time_Consumed is an optional variable that gives information about computer running time
In [ ]:
Final_time = time.clock()
Time_Consumed = (Final_time - Start_time) / 60. # in minutes
print('Time_consumed = {time} minutes'.format(time=Time_Consumed))
Save the outputs using numpy.save(). These files have '.nc' extension, which can be loaded using numpy.load().
In [ ]:
# # Saving
sim = 'VegCA_DEM_26Jul16_'
# Save_(sim, Tb, Tr, P, VegType, yrs, Time_Consumed, Time)
Lets look at outputs.
Plots of the cellular field of PFT at specified year step can be found below where:
GRASS = Green; SHRUB = Red; TREE = Black; BARE = White;
At the end, percentage cover for each PFT is plotted with respect to time.
In [ ]:
Plot_(grid, VegType, yrs, yr_step=10)
In [ ]:
from IPython.display import Image
Image(filename='presentation.png')
If you want to explore this model further, open 'Inputs_Vegetation_CA.txt' and change the input parameters (e.g., initial PFT distribution percentages, storm characteristics, etc..).