Quickstart: the wrapper for VS-Lite in Matlab

Expected time to run through: 5 mins

In this tutorial, we demonstrate how to run the R package of VS-Lite with PyVSL.

%load_ext autoreload
%autoreload 2

import PyVSL
import pandas as pd
import numpy as np

In the below, we check the wrapper function first, and let’s pay attention to the units of the input T and P.

PyVSL.VSL_M?

Signature:
PyVSL.VSL_M(
    syear,
    eyear,
    phi,
    T,
    P,
    T1=8,
    T2=23,
    M1=0.01,
    M2=0.05,
    Mmax=0.76,
    Mmin=0.01,
    alph=0.093,
    m_th=4.886,
    mu_th=5.8,
    rootd=1000,
    M0=0.2,
    substep=0,
    I_0=1,
    I_f=12,
    hydroclim='P',
)
Docstring:
VS-Lite tree-ring PSM

Porting the VSLite Matlab code by Suz Tolwinski-Ward (https://github.com/suztolwinskiward/VSLite/blob/master/VSLite_v2_3.m)

NOTE: Need to install Octave as a substitute of Matlab

Args:
    syear (int): Start year of simulation.
    eyear (int): End year of simulation.
    phi (float): Latitude of site (in degrees N).
    T (array): temperature timeseries in deg C, with length of 12*Nyrs
    P (array): precipitation timeseries in mm/month, with length of 12*Nyrs
    # original param. in R package: T (12 x Nyrs) Matrix of ordered mean monthly temperatures (in degEes C).
    # original param. in R package: P (12 x Nyrs) Matrix of ordered accumulated monthly precipitation (in mm).
    T1: Lower temperature threshold for growth to begin (scalar, deg. C).
    T2: Upper temperature threshold for growth sensitivity to temp (scalar, deg. C).
    M1: Lower moisture threshold for growth to begin (scalar, v.v).
    M2: Upper moisture threshold for growth sensitivity to moisture (scalar, v/v).
    Mmax: Scalar maximum soil moisture held by the soil (in v/v).
    Mmin: Scalar minimum soil moisture (for error-catching) (in v/v).
    alph: Scalar runoff parameter 1 (in inverse months).
    m_th: Scalar runoff parameter 3 (unitless).
    mu_th: Scalar runoff parameter 2 (unitless).
    rootd: Scalar root/"bucket" depth (in mm).
    M0: Initial value for previous month's soil moisture at t = 1 (in v/v).
    substep: Use leaky bucket code with sub-monthly time-stepping? (TRUE/FALSE)
    I_0: lower bound of integration window (months before January in NH)
    I_f: upper bound of integration window (months after January in NH)
    hydroclim: Switch; value is either "P" (default) or "M" depending on whether the second input climate variable
        is precipitation, in which case soil moisture isestimated using the Leaky Bucket model of the CPC,
        or soil moisture, in which case the inputs are used directly to compute the growth response.

Returns:
    res (dict): a dictionary with several lists, keys including:
        trw: tree ring width
        gT, gM, gE, M, potEv, sample.mean.width, sample.std.width

Reference:
    + The original VSLite Matlab code by Suz Tolwinski-Ward (https://github.com/suztolwinskiward/VSLite/blob/master/VSLite_v2_3.m)
    + Tolwinski-Ward, S.E., M.N. Evans, M.K. Hughes, K.J. Anchukaitis, An efficient forward model of the climate controls
        on interannual variation in tree-ring width, Climate Dynamics, doi:10.1007/s00382-010-0945-5, (2011).
    + Tolwinski-Ward, S.E., K.J. Anchukaitis, M.N. Evans, Bayesian parameter estimation and interpretation for an
        intermediate model of tree-ring width , Climate of the Past, doi:10.5194/cpd-9-615-2013, (2013).
    + Tolwinski-Ward, S.E., M.P. Tingley, M.N. Evans, M.K. Hughes, D.W. Nychka, Probabilistic reconstructions of localtemperature and
        soil moisture from tree-ring data with potentially time-varying climatic response, Climate Dynamics, doi:10.1007/s00382-014-2139-z, (2014).
File:      ~/Github/PyVSL/PyVSL/wrapper.py
Type:      function

Load test data

First of all, we load the test data (can be downloaded with this link), which contains T and P at a grid point of the CCSM4 past1000 simulation. Note that the unit for T is degree C and that for P is accumulated monthly precipitation (mm/month).

data_dict = pd.read_pickle('./data/test_T_P.pkl')
time = data_dict['time']
T = data_dict['T'] # unit: degC
P = data_dict['P'] # unit: mm/month

print(time)
print(np.min(T), np.max(T))
print(np.min(P), np.max(P))

[ 850.04109589  850.12328767  850.20273973 ... 1850.7890411  1850.8739726
 1850.95616438]
-31.395203 16.369354
1.3947426 362.0332

Run the wrapper

Now we are ready to run the wrapper to generate simulated TRW.

%%time

res = PyVSL.VSL_M(
        850, 1850, 45,                  # the starting year and ending year of the input T & P, along with the latitude
        T, P,
        T1=1, T2=15, M1=0.01, M2=0.05,  # parameters of the thresholds
)

CPU times: user 7.14 ms, sys: 11.4 ms, total: 18.6 ms
Wall time: 1.64 s

print(res.keys())

dict_keys(['trw', 'gT', 'gM', 'gE', 'M', 'potEv', 'width', 'width_mean', 'width_std', 'Gr'])

print(res['trw'])  # this is the simulated TRW

[ 0.03242296  1.62074876  2.00288841 ... -0.87266278  0.07929202
 -1.20604453]