|
|
|
Learning objective:
Introducing stoPET model
learn the input parameters used in the model
changing parameter values
generating single point PET value
learn data visulaization from stoPET output
est. time 2 hour
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.
1. stoPET
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.
how hourly PET (hPET) was generated?
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 hPET paper.
Estimating PET using Penmann-Monteith method requires 7 input variables:
10 m u-component (zonal) of wind speed [m s−1]
10 m v-component (meridional) of wind speed [m s−1]
2 m dew point temperature [K]
2 m air temperature [K]
surface net solar radiation [J m−2]
surface net thermal radiation [J m−2]
atmospheric surface pressure [Pa]
The equation used is given as:
hPET = [0.408 * ∆ (Rn - G) + γ(37/Ta + 273)* u2(es - ea)] / [∆ + γ(1 + 0.34u2)]
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.
how stoPET was generated? 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.
Y = A sin (B × t + C) + D
where A represents the diurnal amplitude (mm h−1), B is thefrequency (h−1), C is the phase shift (–), and D is the vertical shift (mm h−1). t is time (h), and Y is the new PET value(mm h−1) generated from the sine function.
The overall equation is given as follows:
Stochastic PET = (average diurnal cycle of PET using a sine function × a random noise ratio) + user-defined annual PET variability.
It has three components:
the sine function
the random noise ratio
the adjustment for temperature.
[ ]:
# This is for showing plots
# interactively in Jupiter notebook
%matplotlib inline
[ ]:
import sys
import numpy as np
from stoPET_v1 import *
from mpl_toolkits.basemap import Basemap
[ ]:
def run_stoPET():
## ----- CHANGE THE INPUT VARIABLES HERE -----##
datapath = 'stopet_parameter_files/'
outputpath = 'results/'
runtype = 'single' #'regional' #
startyear = 2000
endyear = 2010
# Single point stoPET run
latval = 1.73
lonval = 40.09
# Regional stoPET run
latval_min = -5.5
latval_max = -4.5 #5.5
lonval_min = 33.0
lonval_max = 34.5 #42.0
locname = 'Wajir' #'Turkana1' #
number_ensm = 2
tempAdj = 3
deltat = 1.5
udpi_pet = 5
## ------ NO CHANGES BELLOW THIS -------------##
if runtype == 'single':
for ens_num in np.arange(0,number_ensm):
stoPET_wrapper_singlepoint(startyear, endyear,
latval, lonval, locname,
ens_num,datapath, outputpath,
tempAdj, deltat, udpi_pet)
elif runtype == 'regional':
for ens_num in np.arange(0,number_ensm):
stoPET_wrapper_regional(startyear, endyear,
latval_min, latval_max,
lonval_min, lonval_max,
locname, ens_num, datapath,
outputpath,
tempAdj, deltat, udpi_pet)
else:
raise ValueError('runtype only takes single and regional ... please check!')
1.1 Changing input parameters
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.
1.2 Adjusting for temperature increase
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).
1. Method 1: user-defined percentage step change in annual PET
2. Method 2: step change in PET based on a user-defined change in atmospheric temperature
3. Method 3: progressive change in PET based on the historical trend in hPET
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.
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.
[ ]:
run_stoPET()
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.
1.3 Model Output
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/.
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.
Example: 1994_3.8_36.6_3_stoPET.txt and 1994_3.8_36.6_3_AdjstoPET.txt
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.
1.4 Post processing and visualization
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.
1. leap_remove(timeseries)
2. running_mean(timeseries, n)
3. aggregate_data(timeseries, period)
4. timeseries_plot(data, xlabel, ylabel, title, plotpath, fname)
5. comparison_timeseries_plot(data_1, data_2, label_1, label_2, xlabel, ylabel, title, plotpath, fname)
6. comparison_density_plot(data_1, data_2, label_1, label_2, xlabel, ylabel, title, plotpath, fname)
7. plot_spatial(data, lats, lons, cmap , title, cbar_label, climin, climax, plotpath, figfname)
You can get the details of the function input parameters by using the folowing help function after importing the script. help(script.function)
Example: help(leap_remove)
Note: If the user wants multi year values comparison one must write a simple script to append the time series of each year before using ploting function.
[ ]:
from post_processing_stopet import *
help(leap_remove)
Example: let us remove the leap year data from 1996 timeseries first read the 1996 PET data and check the array length
Kenya 1996 data doesn’t exist
[ ]:
data = np.genfromtxt('results/Kenya_E0_StoPET/1996_3.8_36.6_3_stoPET.txt')
print(len(data))
new_data = leap_remove(data)
print(len(new_data))
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.
Example
Now lets do a simple exercise based on a 10 year data for a single location in Eastern Kenya (Wajir).
* lat = 1.73
* lon = 40.09
* start year = 2000
* end year = 2010
* two ensembles
* using Method 3 for temperature adjustment
* temperature increase of 1.5 degrees
Adjust the run_stoPET.py acordingly and run to generate the PET data
[ ]:
def run_stoPET():
## ----- CHANGE THE INPUT VARIABLES HERE -----##
datapath = 'stopet_parameter_files/'
outputpath = 'results/'
runtype = 'single' #'regional' #
startyear = 2000
endyear = 2010
# Single point stoPET run
latval = 1.73
lonval = 40.09
# Regional stoPET run
latval_min = -5.5
latval_max = -4.5 #5.5
lonval_min = 33.0
lonval_max = 34.5 #42.0
locname = 'Wajir' #'Turkana1' #
number_ensm = 2
tempAdj = 3
deltat = 1.5
udpi_pet = 5
## ------ NO CHANGES BELLOW THIS -------------##
if runtype == 'single':
for ens_num in np.arange(0,number_ensm):
stoPET_wrapper_singlepoint(startyear, endyear,
latval, lonval, locname,
ens_num,datapath, outputpath,
tempAdj, deltat, udpi_pet)
elif runtype == 'regional':
for ens_num in np.arange(0,number_ensm):
stoPET_wrapper_regional(startyear, endyear,
latval_min, latval_max,
lonval_min, lonval_max,
locname, ens_num, datapath,
outputpath,
tempAdj, deltat, udpi_pet)
else:
raise ValueError('runtype only takes single and regional ... please check!')
[ ]:
run_stoPET()
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
[ ]:
# read the PET data for the year without temperature adjustment
data_1 = np.genfromtxt('results/Wajir_E0_StoPET/2000_1.73_40.09_3_stoPET.txt')
# read the PET data for the year with adjusted temeperature
data_2 = np.genfromtxt('results/Wajir_E0_StoPET/2009_1.73_40.09_3_AdjstoPET.txt')
# lets remove the leap year day data
data_1 = leap_remove(data_1)
data_2 = leap_remove(data_2)
# check the length of the array
print(len(data_1))
print(len(data_2))
[ ]:
# Density plot
label_1 = 'stoPET'
label_2 = 'Adj. stoPET'
xlabel = 'PET ($\mathbf{mm\,hour^{-1}}$)'
ylabel = 'Density'
title = 'Hourly PET value'
plotpath = './plots/'
fig=plt.figure()
# Draw the density plot
sns.kdeplot(data_1, color = 'k',label = label_1)
# Draw the density plot
sns.kdeplot(data_2, color = 'orange',label = label_2)
plt.xlabel(xlabel)
plt.ylabel(ylabel)
plt.title(title,fontweight='bold')
plt.legend(loc='best')
plt.tight_layout()
[ ]:
# Now lets plot the timeseres for both PET
# values to visulize the daily timeseris of
# non-adjusted and adjusted PET
label_1 = 'stoPET'
label_2 = 'Adj. stoPET'
xlabel = 'PET ($\mathbf{mm\,hour^{-1}}$)'
ylabel = 'Density'
title = 'Hourly PET value'
plotpath = './plots/'
fname = 'Wajir_2000_hourly_density_stoPET.png'
# by using the function for plotting from the
# post_processing_stopet.py the figures will be
# saved in the folder user provided.
comparison_density_plot(data_1, data_2, label_1, label_2,
xlabel, ylabel, title, plotpath, fname)
Now lets make the daily and monthly aggregate for the two datasets using the aggregate_data(timeseries, period) function.
[ ]:
# daily
data_1_day = aggregate_data(data_1, 'day')
data_2_day = aggregate_data(data_2, 'day')
print(len(data_1_day))
print(len(data_2_day))
# check for the array length
# monthly
data_1_month = aggregate_data(data_1, 'month')
data_2_month = aggregate_data(data_2, 'month')
# check for the array length
print(len(data_1_month))
print(len(data_2_month))
Now lets plot the timeseres for both PET values to visulize the daily timeseris of non-adjusted and adjusted PET.
[ ]:
data = data_1_day
xlabel = 'Day of year'
ylabel = 'PET ($\mathbf{mm\,day^{-1}}$)'
title = 'Daily PET value'
plotpath = './plots/'
fname = 'Wajir_2000_daily_stoPET.png'
timeseries_plot(data, xlabel, ylabel, title, plotpath, fname)
[ ]:
# This is just to show the plots in the document the function
# already saved the plot in the plots folder.
fig=plt.figure()
plt.plot(data,'k')
plt.ylabel(ylabel)
plt.xlabel(xlabel)
plt.title(title,fontweight='bold')
plt.tight_layout()
[ ]:
# Now lets plot the timeseres for the monthly PET values of the year 2000.
data = data_1_month
xlabel = 'Month'
ylabel = 'PET ($\mathbf{mm\,month^{-1}}$)'
title = 'Monthly PET value'
plotpath = 'plots/'
fname = 'Wajir_2000_monthly_stoPET.png'
timeseries_plot(data, xlabel, ylabel, title, plotpath, fname)
[ ]:
# This is just to show the plots in the document the function
# already saved the plot in the plots folder.
fig=plt.figure()
plt.plot(data,'k')
plt.ylabel(ylabel)
plt.xlabel(xlabel)
plt.title(title,fontweight='bold')
plt.tight_layout()
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.
[ ]:
# first lets make the year list
years = np.arange(2000,2011)
# create empty array
data_1_year = [] # this is non adjusted PET
data_2_year = [] # this is the adjusted PET
# make a loop and estimate the annual PET value for both data
for i in range (0, len(years)):
year = years[i]
data_1 = np.genfromtxt('results/Wajir_E0_StoPET/%s_1.73_40.09_3_stoPET.txt'
%year)
data_2 = np.genfromtxt('results/Wajir_E0_StoPET/%s_1.73_40.09_3_AdjstoPET.txt'
%year)
# make the annual aggregate
data_1_val = aggregate_data(data_1, 'year')
data_2_val = aggregate_data(data_2, 'year')
# append the data to the empty array
data_1_year = np.append(data_1_year, data_1_val)
data_2_year = np.append(data_2_year, data_2_val)
print(year)
[ ]:
# now let us plot the two datasets for comparison
data_1 = data_1_year
data_2 = data_2_year
label_1 = 'stoPET'
label_2 = 'Adj. stoPET'
xlabel = 'Year'
ylabel = 'PET ($\mathbf{mm\,year^{-1}}$)'
title = 'Annual PET value'
plotpath = 'plots/'
fname = 'Wajir_annual_stoPET_comparison.png'
comparison_timeseries_plot(data_1, data_2, label_1, label_2,
xlabel, ylabel, title, plotpath, fname)
[ ]:
# This is just to show the plots in the document the function
# already saved the plot in the plots folder.
fig=plt.figure()
plt.plot(data_1,'k', label=label_1)
plt.plot(data_2,'k--', label=label_2)
plt.ylabel(ylabel)
plt.xlabel(xlabel)
plt.title(title,fontweight='bold')
plt.legend(loc='best')
plt.tight_layout()
[ ]: