{ "cells": [ { "cell_type": "markdown", "id": "bd760e84-8f80-473a-9d92-edda06148e4a", "metadata": {}, "source": [ "# VS-Lite in pure Python" ] }, { "cell_type": "markdown", "id": "241a51b1-1b35-4cd1-b2fa-45632c49d63e", "metadata": {}, "source": [ "**Expected time to run through: 5 mins**\n", "\n", "In this tutorial, we demonstrate how to run VS-Lite in pure Python with PyVSL." ] }, { "cell_type": "code", "execution_count": 1, "id": "09954659-a552-4b7a-ab49-dad49544d41a", "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload 2\n", "\n", "import PyVSL\n", "import pandas as pd\n", "import numpy as np" ] }, { "cell_type": "markdown", "id": "a15cf423-d3e0-4f86-9333-cba2a629ddd6", "metadata": {}, "source": [ "In the below, we check the function first, and let's pay attention to the units of the input T and P." ] }, { "cell_type": "code", "execution_count": 2, "id": "7283417f-4b4d-4ea0-8e0f-1921c0da2c92", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[0;31mSignature:\u001b[0m\n", "\u001b[0mPyVSL\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mVSL\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0msyear\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0meyear\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mphi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mT1\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mT2\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m23\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mM1\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.01\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mM2\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.05\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mMmax\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.76\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mMmin\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.01\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0malph\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.093\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mm_th\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m4.886\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mmu_th\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m5.8\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mrootd\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mM0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0msubstep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mI_0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mI_f\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m12\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mhydroclim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'P'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mDocstring:\u001b[0m\n", "Translated from VSLite_v2_3.m - Simulate tree ring width index given monthly climate inputs.\n", "\n", "Basic Usage:\n", " trw = VSLite_v2_3(syear,eyear,phi,T1,T2,M1,M2,T,P)\n", " gives just simulated tree ring as ouput.\n", "\n", " [trw,gT,gM,gE,M] = VSLite_v2_3(syear,eyear,phi,T1,T2,M1,M2,T,P))\n", " also includes growth response to temperature, growth response to soil\n", " moisture, scaled insolation index, and soil moisture estimate in outputs.\n", "\n", "Basic Inputs:\n", " syear = start year of simulation.\n", " eyear = end year of simulation.\n", " phi = latitude of site (in degrees N)\n", " T1 = scalar temperature threshold below which temp. growth response is zero (in deg. C)\n", " T2 = scalar temperature threshold above which temp. growth response is one (in deg. C)\n", " M1 = scalar soil moisture threshold below which moist. growth response is zero (in v/v)\n", " M2 = scalar soil moisture threshold above which moist. growth response is one (in v/v)\n", " (Note that optimal growth response parameters T1, T2, M1, M2 may be estimated\n", " using code estimate_vslite_params_v2_3.m also freely available at\n", " the NOAA NCDC Paleoclimatology software library.)\n", " T = (12 x Nyrs) matrix of ordered mean monthly temperatures (in degEes C)\n", " P = (12 x Nyrs) matrix of ordered accumulated monthly precipitation (in mm)\n", "\n", "Advanced Inputs (must be specified as property/value pairs):\n", " 'lbparams': Parameters of the Leaky Bucket model of soil moisture.\n", " These may be specified in an 8 x 1 vector in the following\n", " order (otherwise the default values are read in):\n", " Mmax: scalar maximum soil moisture content (in v/v),\n", " default value is 0.76\n", " Mmin: scalar minimum soil moisture (in v/v), default\n", " value is 0.01\n", " alph: scalar runoff parameter 1 (in inverse months),\n", " default value is 0.093\n", " m_th: scalar runoff parameter 3 (unitless), default\n", " value is 4.886\n", " mu_th: scalar runoff parameter 2 (unitless), default\n", " value is 5.80\n", " rootd: scalar root/\"bucket\" depth (in mm), default\n", " value is 1000\n", " M0: initial value for previous month's soil moisture at\n", " t = 1 (in v/v), default value is 0.2\n", " substep: logical 1 or 0; perform monthly substepping in\n", " leaky bucket (1) or not (0)? Default value is 0.\n", " 'intwindow': Integration window. Which months' growth responses should\n", " be intregrated to compute the annual ring-width index?\n", " Specified as a 2 x 1 vector of integer values. Both\n", " elements are given in integer number of months since January\n", " (July) 1st of the current year in the Northern (Southern)\n", " hemisphere, and specify the beginning and end of the integration\n", " window, respectively. Defaults is [1 ; 12] (eg. integrate\n", " response to climate over the corresponding calendar year,\n", " assuming location is in the northern hemisphere).\n", " 'hydroclim': Value is a single character either taking value ['P'] or ['M'].\n", " If ['M'], then 9th input is interpreted as an estimate of\n", " soil moisture content (in v/v) rather than as precipitation.\n", " Model default is to read in precipitation and use the CPC's\n", " Leaky Bucket model of hydrology to estimate soil moisture,\n", " however if soil moisture observations or a more sophisticated\n", " estimate of moisture accounting for snow-related processes\n", " is available, then using these data directly are recommended\n", " (and will also speed up code).\n", "\n", " For more detailed documentation, see:\n", " 1) Tolwinski-Ward et al., An efficient forward model of the climate\n", " controls on interannual variation in tree-ring width, Climate Dynamics (2011)\n", " DOI: 10.1007/s00382-010-0945-5\n", "\n", " 2) Tolwinski-Ward et al., Erratum to: An efficient forward model of the climate\n", " controls on interannual variation in tree-ring width, Climate Dynamics (2011)\n", " DOI: 10.1007/s00382-011-1062-9\n", "\n", " 3) Tolwinski-Ward et al., Bayesian parameter estimation and\n", " interpretation for an intermediate model of tree-ring width, Clim. Past\n", " (2013), DOI: 10.5194/cp-9-1-2013\n", "\n", " 4) Documentation available with the model at http://www.ncdc.noaa.gov/paleo/softlib/softlib.html\n", "\u001b[0;31mFile:\u001b[0m ~/Github/PyVSL/PyVSL/core.py\n", "\u001b[0;31mType:\u001b[0m function\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "PyVSL.VSL?" ] }, { "cell_type": "markdown", "id": "75577926-587f-477b-ad60-2e8643d8d0af", "metadata": {}, "source": [ "## Load test data\n", "\n", "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)." ] }, { "cell_type": "code", "execution_count": 3, "id": "c778b014-f1ff-455d-a5d4-dd9076dbbebb", "metadata": {}, "outputs": [], "source": [ "data_dict = pd.read_pickle('./data/test_T_P.pkl')\n", "time = data_dict['time']\n", "T = data_dict['T'] # unit: degC\n", "P = data_dict['P'] # unit: mm/month" ] }, { "cell_type": "markdown", "id": "446dd406-bcd7-4f4f-aa07-a5f62b7c62d4", "metadata": {}, "source": [ "## Run VS-Lite in pure Python\n", "\n", "Now we are ready to run the wrapper to generate simulated TRW." ] }, { "cell_type": "code", "execution_count": 4, "id": "3bcfbe86-36c9-4e1f-8e54-aa114ab37698", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "dict_keys(['trw', 'gT', 'gM', 'gE', 'Gr', 'M', 'potEv', 'width', 'width_mean', 'width_std'])\n", "[ 0.03243917 1.62155893 2.00388961 ... -0.873099 0.07933166\n", " -1.2066474 ]\n", "CPU times: user 306 ms, sys: 2.23 ms, total: 308 ms\n", "Wall time: 308 ms\n" ] } ], "source": [ "%%time\n", "\n", "T1, T2 = 1, 15\n", "M1, M2 = 0.01, 0.05\n", "\n", "res = PyVSL.VSL(\n", " 850, 1850, 45, # the starting year and ending year of the input T & P, along with the latitude\n", " T, P,\n", " T1=T1, T2=T2, M1=M1, M2=M2, # parameters of the thresholds\n", ")\n", "print(res.keys())\n", "TRW = res['trw'] # this is the simulated TRW\n", "print(TRW)" ] }, { "cell_type": "code", "execution_count": null, "id": "14ab812f-fab7-407b-9f07-29cc00a0b8f2", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.7" } }, "nbformat": 4, "nbformat_minor": 5 }