{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "#### Contents" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following course cover the following content:\n", "\n", "* **Convert climate seasonal forecasting into model forcing datasets**\n", "* Translate climate forecasting into hydrological forecasting\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4. Convert climate seasonal forecasting into model forcing datasets" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Create an ensamble of precipitation for the Tana basin\n", "* Create ensamble of potential evapotranspiration for the Tana basin" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 4.1. Create an ensamble of precipitation for the Tana basin" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Create stochastic realizations of precipitation for the MAM season using STORM\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 4.2. Create ensamble of potential evapotranspiration for the Tana basin" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Create stochastic realizations of potential evapotranspiration for the MAM season using stoPET.\n", "Forecasting the seasonal Potential Evapotranspiration (PET) values uses the stoPET model and the ICPAC seasonal temperature forecast (ICPAC seasonal forecast).\n", "\n", "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", "Here, we are assuming that PET is related to temperature. In areas where the seasonal temperature is expected to be high, we also expect higher PET and vice versa. \n", "\n", "# How does it work?\n", "The seasonal PET forecast works based on the PET that will be generated using the stochastic PET generator model (stoPET). This model allows PET to be generated based on increasing or decreasing temperature. \n", "\n", "Since the ICPAC seasonal temperature forecasts are provided in tercile probabilistic terms, we need to follow a similar approach in producing the PET that reflects the seasonal probabilistic temperature forecast. \n", "\n", "We need to generate an ensemble of PET values to represent all possibilities of the future seasonal PET based on the temperature forecast. Hence, choosing a minimum of 30 ensembles is preferable. \n", "\n", "To do this, we follow a two-step approach: \n", "\n", " 1. Generate a pool of PET values for each of the tercile categories based on changes in high, medium, and low temperatures.\n", " \n", " Here, we use 1.5 degrees for the above-average (A), 0.0 degrees for the normal (N) and \n", " -1.5 degrees for the below-average (B) tercile categories. \n", " \n", " We generate an equal number of ensembles for each tercile category. (e.g. A = 30, N=30, B=30)\n", " \n", " 2. For each grid cell, we will choose a similar percentage of PET values from the pools generated based on what is \n", " provided as a seasonal probability for temperature from ICPAC. \n", " \n", " Here, we will take the ICPAC probabilities of Above, Normal, and Below and identify the number of ensembles in each \n", " category and choose PET values from each category based on the result. \n", " \n", " Example:\n", " If we have 30 ensembles of PET and the ICPAC temperature forecast is given as A=50, N=30, B=20, then we take \n", " A => 30 * 0.5 = 15 values\n", " N => 30 * 0.3 = 9 values\n", " B => 30 * 0.2 = 6 values\n", " from this, we will get ensembles of forecast PET where 15 ensembles are above average, 9 ensembles are normal, and \n", " 6 ensembles are below average.\n", " \n", "# How to do it?\n", "To perform all the above steps and generate a seasonal PET forecast, we need the stoPET model (which will be used for step 1) and the forecast generation model (which will be used in step 2).\n", "\n", "Four Python scripts, four parameter NetCDF files, and one NetCDF file containing the ICPAC seasonal temperature forecast are required.\n", "\n", "Python scripts:\n", "\n", " a. stoPET_v1_4dryp.py\n", " b. config_stoPET.py\n", " c. run_stoPET_4dryp.py\n", " d. forecast_generation.py\n", " \n", "Parameter files:\n", "\n", " a. dpetdt.nc\n", " b. hpet_slope.nc\n", " c. monthly_cont_percentage.nc\n", " d. stopet_parameters.nc\n", "\n", "ICPAC seasonal temperature forecast:\n", "\n", " a. TanaBasin_TempF_MAM2022.nc\n", "\n", "Put all the Python scripts in one folder and the parameter files separately in another folder. \n", "Then open the **config_stoPET.py** script and provide all the adjustments you want to make in that script. This script is where you provide 20 variables required to run the stoPET model and forecast generation. \n", "\n", "The **config_stoPET.py** looks like the following. Basically, all 20 parameters are self-explanatory. You can make any change you want based on each variable. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# How to edit the Config file\n", "\n", "This config file (config_stoPET.py) is used to provide all the required \n", "data and parameters to run the stoPET model and generate the \n", "required seasonal forecasts accounting for the ICPAC temperature forecast\n", "given as a seasonal probability. \n", "\n", "This will be used in running the **run_stoPET.py** to generate multiple files\n", "of forecast pools with above-average, normal and below-average forecasts.\n", "Finally, these generated PET values will be used in the **forecast_generation.py** to produce the final PET \n", "forecasts for the season accounting ICPAC temperature forecast.\n", "\n", "**CHANGE THE INPUT VARIABLES HERE** \n", "\n", "1. This is the data path where the stoPET parameter files are kept.\n", "\n", "**datapath = '/user/work/fp20123/stoPET/stopet_parameter_files/'**\n", "\n", "2. This is the output path where you want to keep the generated PET files\n", "\n", "**outputpath = '/user/home/fp20123/new_stopet/result/'**\n", "\n", "3. This is where you decide whether to run a 'regional' model or a 'single' point model.\n", "\n", "**runtype = 'regional'**\n", "\n", "4. The number of years to generate data for\n", "\n", "**startyear = 2022**\n", "\n", "**endyear = 2022**\n", " \n", "5. These are the seasonal Julian dates required as start and end dates.\n", "\n", "**seasonswitch = 1** \n", "\n", "**startdate = 60** \n", "\n", "**enddate = 151**\n", "\n", "**seasonName = 'MAM'**\n", "\n", "6. Single-point stoPET run\n", "\n", "**latval = 1.0**\n", "\n", "**lonval = 35.0**\n", "\n", "7. Regional stoPET run\n", "\n", "**latval_min = -1.6**\n", "\n", "**latval_max = 1.1**\n", "\n", "**lonval_min = 35.9**\n", "\n", "**lonval_max = 40.1** \n", "\n", "8. This is the number of ensembles you want to generate.\n", " \n", "**number_ensm = 10**\n", "\n", "9. This is the method you use for adjusting the PET to account for temperature changes.\n", "\n", "**tempAdj = 2**\n", "\n", "10. This is the change in temperature you prefer to have if you use tempAdj = 2\n", "\n", "**deltat = 0.0**\n", "\n", "11. This is where you provide a sting of name to identify your region or point, e.g. 'tana_basin'.\n", " \n", "**locname = 'TanaBasin'** \n", "\n", "12. This is the user-defined percentage change in PET if you use tempAdj = 1\n", "\n", "**udpi_pet = 5**\n", "\n", "13. This is the file containing the ICPAC seasonal tercile temperature forecast.\n", "Provide the full path where it is located.\n", "\n", "**tercile_forecast_file = '/user/home/fp20123/new_stopet/TanaBasin_TempF_MAM2022.nc'**\n", "\n", "14. This is the file containing masks required for plotting and filling gaps\n", "\n", "**tercile_forecast_file = '/user/home/fp20123/new_stopet/TanaBasin_mask.nc'**" ] }, { "cell_type": "markdown", "metadata": { "vscode": { "languageId": "plaintext" } }, "source": [ "# Generating PET forecast\n", "\n", "Once you make all the required changes to the **config_stoPET.py** file, you can run the **run_stoPET_4dryp.py** script, which will generate the PET ensembles to be used as a pool. Here, you have to repeat the process by changing the **deltat** variable to 1.5 for Above, 0.0 for Normal, and -1.5 for Below categories.\n", "\n", "This means you modify the **config_stoPET.py** file and run **run_stoPET_4dryp.py**, and when one is finished, you repeat the process by modifying the **deltat** accordingly. \n", "\n", "## NB\n", "### Please restart the Kernel before running the code, as it saves all the variables from the previous run." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "# This is the path where you put all the python scripts and the required input data files\n", "# path = 'C:/Users/fp20123/Dropbox/CUWALID/new_stopet/'\n", "path = '/home/cuwalid/CUWALID/CUWALID_training/stoPET/'" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "import os\n", "# create a folder to save the data\n", "if not os.path.isdir(path + 'result/'):\n", " os.mkdir(path + 'result/')\n", "\n", "# create a folder to save plots\n", "if not os.path.isdir(path + 'plots/'):\n", " os.mkdir(path + 'plots/')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "# Let's run 10 ensembles for the Above-average category pool\n", "# make sure you change deltat = 1.5 in the config_stoPET.py\n", "%run /home/cuwalid/CUWALID/CUWALID_training/stoPET/run_stoPET_4dryp.py" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "# Let's run 10 ensembles for the Below-average category pool\n", "# make sure you change deltat = -1.5 in the config_stoPET.py\n", "%run /home/cuwalid/CUWALID/CUWALID_training/stoPET/run_stoPET_4dryp.py" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "# Let's run 10 ensembles for the Normal category pool\n", "# make sure you change deltat = 0.0 in the config_stoPET.py\n", "%run /home/cuwalid/CUWALID/CUWALID_training/stoPET/run_stoPET_4dryp.py" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "# Let's run 10 ensemble forecasts of PET\n", "%run /home/cuwalid/CUWALID/CUWALID_training/stoPET/forecast_generation.py" ] }, { "cell_type": "markdown", "metadata": { "vscode": { "languageId": "plaintext" } }, "source": [ "## Visualization (Plotting PET output)\n", "Let's plot the values of the PET for visualization." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "import warnings\n", "warnings.filterwarnings(\"ignore\", \"is_categorical_dtype\")\n", "warnings.filterwarnings(\"ignore\", \"use_inf_as_na\")\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from netCDF4 import Dataset, num2date\n", "from mpl_toolkits.basemap import Basemap, maskoceans\n", "import seaborn as sns\n", "# Import libraries from local repository\n", "import sys\n", "\n", "#sys.path.append('/home//CUWALID/CUWALID_training/stoPET/')\n", "sys.path.append('/home/aquichimbo/CUWALID/CUWALID_training/stoPET/')\n", "from post_processing_stopet import *\n", "#from home.aquichimbo.CUWALID.CUWALID_training.stoPET.post_processing_stopet import *" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from h" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "# This is a plotting function \n", "def trend_plot(data, lats, lons, title, label, vmin, vmax, fname, creverse, title2, shapefile):\n", " \n", " if creverse == 0:\n", " cmap = plt.cm.bwr_r #Spectral\n", " else:\n", " cmap = plt.cm.bwr #Spectral_r\n", " \n", " fig = plt.figure(figsize=(6, 6))\n", " ax=fig.add_axes([0.1,0.1,0.8,0.8])\n", " m = Basemap(projection='cyl', llcrnrlat=min(lats), urcrnrlat=max(lats), llcrnrlon=min(lons),\n", " urcrnrlon=max(lons), resolution='c')\n", " cs4 = plt.imshow(data[:,:], interpolation='nearest', cmap=cmap,\n", " extent=[min(lons), max(lons), min(lats), max(lats)],vmin=vmin,vmax=vmax) \n", " m.drawcoastlines(linewidth=0.75)\n", "# m.drawcountries(linewidth=0.75)\n", " parallels=np.arange(-90.,90.,0.8)\n", " meridians=np.arange(-180.,180.,0.8)\n", " m.drawparallels(parallels,labels=[True, False, False, False], size=8)\n", " m.drawmeridians(meridians,labels=[False, False, False, True], size=8)\n", " cb4 = plt.colorbar(cs4, label=label, shrink=0.75, pad=0.05, extend='both', orientation='horizontal')#ticks=tickz,\n", " m.readshapefile(shapefile,'tana_basins', linewidth=2) \n", " plt.title(title, fontweight='bold')\n", " plt.title(title2,loc='left', fontweight='bold')\n", "# plt.tight_layout() \n", "# fig.savefig(fname,bbox_inches='tight', dpi=300)\n", "# plt.close() " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "# read the data\n", "filename = 'result/ensemble_forecast/Forecast_PET_TanaBasin_ens_0_MAM_2022.nc'\n", "\n", "nc = Dataset(path + filename)\n", "lats = nc.variables['latitude'][:]\n", "lons = nc.variables['longitude'][:]\n", "pet = nc.variables['pet'][:,:,:]\n", "\n", "# This checks the shape of the data \n", "print(pet.shape)\n", "\n", "# read the data\n", "filename = 'result/ensemble_forecast/Forecast_PET_TanaBasin_ens_6_MAM_2022.nc'\n", "\n", "nc = Dataset(path + filename)\n", "lats = nc.variables['latitude'][:]\n", "lons = nc.variables['longitude'][:]\n", "pet2 = nc.variables['pet'][:,:,:]\n", "\n", "# This checks the shape of the data \n", "print(pet2.shape)\n", "\n", "# read the data\n", "filename = 'result/ensemble_forecast/Forecast_PET_TanaBasin_ens_9_MAM_2022.nc'\n", "\n", "nc = Dataset(path + filename)\n", "lats = nc.variables['latitude'][:]\n", "lons = nc.variables['longitude'][:]\n", "pet3 = nc.variables['pet'][:,:,:]\n", "\n", "# This checks the shape of the data \n", "print(pet3.shape)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "# let's plot hourly PET\n", "# input for plot\n", "plotpath = path + 'plots/'\n", "shapefile = path + 'shapefiles/Tana_basin_wgs84' # Kenya\n", "#shapefile = path + 'shapefiles/UAVBasin' # Ethiopia\n", "#shapefile = path + 'shapefiles/OD_domain' # Somalia\n", "# ------------------------------- #\n", "hour = 12\n", "vmin=0.0\n", "vmax=1.0\n", "label = '$\\mathrm{mm\\,h^{-1}}$'\n", "title = 'hourly PET (hour-%s)'%(hour+1)\n", "fname = plotpath + 'hourlyPET_plot_%s.png'%(hour+1)\n", "\n", "data1 = pet[hour,:,:]\n", "trend_plot(data1, lats, lons, title, label, vmin, vmax, fname, 1, 'a)', shapefile)\n", "\n", "data2 = pet2[hour,:,:]\n", "trend_plot(data2, lats, lons, title, label, vmin, vmax, fname, 1, 'b)', shapefile)\n", "\n", "data3 = pet3[hour,:,:]\n", "trend_plot(data3, lats, lons, title, label, vmin, vmax, fname, 1, 'c)', shapefile)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "# let's plot seasonal PET\n", "# input for plot\n", "plotpath = path + 'plots/'\n", "shapefile = path + 'shapefiles/Tana_basin_wgs84' # Kenya\n", "#shapefile = path + 'shapefiles/UAVBasin' # Ethiopia\n", "#shapefile = path + 'shapefiles/OD_domain' # Somalia\n", "# -------------------------------------- #\n", "\n", "\n", "vmin=300.0\n", "vmax=500.0\n", "label = '$\\mathrm{mm}$'\n", "title = 'Seasonal PET (MAM)'\n", "fname = plotpath + 'seasonalPET_plot.png'\n", "\n", "data1 = np.sum(pet[:,:,:], axis=0)\n", "trend_plot(data1, lats, lons, title, label, vmin, vmax, fname, 1, 'a)', shapefile)\n", "\n", "data2 = np.sum(pet2[:,:,:], axis=0)\n", "trend_plot(data2, lats, lons, title, label, vmin, vmax, fname, 1, 'b)', shapefile)\n", "\n", "data3 = np.sum(pet3[:,:,:], axis=0)\n", "trend_plot(data3, lats, lons, title, label, vmin, vmax, fname, 1, 'b)', shapefile)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "# This function plots daily PET density plots.\n", "# The plot shows how the PET values include the ICPAC forecast\n", "\n", "def density_plot_seasonal(data, title):\n", " for i in range(data.shape[0]):\n", " # Draw the density plot\n", " if i == 0:\n", " sns.kdeplot(data[i,:], lw=3, color = 'grey', label = 'ensem-PET')\n", " else:\n", " sns.kdeplot(data[i,:], lw=3, color = 'grey')\n", " \n", " # Plot formatting\n", " plt.legend(prop={'size': 10})\n", " plt.title(title, weight='bold')\n", " plt.xlabel('$\\mathrm{PET\\,(mm)}$')\n", " plt.ylabel('Density')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "# read the data\n", "# choose the point you want to plot\n", "lat_var = 0.4\n", "lon_var = 39.2\n", "filename = 'TanaBasin_TempF_MAM2022.nc'\n", "nc = Dataset(path + filename)\n", "lats = nc.variables['lat'][:]\n", "lons = nc.variables['lon'][:]\n", "above = nc.variables['above'][:,:]\n", "normal = nc.variables['normal'][:,:]\n", "below = nc.variables['below'][:,:]\n", "lati, loni = nearest_point(lat_var, lon_var, lats, lons)\n", "print(lati,loni)\n", "print('Above = %s'%(above[lati,loni]))\n", "print('Normal = %s'%(normal[lati,loni]))\n", "print('Below = %s'%(below[lati,loni]))\n", "\n", "\n", "data = []\n", "for i in range(0,10):\n", " filename = 'result/ensemble_forecast/Forecast_PET_TanaBasin_ens_%s_MAM_2022.nc'%i\n", " nc = Dataset(path + filename)\n", " lats = nc.variables['latitude'][:]\n", " lons = nc.variables['longitude'][:]\n", " pet = nc.variables['pet'][:,lati,loni]\n", " out = aggregate_data(pet, 'day')\n", " data.append(out)\n", "data = np.asarray(data)\n", "print(data.shape)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "title = 'Daily PET'\n", "density_plot_seasonal(data, title)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "# This is to plot the daily timeserise of an above normal and below normal PET\n", "xlabel='Days'\n", "ylabel='$\\mathrm{PET (mm\\,d^{-1})}$'\n", "title='Daily PET timeseries (MAM)'\n", "plotpath= path + 'plots/'\n", "fname=plotpath + 'daily_PET_plot.png'\n", "label_1 = 'above normal' \n", "label_2 = 'below normal'\n", "\n", "data_1 = np.array(data[0,:])\n", "data_2 = np.array(data[9,:])\n", "\n", "comparison_timeseries_plot(data_1, data_2, label_1, label_2, xlabel, ylabel, title, plotpath, fname)" ] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:cwld]", "language": "python", "name": "conda-env-cwld-py" }, "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.11.4" } }, "nbformat": 4, "nbformat_minor": 4 }