<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <a href="http
Step1:  Note
Step2:  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.
Step3:  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
Step4:  Lets look at the initial organization of PFTs
Step5:  Specify an approximate number of years for the model to run.
IMPORTANT
Step6:  Create empty arrays to store spatio-temporal data over multiple iterations. The captured data can be used for plotting model outputs.
Step7:  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.
Step8:  Specify current_time (in years). current_time is the current time in the simulation.
Step9:  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
Step10:  Time_Consumed is an optional variable that gives information about computer running time
Step11:  Save the outputs using numpy.save(). These files have '.nc' extension, which can be loaded using numpy.load().
Step12:  Lets look at outputs.
Plots of the cellular field of PFT at specified year step can be found below where
Step13:  If you run this model for around 900 years, you will observe patterns of PFTs. For example, you will find more trees on north facing slopes and mostly shrubs and grass on south facing slopes, as shown below
<ASSISTANT_TASK:>
Python Code:
from __future__ import print_function
%matplotlib inline
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_,
)
Explanation: <a href="http://landlab.github.io"><img style="float: left" src="../../landlab_header.png"></a>
WARNING: This tutorial has not been updated to work with Landlab 2.0 and is thus not tested to verify that it will run. 
Tutorial For Cellular Automaton Vegetation Model Coupled With Ecohydrologic Model
<hr>
<small>For more Landlab tutorials, click here: <a href="https://landlab.readthedocs.io/en/v2_dev/user_guide/tutorials.html">https://landlab.readthedocs.io/en/v2_dev/user_guide/tutorials.html</a></small>
<hr>
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:
    - PrecipitationDistribution
    - Radiation
    - PotentialEvapotranspiration
    - SoilMoisture
    - Vegetation
A spatially explicit probabilistic cellular automaton component that simulates plant competition by tracking establishment and mortality of plants based on soil moisture stress:
VegCA
To run this Jupyter notebook, please make sure that the following files are in the same folder:
 - cellular_automaton_vegetation_DEM.ipynb (this notebook)
 - Inputs_Vegetation_CA.txt (Input parameters for the model)
 - Ecohyd_functions_DEM.py (Utility functions)
[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 a landscape 
- Initialize the landscape with random distribution of PFTs
- Run the coupled Ecohydrology and cellular automata plant competition model for 50 years
- Visualize and examine outputs
Let's walk through the code:
Import the required libraries:
End of explanation
(grid, elevation) = read_esri_ascii("DEM_10m.asc")  # Read the DEM
grid1 = rmg((5, 4), xy_spacing=(5.0, 5.0))  # Representative grid
Explanation: 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'.
End of explanation
InputFile = "Inputs_Vegetation_CA_DEM.txt"
data = load_params(InputFile)  # Creates dictionary that holds the inputs
Explanation: 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.
End of explanation
PD_D, PD_W, Rad, Rad_PET, PET_Tree, PET_Shrub, PET_Grass, SM, VEG, vegca = Initialize_(
    data, grid, grid1, elevation
)
Explanation: 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:
- PD_D: object for dry season rainfall, 
- PD_W: object for wet season rainfall, 
- Rad: Radiation object computes radiation factor defined as the ratio of total shortwave radiation incident on a sloped surface to total shortwave radiation incident on a flat surface. 
- Rad_PET: Representative radiation object which is used only as an input for PET.
- PET_PFT: Plant specific PET objects (we use a cosine function, fitted to calculated PET, as a function of Day Of the Year (DOY) to reduce computation overhead). This value is spatially distributed by using a radiation factor.
- SM: Soil Moisture object simulates root-zone average soil moisture at each cell using inputs of potential evapotranspiration, live leaf area index, and vegetation cover.
- VEG: Vegetation dynamics object simulates net primary productivity and biomass and thus leaf area index at each cell based on inputs of root-zone average soil moisture.
- vegca: Cellular Automaton plant competition object. This object simulates the spatial dynamics of PFTs. It is run once every year at the end of the growing season. This object is initialized with a random cellular field of PFT. Each year, this object updates the field of PFTs based on probabilistic establishment and mortality rules employed at each cell of the modeled DEM.
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:
End of explanation
import matplotlib.pyplot as plt
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)
grid.imshow(
    "vegetation__plant_functional_type",
    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)
Explanation: Lets look at the initial organization of PFTs
End of explanation
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.0
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))
Explanation: 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.
End of explanation
P, Tb, Tr, Time, VegType, PET_, Rad_Factor, EP30, PET_threshold = Empty_arrays(
    n, n_years, grid, grid1
)
Explanation: Create empty arrays to store spatio-temporal data over multiple iterations. The captured data can be used for plotting model outputs.
End of explanation
Create_PET_lookup(
    Rad, PET_Tree, PET_Shrub, PET_Grass, PET_, Rad_Factor, EP30, Rad_PET, grid
)
Explanation: 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.
End of explanation
# # 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.0  # Buffer to store current_time at previous storm
yrs = 0  # Keep track of number of years passed
WS = 0.0  # Buffer for Water Stress
Tg = 365  # Growing season in days
Explanation: Specify current_time (in years). current_time is the current time in the simulation.
End of explanation
# # Run storm Loop
for i in range(0, n):
    # # Update objects
    # Calculate Day of Year (DOY)
    Julian = int(np.floor((current_time - np.floor(current_time)) * 365.0))
    # 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.0
    # Record time (optional)
    Time[i] = current_time
    # Cellular Automata
    if (current_time - time_check) >= 1.0:
        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"]
Explanation: 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.
End of explanation
Final_time = time.clock()
Time_Consumed = (Final_time - Start_time) / 60.0  # in minutes
print("Time_consumed = {time} minutes".format(time=Time_Consumed))
Explanation: Time_Consumed is an optional variable that gives information about computer running time
End of explanation
# # Saving
sim = "VegCA_DEM_26Jul16_"
# Save_(sim, Tb, Tr, P, VegType, yrs, Time_Consumed, Time)
Explanation: Save the outputs using numpy.save(). These files have '.nc' extension, which can be loaded using numpy.load().
End of explanation
Plot_(grid, VegType, yrs, yr_step=10)
Explanation: 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.
End of explanation
from IPython.display import Image
Image(filename="presentation.png")
Explanation: If you run this model for around 900 years, you will observe patterns of PFTs. For example, you will find more trees on north facing slopes and mostly shrubs and grass on south facing slopes, as shown below:
End of explanation
<END_TASK>