Contents

The following course cover the following content:

• Convert climate seasonal forecasting into model forcing datasets

• Translate climate forecasting into hydrological forecasting

4. Convert climate seasonal forecasting into model forcing datasets

1. Create an ensamble of precipitation for the Tana basin

• Create ensamble of potential evapotranspiration for the Tana basin

4.1. Create an ensamble of precipitation for the Tana basin

• Create stochastic realizations of precipitation for the MAM season using STORM

4.2. Create ensamble of potential evapotranspiration for the Tana basin

• Create stochastic realizations of potential evapotranspiration for the MAM season using stoPET. Forecasting the seasonal Potential Evapotranspiration (PET) values uses the stoPET model and the ICPAC seasonal temperature forecast (ICPAC seasonal forecast).

Details of the model description can be found in the following paper stoPET paper and you can download the model from this link.

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.

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.

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.

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.

To do this, we follow a two-step approach:

1. Generate a pool of PET values for each of the tercile categories based on changes in high, medium, and low temperatures.

    Here, we use 1.5 degrees for the above-average (A), 0.0 degrees for the normal (N) and
    -1.5 degrees for the below-average (B) tercile categories.

    We generate an equal number of ensembles for each tercile category. (e.g. A = 30, N=30, B=30)

2. For each grid cell, we will choose a similar percentage of PET values from the pools generated based on what is
    provided as a seasonal probability for temperature from ICPAC.

    Here, we will take the ICPAC probabilities of Above, Normal, and Below and identify the number of ensembles in each
    category and choose PET values from each category based on the result.

    Example:
    If we have 30 ensembles of PET and the ICPAC temperature forecast is given as A=50, N=30, B=20, then we take
        A => 30 * 0.5 = 15 values
        N => 30 * 0.3 = 9 values
        B => 30 * 0.2 = 6 values
    from this, we will get ensembles of forecast PET where 15 ensembles are above average, 9 ensembles are normal, and
    6 ensembles are below average.

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).

Four Python scripts, four parameter NetCDF files, and one NetCDF file containing the ICPAC seasonal temperature forecast are required.

Python scripts:

a. stoPET_v1_4dryp.py
b. config_stoPET.py
c. run_stoPET_4dryp.py
d. forecast_generation.py

Parameter files:

a. dpetdt.nc
b. hpet_slope.nc
c. monthly_cont_percentage.nc
d. stopet_parameters.nc

ICPAC seasonal temperature forecast:

a. TanaBasin_TempF_MAM2022.nc

Put all the Python scripts in one folder and the parameter files separately in another folder. 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.

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.

This config file (config_stoPET.py) is used to provide all the required data and parameters to run the stoPET model and generate the required seasonal forecasts accounting for the ICPAC temperature forecast given as a seasonal probability.

This will be used in running the run_stoPET.py to generate multiple files of forecast pools with above-average, normal and below-average forecasts. Finally, these generated PET values will be used in the forecast_generation.py to produce the final PET forecasts for the season accounting ICPAC temperature forecast.

CHANGE THE INPUT VARIABLES HERE

1. This is the data path where the stoPET parameter files are kept.

datapath = ‘/user/work/fp20123/stoPET/stopet_parameter_files/’

2. This is the output path where you want to keep the generated PET files

outputpath = ‘/user/home/fp20123/new_stopet/result/’

3. This is where you decide whether to run a ‘regional’ model or a ‘single’ point model.

runtype = ‘regional’

4. The number of years to generate data for

startyear = 2022

endyear = 2022

5. These are the seasonal Julian dates required as start and end dates.

seasonswitch = 1

startdate = 60

enddate = 151

seasonName = ‘MAM’

6. Single-point stoPET run

latval = 1.0

lonval = 35.0

7. Regional stoPET run

latval_min = -1.6

latval_max = 1.1

lonval_min = 35.9

lonval_max = 40.1

8. This is the number of ensembles you want to generate.

number_ensm = 10

9. This is the method you use for adjusting the PET to account for temperature changes.

tempAdj = 2

10. This is the change in temperature you prefer to have if you use tempAdj = 2

deltat = 0.0

11. This is where you provide a sting of name to identify your region or point, e.g. ‘tana_basin’.

locname = ‘TanaBasin’

12. This is the user-defined percentage change in PET if you use tempAdj = 1

udpi_pet = 5

13. This is the file containing the ICPAC seasonal tercile temperature forecast. Provide the full path where it is located.

tercile_forecast_file = ‘/user/home/fp20123/new_stopet/TanaBasin_TempF_MAM2022.nc’

14. This is the file containing masks required for plotting and filling gaps

