{ "cells": [ { "cell_type": "markdown", "id": "ac04081b", "metadata": {}, "source": [ "\n", " \n", " \n", " \n", " \n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", "
" ] }, { "cell_type": "markdown", "id": "88f75198", "metadata": {}, "source": [ "## Learning objective:\n", "* Introducing stoPET model\n", "* learn the input parameters used in the model\n", "* changing parameter values\n", "* generating single point PET value \n", "* learn data visulaization from stoPET output\n", "\n", "est. time 2 hour" ] }, { "cell_type": "markdown", "id": "68bc896b", "metadata": {}, "source": [ "CUWALID is an integrated model used to obtain actionable forecasts for hydrological components in drylands. It consists of a main hydrological model (DRYP) which allows to estimate the partioning of the water balance. This hydrological model is driven by two major climate inputs **precipitation** and **potential evapotranspiration (PET)**. The driving climate variables are obtained based on stochastic models integrated in the system. **STORM** is the precipitation model and **stoPET** is the PET model that generates the required driving climate variables. Here, we disscuss the stoPET model and how it works in the CUWALID system." ] }, { "cell_type": "markdown", "id": "37dbb309", "metadata": {}, "source": [ "## 1. stoPET\n", "stoPET is a stochastic PET generator over the globe (55N - 55S). The model consists of two scripts **run_stoPET.py** and **stoPET_v1.py** where the first one is used to provide the nucessary inputs and run the model. The second one is the script containing the functions to generate the PET values. For the trainning purpose the function in the run_stoPET.py is provided below. Details of the model description can be found in the following paper stoPET paper and you can download the model from this link.\n", "\n", "**how hourly PET (hPET) was generated?**\n", "\n", "PET was estimated using the ERA5-Land dataset as a driving variable. The PET dataset generated using the Penmann-Monteith method is provided as an hourly and daily data (hPET and dPET). The paper explaining details of the method and the data is given in\n", "hPET paper. \n", "\n", "Estimating PET using Penmann-Monteith method requires 7 input variables:\n", "\n", "* 10 m u-component (zonal) of wind speed [m s−1]\n", "\n", "* 10 m v-component (meridional) of wind speed [m s−1]\n", "\n", "* 2 m dew point temperature [K]\n", "\n", "* 2 m air temperature [K]\n", "\n", "* surface net solar radiation [J m−2]\n", "\n", "* surface net thermal radiation [J m−2]\n", "\n", "* atmospheric surface pressure [Pa]\n", "\n", "The equation used is given as: \n", "\n", "## hPET = [0.408 * ∆ (Rn - G) + γ(37/Ta + 273)* u2(es - ea)] / [∆ + γ(1 + 0.34u2)]\n", "\n", "where Rn is hourly net radiation (MJ m−2), G the soil heat flux (MJ m−2), γ is the psychrometric constant (kPa °C−1), ∆ is slope of saturation vapour pressure curve (kPa °C−1), Ta is hourly air temperature (°C) after converting from ERA5-Land temperature in K, es is hourly saturation vapour pressure (kPa), ea is actual hourly vapour pressure (kPa), and u2 is the hourly wind speed (m s−1) at 2 m above the land surface. \n", "\n", "**how stoPET was generated?**\n", "The stoPET model generates hourly PET values based on sine function parameters estimated from hPET. The resulting PET generated from stoPET retains the diurnal and seasonal variations in PET contained within the hPET dataset, but notably, stoPET injects randomness (stochasticity) in the simulated series via a noise factor.\n", "\n", "## Y = A sin (B × t + C) + D\n", "\n", "where A represents the diurnal amplitude (mm h−1), B is thefrequency (h−1), C is the phase shift (–), and D is the vertical\n", "shift (mm h−1). t is time (h), and Y is the new PET value(mm h−1) generated from the sine function.\n", "\n", "The overall equation is given as follows:\n", "\n", "## Stochastic PET = (average diurnal cycle of PET using a sine function × a random noise ratio) + user-defined annual PET variability.\n", "\n", "It has three components:\n", "* 1. the sine function \n", "* 2. the random noise ratio \n", "* 3. the adjustment for temperature." ] }, { "cell_type": "code", "execution_count": null, "id": "276ff233", "metadata": {}, "outputs": [], "source": [ "# This is for showing plots \n", "# interactively in Jupiter notebook\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "id": "8764b0fe", "metadata": {}, "outputs": [], "source": [ "import sys\n", "import numpy as np\n", "from stoPET_v1 import *\n", "from mpl_toolkits.basemap import Basemap" ] }, { "cell_type": "code", "execution_count": null, "id": "ceb02d7e", "metadata": {}, "outputs": [], "source": [ "def run_stoPET():\n", " ## ----- CHANGE THE INPUT VARIABLES HERE -----##\n", " datapath = 'stopet_parameter_files/'\n", " outputpath = 'results/'\n", " runtype = 'single' #'regional' #\n", " startyear = 2000\n", " endyear = 2010\n", "\n", " # Single point stoPET run\n", " latval = 1.73\n", " lonval = 40.09\n", "\n", " # Regional stoPET run\n", " latval_min = -5.5\n", " latval_max = -4.5 #5.5\n", " lonval_min = 33.0\n", " lonval_max = 34.5 #42.0\n", " locname = 'Wajir' #'Turkana1' #\n", "\n", " number_ensm = 2\n", " tempAdj = 3\n", " deltat = 1.5\n", " udpi_pet = 5\n", "\n", " ## ------ NO CHANGES BELLOW THIS -------------##\n", " if runtype == 'single':\n", " for ens_num in np.arange(0,number_ensm):\n", " stoPET_wrapper_singlepoint(startyear, endyear, \n", " latval, lonval, locname,\n", " ens_num,datapath, outputpath, \n", " tempAdj, deltat, udpi_pet)\n", " elif runtype == 'regional':\n", " for ens_num in np.arange(0,number_ensm):\n", " stoPET_wrapper_regional(startyear, endyear, \n", " latval_min, latval_max, \n", " lonval_min, lonval_max,\n", " locname, ens_num, datapath, \n", " outputpath, \n", " tempAdj, deltat, udpi_pet)\n", " else:\n", " raise ValueError('runtype only takes single and regional ... please check!')\n", "\n" ] }, { "cell_type": "markdown", "id": "40e3feda", "metadata": {}, "source": [ "### 1.1 Changing input parameters\n", "The above function is used to run the stoPET model and generate the required PET values. There are input parameters that are required to run it. The user can change these parameters in this function. stoPET can run in two types: one is a **single point run** and two is a **regional run**. If a single point run is selected, the user will need to provide the **lat/lon** of the specific location and the model will choose the nearest grid point to generate the PET value. If a regional run is selected, the user will need to provide the **four corners of a rectangle** that contain the region of interest. The user also needs to provide the **start year** and **end year** for the data and a **local name** to differentiate the region or location.\n", "\n", "### 1.2 Adjusting for temperature increase\n", "If the user wants to adjust the PET values to account for increasing temperature, they can provide one of the three mothods included in the model **(tempAdj)**. This is an integer number with values 1, 2, or 3. Each of these numbers represents what method to use for the model to account for temperature adjustment on future PET. Method 1 = 1, Method 2 = 2, Method 3 = 3 (Refer to the stoPET paper for the description of each method).\n", "\n", " 1. Method 1: user-defined percentage step change in annual PET\n", " 2. Method 2: step change in PET based on a user-defined change in atmospheric temperature\n", " 3. Method 3: progressive change in PET based on the historical trend in hPET\n", "\n", "If the user chooses Method 1, the value used to increase the PET given as a percentage **(udpi_pet)**, will be used. If Method 2 is chosen, the value of the increased temperature given as **(deltat)** will be used. If Method 3 is chosen, **(deltat)** will be used and the **(udpi_pet)** will be ignored. \n", "\n", "The user also needs to provide how many ensembles of run are needed **(number_ensm)**. Once the user provides the required input values, the model can be executed using the following function." ] }, { "cell_type": "code", "execution_count": null, "id": "707ff88f", "metadata": {}, "outputs": [], "source": [ "run_stoPET()" ] }, { "cell_type": "markdown", "id": "f0422346", "metadata": {}, "source": [ "As seen in the output value, the model generated eleven years (2000 - 2010) of single point PET. stoPET runs twice as the number of ensembles given is two." ] }, { "cell_type": "markdown", "id": "68b85bfa", "metadata": {}, "source": [ "### 1.3 Model Output\n", "Once we run the model the outputs will be saved in the output folder user provided. The output of these functions will be written in the output directory provided, where the model creates a new named directory `outputpath+locname+_E + number_ensm +_stoPET/`. \n", " \n", "For the single point run, stoPET generates text files `year_latval+_+lonval+_+tempAdj+stoPET.txt` and `year_latval+_+lonval+_+tempAdj+AdjstoPET.txt`. The first file is the PET generated without accounting for any adjustment for temperature. The second file is the adjusted PET based on the user’s choice (tempAdj). If one wants to avoid any temperature change adjustment in the model, just use the first file output and ignore the second file.\n", "\n", "**Example:\n", "1994_3.8_36.6_3_stoPET.txt and 1994_3.8_36.6_3_AdjstoPET.txt**\n", "\n", "stoPET output are hourly PET values of the year which has 8760 hours for normal year and 8784 for leap years. Hence, in procesing the values notice these array length for each year." ] }, { "cell_type": "markdown", "id": "2eebcfff", "metadata": {}, "source": [ "### 1.4 Post processing and visualization\n", "Here we have a basic script to analyse the generated PET data and visualize it through simple plots. The script is named as **post_processing_stopet.py** and it consists of functions that helps to process the hourly data and plot the values. The following functions are available currently but you can add any additional functions you would like to have.\n", "\n", " 1. leap_remove(timeseries) \n", " 2. running_mean(timeseries, n)\n", " 3. aggregate_data(timeseries, period)\n", " 4. timeseries_plot(data, xlabel, ylabel, title, plotpath, fname)\n", " 5. comparison_timeseries_plot(data_1, data_2, label_1, label_2, xlabel, ylabel, title, plotpath, fname)\n", " 6. comparison_density_plot(data_1, data_2, label_1, label_2, xlabel, ylabel, title, plotpath, fname)\n", " 7. plot_spatial(data, lats, lons, cmap , title, cbar_label, climin, climax, plotpath, figfname)\n", "\n", "You can get the details of the function input parameters by using the folowing help function after importing the script.\n", "**help(script.function)** \n", "\n", "**Example: help(leap_remove)**\n", "\n", "**Note:**\n", "If the user wants multi year values comparison one must write a simple script to append the time series of each year\n", "before using ploting function. " ] }, { "cell_type": "code", "execution_count": null, "id": "f14b2e17", "metadata": {}, "outputs": [], "source": [ "from post_processing_stopet import *\n", "help(leap_remove)" ] }, { "cell_type": "markdown", "id": "388e9dee", "metadata": {}, "source": [ "**Example:** let us remove the leap year data from 1996 timeseries\n", "first read the 1996 PET data and check the array length" ] }, { "cell_type": "markdown", "id": "1d873308", "metadata": {}, "source": [ "# Kenya 1996 data doesn't exist" ] }, { "cell_type": "code", "execution_count": null, "id": "a33a558f", "metadata": {}, "outputs": [], "source": [ "data = np.genfromtxt('results/Kenya_E0_StoPET/1996_3.8_36.6_3_stoPET.txt')\n", "print(len(data))\n", "new_data = leap_remove(data)\n", "print(len(new_data))" ] }, { "cell_type": "markdown", "id": "0cfa6df4", "metadata": {}, "source": [ "As seen above the origional data has a length of 8784 which is 366 days of hourly data. After we use the remove leap year function the new data has a length of 8760 which is 365 days of hourly values." ] }, { "cell_type": "markdown", "id": "8af673cd", "metadata": {}, "source": [ "### Example \n", "Now lets do a simple exercise based on a 10 year data for a single location in Eastern Kenya (Wajir).\n", "\n", " * lat = 1.73\n", " * lon = 40.09\n", " * start year = 2000\n", " * end year = 2010\n", " * two ensembles\n", " * using Method 3 for temperature adjustment\n", " * temperature increase of 1.5 degrees\n", "\n", "Adjust the **run_stoPET.py** acordingly and run to generate the PET data" ] }, { "cell_type": "code", "execution_count": null, "id": "99506dbe", "metadata": {}, "outputs": [], "source": [ "def run_stoPET():\n", " ## ----- CHANGE THE INPUT VARIABLES HERE -----##\n", " datapath = 'stopet_parameter_files/'\n", " outputpath = 'results/' \n", " runtype = 'single' #'regional' #\n", " startyear = 2000\n", " endyear = 2010\n", "\n", " # Single point stoPET run\n", " latval = 1.73\n", " lonval = 40.09\n", "\n", " # Regional stoPET run\n", " latval_min = -5.5\n", " latval_max = -4.5 #5.5\n", " lonval_min = 33.0\n", " lonval_max = 34.5 #42.0\n", " locname = 'Wajir' #'Turkana1' #\n", "\n", " number_ensm = 2\n", " tempAdj = 3\n", " deltat = 1.5\n", " udpi_pet = 5\n", "\n", " ## ------ NO CHANGES BELLOW THIS -------------##\n", " if runtype == 'single':\n", " for ens_num in np.arange(0,number_ensm):\n", " stoPET_wrapper_singlepoint(startyear, endyear, \n", " latval, lonval, locname,\n", " ens_num,datapath, outputpath, \n", " tempAdj, deltat, udpi_pet)\n", " elif runtype == 'regional':\n", " for ens_num in np.arange(0,number_ensm):\n", " stoPET_wrapper_regional(startyear, endyear, \n", " latval_min, latval_max, \n", " lonval_min, lonval_max,\n", " locname, ens_num, datapath, \n", " outputpath, \n", " tempAdj, deltat, udpi_pet)\n", " else:\n", " raise ValueError('runtype only takes single and regional ... please check!')\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f21cb925", "metadata": {}, "outputs": [], "source": [ "run_stoPET()" ] }, { "cell_type": "markdown", "id": "be625eb1", "metadata": {}, "source": [ "Now let us make the visulization of the daily time series of PET for the first year (2000). Notice 2000 is a leap year so first we need to remove the leap year values since it won't be necessary for now" ] }, { "cell_type": "code", "execution_count": null, "id": "37415a3e", "metadata": {}, "outputs": [], "source": [ "# read the PET data for the year without temperature adjustment\n", "data_1 = np.genfromtxt('results/Wajir_E0_StoPET/2000_1.73_40.09_3_stoPET.txt')\n", "# read the PET data for the year with adjusted temeperature\n", "data_2 = np.genfromtxt('results/Wajir_E0_StoPET/2009_1.73_40.09_3_AdjstoPET.txt')\n", "# lets remove the leap year day data\n", "data_1 = leap_remove(data_1)\n", "data_2 = leap_remove(data_2)\n", "# check the length of the array\n", "print(len(data_1))\n", "print(len(data_2))" ] }, { "cell_type": "code", "execution_count": null, "id": "d201a00e", "metadata": {}, "outputs": [], "source": [ "# Density plot\n", "label_1 = 'stoPET'\n", "label_2 = 'Adj. stoPET'\n", "xlabel = 'PET ($\\mathbf{mm\\,hour^{-1}}$)' \n", "ylabel = 'Density' \n", "title = 'Hourly PET value' \n", "plotpath = './plots/' \n", "fig=plt.figure()\n", "# Draw the density plot\n", "sns.kdeplot(data_1, color = 'k',label = label_1)\n", "# Draw the density plot\n", "sns.kdeplot(data_2, color = 'orange',label = label_2)\n", "plt.xlabel(xlabel)\n", "plt.ylabel(ylabel)\n", "plt.title(title,fontweight='bold') \n", "plt.legend(loc='best')\n", "plt.tight_layout() " ] }, { "cell_type": "code", "execution_count": null, "id": "e0871bb1", "metadata": {}, "outputs": [], "source": [ "# Now lets plot the timeseres for both PET \n", "# values to visulize the daily timeseris of \n", "# non-adjusted and adjusted PET \n", "label_1 = 'stoPET'\n", "label_2 = 'Adj. stoPET'\n", "xlabel = 'PET ($\\mathbf{mm\\,hour^{-1}}$)' \n", "ylabel = 'Density' \n", "title = 'Hourly PET value' \n", "plotpath = './plots/' \n", "fname = 'Wajir_2000_hourly_density_stoPET.png'\n", "# by using the function for plotting from the \n", "# post_processing_stopet.py the figures will be \n", "# saved in the folder user provided.\n", "comparison_density_plot(data_1, data_2, label_1, label_2, \n", " xlabel, ylabel, title, plotpath, fname)" ] }, { "cell_type": "markdown", "id": "4beed528", "metadata": {}, "source": [ "Now lets make the daily and monthly aggregate for the two datasets using the **aggregate_data(timeseries, period)** function." ] }, { "cell_type": "code", "execution_count": null, "id": "99ce096f", "metadata": {}, "outputs": [], "source": [ "# daily\n", "data_1_day = aggregate_data(data_1, 'day')\n", "data_2_day = aggregate_data(data_2, 'day')\n", "print(len(data_1_day))\n", "print(len(data_2_day))\n", "# check for the array length\n", "# monthly\n", "data_1_month = aggregate_data(data_1, 'month')\n", "data_2_month = aggregate_data(data_2, 'month')\n", "# check for the array length\n", "print(len(data_1_month))\n", "print(len(data_2_month))" ] }, { "cell_type": "markdown", "id": "9bd06a41", "metadata": {}, "source": [ "Now lets plot the timeseres for both PET values to visulize the daily timeseris of non-adjusted and adjusted PET. " ] }, { "cell_type": "code", "execution_count": null, "id": "42ac1705", "metadata": {}, "outputs": [], "source": [ "data = data_1_day\n", "xlabel = 'Day of year'\n", "ylabel = 'PET ($\\mathbf{mm\\,day^{-1}}$)' \n", "title = 'Daily PET value' \n", "plotpath = './plots/' \n", "fname = 'Wajir_2000_daily_stoPET.png'\n", "timeseries_plot(data, xlabel, ylabel, title, plotpath, fname)" ] }, { "cell_type": "code", "execution_count": null, "id": "09180d02", "metadata": {}, "outputs": [], "source": [ "# This is just to show the plots in the document the function \n", "# already saved the plot in the plots folder.\n", "fig=plt.figure()\n", "plt.plot(data,'k')\n", "plt.ylabel(ylabel) \n", "plt.xlabel(xlabel) \n", "plt.title(title,fontweight='bold') \n", "plt.tight_layout() " ] }, { "cell_type": "code", "execution_count": null, "id": "92bc1c25", "metadata": {}, "outputs": [], "source": [ "# Now lets plot the timeseres for the monthly PET values of the year 2000.\n", "data = data_1_month\n", "xlabel = 'Month'\n", "ylabel = 'PET ($\\mathbf{mm\\,month^{-1}}$)' \n", "title = 'Monthly PET value' \n", "plotpath = 'plots/' \n", "fname = 'Wajir_2000_monthly_stoPET.png'\n", "timeseries_plot(data, xlabel, ylabel, title, plotpath, fname)" ] }, { "cell_type": "code", "execution_count": null, "id": "957957b0", "metadata": {}, "outputs": [], "source": [ "# This is just to show the plots in the document the function \n", "# already saved the plot in the plots folder.\n", "fig=plt.figure()\n", "plt.plot(data,'k')\n", "plt.ylabel(ylabel) \n", "plt.xlabel(xlabel) \n", "plt.title(title,fontweight='bold') \n", "plt.tight_layout() " ] }, { "cell_type": "markdown", "id": "0536afcf", "metadata": {}, "source": [ "Now lets plot the annual PET values and compare the temperature adjusted PET with non adjusted PET. Remember stoPET provide one file for each year so we need to loop through each year data and concatenate the values before plotting." ] }, { "cell_type": "code", "execution_count": null, "id": "b8821bc4", "metadata": {}, "outputs": [], "source": [ "# first lets make the year list\n", "years = np.arange(2000,2011)\n", "# create empty array\n", "data_1_year = [] # this is non adjusted PET\n", "data_2_year = [] # this is the adjusted PET\n", "# make a loop and estimate the annual PET value for both data\n", "for i in range (0, len(years)):\n", " year = years[i]\n", " data_1 = np.genfromtxt('results/Wajir_E0_StoPET/%s_1.73_40.09_3_stoPET.txt'\n", " %year)\n", " data_2 = np.genfromtxt('results/Wajir_E0_StoPET/%s_1.73_40.09_3_AdjstoPET.txt'\n", " %year)\n", " # make the annual aggregate\n", " data_1_val = aggregate_data(data_1, 'year')\n", " data_2_val = aggregate_data(data_2, 'year')\n", " # append the data to the empty array\n", " data_1_year = np.append(data_1_year, data_1_val)\n", " data_2_year = np.append(data_2_year, data_2_val)\n", " print(year)" ] }, { "cell_type": "code", "execution_count": null, "id": "ebe8cadf", "metadata": {}, "outputs": [], "source": [ "# now let us plot the two datasets for comparison\n", "data_1 = data_1_year\n", "data_2 = data_2_year\n", "label_1 = 'stoPET'\n", "label_2 = 'Adj. stoPET'\n", "xlabel = 'Year'\n", "ylabel = 'PET ($\\mathbf{mm\\,year^{-1}}$)' \n", "title = 'Annual PET value' \n", "plotpath = 'plots/' \n", "fname = 'Wajir_annual_stoPET_comparison.png'\n", "comparison_timeseries_plot(data_1, data_2, label_1, label_2, \n", " xlabel, ylabel, title, plotpath, fname)" ] }, { "cell_type": "code", "execution_count": null, "id": "56c23512", "metadata": {}, "outputs": [], "source": [ "# This is just to show the plots in the document the function \n", "# already saved the plot in the plots folder.\n", "fig=plt.figure()\n", "plt.plot(data_1,'k', label=label_1)\n", "plt.plot(data_2,'k--', label=label_2)\n", "plt.ylabel(ylabel) \n", "plt.xlabel(xlabel) \n", "plt.title(title,fontweight='bold') \n", "plt.legend(loc='best')\n", "plt.tight_layout() " ] }, { "cell_type": "markdown", "id": "e113e86b", "metadata": {}, "source": [ "---\n", "## EXERCISE \n", "Based on the above examples please try to do the folowing:\n", "\n", " 1. Generate a 15 year dataset for the location of your choice.\n", " 2. Use each method of temperature adjustment (Method 1, Method 2 and Method 3) separatelly.\n", " 3. Generate 3 plots (one for each method) comparing the values of the annual PET and adjusted PET.\n", "--- " ] }, { "cell_type": "code", "execution_count": null, "id": "ddcfe71e", "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.8.16" } }, "nbformat": 4, "nbformat_minor": 5 }