tercile_forecast_file = ‘/user/home/fp20123/new_stopet/TanaBasin_mask.nc’

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.

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.

Please restart the Kernel before running the code, as it saves all the variables from the previous run.

# This is the path where you put all the python scripts and the required input data files
# path = 'C:/Users/fp20123/Dropbox/CUWALID/new_stopet/'
path = '/home/cuwalid/CUWALID/CUWALID_training/stoPET/'

import os
# create a folder to save the data
if not os.path.isdir(path + 'result/'):
    os.mkdir(path + 'result/')

# create a folder to save plots
if not os.path.isdir(path + 'plots/'):
    os.mkdir(path + 'plots/')

# Let's run 10 ensembles for the Above-average category pool
# make sure you change deltat = 1.5 in the config_stoPET.py
%run /home/cuwalid/CUWALID/CUWALID_training/stoPET/run_stoPET_4dryp.py

# Let's run 10 ensembles for the Below-average category pool
# make sure you change deltat = -1.5 in the config_stoPET.py
%run /home/cuwalid/CUWALID/CUWALID_training/stoPET/run_stoPET_4dryp.py

# Let's run 10 ensembles for the Normal category pool
# make sure you change deltat = 0.0 in the config_stoPET.py
%run /home/cuwalid/CUWALID/CUWALID_training/stoPET/run_stoPET_4dryp.py

# Let's run 10 ensemble forecasts of PET
%run /home/cuwalid/CUWALID/CUWALID_training/stoPET/forecast_generation.py

Visualization (Plotting PET output)

Let’s plot the values of the PET for visualization.

import warnings
warnings.filterwarnings("ignore", "is_categorical_dtype")
warnings.filterwarnings("ignore", "use_inf_as_na")

import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset, num2date
from mpl_toolkits.basemap import Basemap, maskoceans
import seaborn as sns
# Import libraries from local repository
import sys

#sys.path.append('/home/<username>/CUWALID/CUWALID_training/stoPET/')
sys.path.append('/home/aquichimbo/CUWALID/CUWALID_training/stoPET/')
from post_processing_stopet import *
#from home.aquichimbo.CUWALID.CUWALID_training.stoPET.post_processing_stopet import *

from h

# This is a plotting function
def trend_plot(data, lats, lons, title, label, vmin, vmax, fname, creverse, title2, shapefile):

    if creverse == 0:
      cmap = plt.cm.bwr_r  #Spectral
    else:
      cmap = plt.cm.bwr  #Spectral_r

    fig = plt.figure(figsize=(6, 6))
    ax=fig.add_axes([0.1,0.1,0.8,0.8])
    m = Basemap(projection='cyl', llcrnrlat=min(lats), urcrnrlat=max(lats), llcrnrlon=min(lons),
                urcrnrlon=max(lons), resolution='c')
    cs4 = plt.imshow(data[:,:], interpolation='nearest', cmap=cmap,
                     extent=[min(lons), max(lons), min(lats), max(lats)],vmin=vmin,vmax=vmax)
    m.drawcoastlines(linewidth=0.75)
#    m.drawcountries(linewidth=0.75)
    parallels=np.arange(-90.,90.,0.8)
    meridians=np.arange(-180.,180.,0.8)
    m.drawparallels(parallels,labels=[True, False, False, False], size=8)
    m.drawmeridians(meridians,labels=[False, False, False, True], size=8)
    cb4 = plt.colorbar(cs4, label=label, shrink=0.75, pad=0.05, extend='both', orientation='horizontal')#ticks=tickz,
    m.readshapefile(shapefile,'tana_basins', linewidth=2)
    plt.title(title, fontweight='bold')
    plt.title(title2,loc='left', fontweight='bold')
#    plt.tight_layout()
#    fig.savefig(fname,bbox_inches='tight', dpi=300)
#    plt.close()

# read the data
filename = 'result/ensemble_forecast/Forecast_PET_TanaBasin_ens_0_MAM_2022.nc'

nc = Dataset(path + filename)
lats = nc.variables['latitude'][:]
lons = nc.variables['longitude'][:]
pet = nc.variables['pet'][:,:,:]

# This checks the shape of the data
print(pet.shape)

# read the data
filename = 'result/ensemble_forecast/Forecast_PET_TanaBasin_ens_6_MAM_2022.nc'

nc = Dataset(path + filename)
lats = nc.variables['latitude'][:]
lons = nc.variables['longitude'][:]
pet2 = nc.variables['pet'][:,:,:]

# This checks the shape of the data
print(pet2.shape)

# read the data
filename = 'result/ensemble_forecast/Forecast_PET_TanaBasin_ens_9_MAM_2022.nc'

nc = Dataset(path + filename)
lats = nc.variables['latitude'][:]
lons = nc.variables['longitude'][:]
pet3 = nc.variables['pet'][:,:,:]

# This checks the shape of the data
print(pet3.shape)

# let's plot hourly PET
# input for plot
plotpath = path + 'plots/'
shapefile = path + 'shapefiles/Tana_basin_wgs84'   # Kenya
#shapefile = path + 'shapefiles/UAVBasin'           # Ethiopia
#shapefile = path + 'shapefiles/OD_domain'           # Somalia
# ------------------------------- #
hour = 12
vmin=0.0
vmax=1.0
label = '$\mathrm{mm\,h^{-1}}$'
title = 'hourly PET (hour-%s)'%(hour+1)
fname = plotpath + 'hourlyPET_plot_%s.png'%(hour+1)

data1 = pet[hour,:,:]
trend_plot(data1, lats, lons, title, label, vmin, vmax, fname, 1, 'a)', shapefile)

data2 = pet2[hour,:,:]
trend_plot(data2, lats, lons, title, label, vmin, vmax, fname, 1, 'b)', shapefile)

data3 = pet3[hour,:,:]
trend_plot(data3, lats, lons, title, label, vmin, vmax, fname, 1, 'c)', shapefile)

# let's plot seasonal PET
# input for plot
plotpath = path + 'plots/'
shapefile = path + 'shapefiles/Tana_basin_wgs84'   # Kenya
#shapefile = path + 'shapefiles/UAVBasin'           # Ethiopia
#shapefile = path + 'shapefiles/OD_domain'           # Somalia
# -------------------------------------- #


vmin=300.0
vmax=500.0
label = '$\mathrm{mm}$'
title = 'Seasonal PET (MAM)'
fname = plotpath + 'seasonalPET_plot.png'

data1 = np.sum(pet[:,:,:], axis=0)
trend_plot(data1, lats, lons, title, label, vmin, vmax, fname, 1, 'a)', shapefile)

data2 = np.sum(pet2[:,:,:], axis=0)
trend_plot(data2, lats, lons, title, label, vmin, vmax, fname, 1, 'b)', shapefile)

data3 = np.sum(pet3[:,:,:], axis=0)
trend_plot(data3, lats, lons, title, label, vmin, vmax, fname, 1, 'b)', shapefile)

# This function plots daily PET density plots.
# The plot shows how the PET values include the ICPAC forecast

def density_plot_seasonal(data, title):
    for i in range(data.shape[0]):
        # Draw the density plot
        if i == 0:
            sns.kdeplot(data[i,:], lw=3, color = 'grey', label = 'ensem-PET')
        else:
            sns.kdeplot(data[i,:], lw=3, color = 'grey')

        # Plot formatting
        plt.legend(prop={'size': 10})
        plt.title(title, weight='bold')
        plt.xlabel('$\mathrm{PET\,(mm)}$')
        plt.ylabel('Density')

# read the data
# choose the point you want to plot
lat_var = 0.4
lon_var = 39.2
filename = 'TanaBasin_TempF_MAM2022.nc'
nc = Dataset(path + filename)
lats = nc.variables['lat'][:]
lons = nc.variables['lon'][:]
above = nc.variables['above'][:,:]
normal = nc.variables['normal'][:,:]
below = nc.variables['below'][:,:]
lati, loni = nearest_point(lat_var, lon_var, lats, lons)
print(lati,loni)
print('Above = %s'%(above[lati,loni]))
print('Normal = %s'%(normal[lati,loni]))
print('Below = %s'%(below[lati,loni]))


data = []
for i in range(0,10):
    filename = 'result/ensemble_forecast/Forecast_PET_TanaBasin_ens_%s_MAM_2022.nc'%i
    nc = Dataset(path + filename)
    lats = nc.variables['latitude'][:]
    lons = nc.variables['longitude'][:]
    pet = nc.variables['pet'][:,lati,loni]
    out = aggregate_data(pet, 'day')
    data.append(out)
data = np.asarray(data)
print(data.shape)

title = 'Daily PET'
density_plot_seasonal(data, title)

# This is to plot the daily timeserise of an above normal and below normal PET
xlabel='Days'
ylabel='$\mathrm{PET (mm\,d^{-1})}$'
title='Daily PET timeseries (MAM)'
plotpath= path + 'plots/'
fname=plotpath + 'daily_PET_plot.png'
label_1 = 'above normal'
label_2 = 'below normal'

data_1 = np.array(data[0,:])
data_2 = np.array(data[9,:])

comparison_timeseries_plot(data_1, data_2, label_1, label_2, xlabel, ylabel, title, plotpath, fname)