This notebook is meant to demonstrate the consistency of the ibicus implementation of ISIMIP3BASD v3.0.1 (most recent version, in the following referred to as ISIMIP) with the reference implementation provided by Lange 2022:
The code in the ibicus implementation is based upon Lange 2019 and the ISIMIP3b bias adjustment fact sheet. More details of the implementation of ISIMIP within ibicus can be found in the documentation.
Numerous checks were conducted during and after the development process to ensure consistency. This notebook provides a final validation and demonstration of this consistency on the testing data published with the ISIMIP reference implementation.
The structure of this notebook is as follows. First the ISIMIP reference implementation is pulled from https://zenodo.org/record/6758997/files/isimip3basd-master.tar.gz and run on the testing data. In a second step debiasers are initialised using ibicus and run for all variables included in ISIMIP. Consistency is demonstrated by plotting and comparing the values computed by the reference implementation and the ibicus implementation.
2.1. Helpers
2.2. tas
2.3. pr
2.4. ps / psl
2.5. rlds
2.6. sfcWind
2.7. tasrange
2.8. tasskew
2.9. hurs
2.10. rsds
2.11. prsnratio
References:
import iris
import numpy as np
import scipy.stats
from cf_units import num2date
import matplotlib.pyplot as plt
We pull v3.0.1 of the ISIMIP reference implementation -- the latest version and the version implemented in ibicus and unpack the code:
!wget https://zenodo.org/record/6758997/files/isimip3basd-master.tar.gz -c
!tar -xf isimip3basd-master.tar.gz
--2022-08-30 12:25:16-- https://zenodo.org/record/6758997/files/isimip3basd-master.tar.gz Resolving zenodo.org (zenodo.org)... 137.138.76.77 Connecting to zenodo.org (zenodo.org)|137.138.76.77|:443... connected. HTTP request sent, awaiting response... 206 Partial Content Length: 12112564 (12M), 3167934 (3,0M) remaining [application/octet-stream] Saving to: ‘isimip3basd-master.tar.gz’ isimip3basd-master. 100%[++++++++++++++=====>] 11,55M 532KB/s in 5,8s 2022-08-30 12:25:25 (532 KB/s) - ‘isimip3basd-master.tar.gz’ saved [12112564/12112564]
We can now let the ISIMIP reference code run through.
The command below is the standard call for ISIMIP debiasing provided application_example.sh
inside the unpacked isimip3basd-master.tar.gz
. Three small modifications were made:
bias_adjustment.py
and to the data (in -o
, -s
, -f
and -b
) were adjusted.--n-quantiles 15000
was added to ensure a high number of quantiles for comparison: this is slower but more exact. With a lower number of quantiles differences between the ibicus implementation and the ISIMIP reference get bigger, due to the ISIMIP linear interpolation being not exact. Those are not directly used, but the number gets reduced to a high isimip default.2> isimip_output.txt
was added to write stderr
to file.!python -u isimip3basd-master/code/bias_adjustment.py \
--n-quantiles 15000 \
--n-processes 5 \
--randomization-seed 0 \
--step-size 1 \
-v hurs,pr,prsnratio,ps,rlds,rsds,sfcWind,tas,tasrange,tasskew \
--lower-bound 0,0,0,,,0,0,,0,0 \
--lower-threshold .01,.0000011574,.0001,,,.0001,.01,,.01,.0001 \
--upper-bound 100,,1,,,1,,,,1 \
--upper-threshold 99.99,,.9999,,,.9999,,,,.9999 \
--distribution ,gamma,,normal,normal,,weibull,normal,weibull, \
-t bounded,mixed,bounded,additive,additive,bounded,mixed,additive,mixed,bounded \
--unconditional-ccs-transfer 1,,,,,,,,, \
--trendless-bound-frequency 1,,,,,,,,, \
-d ,,,1,1,,,1,, \
-w 0,0,0,0,0,15,0,0,0,0 \
--if-all-invalid-use ,,0.,,,,,,, \
-o isimip3basd-master/data/hurs_obs-hist_coarse_1979-2014.nc,isimip3basd-master/data/pr_obs-hist_coarse_1979-2014.nc,isimip3basd-master/data/prsnratio_obs-hist_coarse_1979-2014.nc,isimip3basd-master/data/ps_obs-hist_coarse_1979-2014.nc,isimip3basd-master/data/rlds_obs-hist_coarse_1979-2014.nc,isimip3basd-master/data/rsds_obs-hist_coarse_1979-2014.nc,isimip3basd-master/data/sfcWind_obs-hist_coarse_1979-2014.nc,isimip3basd-master/data/tas_obs-hist_coarse_1979-2014.nc,isimip3basd-master/data/tasrange_obs-hist_coarse_1979-2014.nc,isimip3basd-master/data/tasskew_obs-hist_coarse_1979-2014.nc \
-s isimip3basd-master/data/hurs_sim-hist_coarse_1979-2014.nc,isimip3basd-master/data/pr_sim-hist_coarse_1979-2014.nc,isimip3basd-master/data/prsnratio_sim-hist_coarse_1979-2014.nc,isimip3basd-master/data/ps_sim-hist_coarse_1979-2014.nc,isimip3basd-master/data/rlds_sim-hist_coarse_1979-2014.nc,isimip3basd-master/data/rsds_sim-hist_coarse_1979-2014.nc,isimip3basd-master/data/sfcWind_sim-hist_coarse_1979-2014.nc,isimip3basd-master/data/tas_sim-hist_coarse_1979-2014.nc,isimip3basd-master/data/tasrange_sim-hist_coarse_1979-2014.nc,isimip3basd-master/data/tasskew_sim-hist_coarse_1979-2014.nc \
-f isimip3basd-master/data/hurs_sim-fut_coarse_2065-2100.nc,isimip3basd-master/data/pr_sim-fut_coarse_2065-2100.nc,isimip3basd-master/data/prsnratio_sim-fut_coarse_2065-2100.nc,isimip3basd-master/data/ps_sim-fut_coarse_2065-2100.nc,isimip3basd-master/data/rlds_sim-fut_coarse_2065-2100.nc,isimip3basd-master/data/rsds_sim-fut_coarse_2065-2100.nc,isimip3basd-master/data/sfcWind_sim-fut_coarse_2065-2100.nc,isimip3basd-master/data/tas_sim-fut_coarse_2065-2100.nc,isimip3basd-master/data/tasrange_sim-fut_coarse_2065-2100.nc,isimip3basd-master/data/tasskew_sim-fut_coarse_2065-2100.nc \
-b isimip3basd-master/data/hurs_sim-fut-basd_coarse_2065-2100.nc,isimip3basd-master/data/pr_sim-fut-basd_coarse_2065-2100.nc,isimip3basd-master/data/prsnratio_sim-fut-basd_coarse_2065-2100.nc,isimip3basd-master/data/ps_sim-fut-basd_coarse_2065-2100.nc,isimip3basd-master/data/rlds_sim-fut-basd_coarse_2065-2100.nc,isimip3basd-master/data/rsds_sim-fut-basd_coarse_2065-2100.nc,isimip3basd-master/data/sfcWind_sim-fut-basd_coarse_2065-2100.nc,isimip3basd-master/data/tas_sim-fut-basd_coarse_2065-2100.nc,isimip3basd-master/data/tasrange_sim-fut-basd_coarse_2065-2100.nc,isimip3basd-master/data/tasskew_sim-fut-basd_coarse_2065-2100.nc \
2> isimip_output.txt
checking inputs ... adjusting at location (lon, lat) ... (0, 0) (0, 1) (1, 0) (1, 1)
A number of warnings are generated and written to isimip_output.txt
indicating that the number of quantiles get reduced to a high ISIMIP default. This is not a problem for the values.
After having executed the reference implementation we can now come to the ibicus implementation. Let' s import it:
from ibicus.debias import ISIMIP
We first define some helpers to read in the ISIMIP-testing and debiased data and gets the dates (last coordinate). Those are needed for the debiasing below.
# Given an iris-cube this returns the dates stored in the last time-dimension
def get_dates(x):
time_dimension = x.coords()[2]
dates = time_dimension.units.num2date(time_dimension.points)
return dates
get_dates = np.vectorize(get_dates)
# This reads in the testing-data from ISIMIP stored in isimip3basd-master/data
def read_in_and_preprocess_isimip_testing_data_with_dates(variable, isimip_data_path = "isimip3basd-master/data/"):
# Load in data
obs = iris.load_cube(isimip_data_path+variable+"_obs-hist_coarse_1979-2014.nc")
cm_hist = iris.load_cube(isimip_data_path+variable+"_sim-hist_coarse_1979-2014.nc")
cm_future = iris.load_cube(isimip_data_path+variable+"_sim-fut_coarse_2065-2100.nc")
# Extract dates
dates = {
"time_obs": get_dates(obs),
"time_cm_hist": get_dates(cm_hist),
"time_cm_future": get_dates(cm_future)
}
# Convert to np.array (from masked-array)
obs = obs.data
cm_hist = cm_hist.data
cm_future = cm_future.data
# Move time to first axis (our convention)
obs = np.moveaxis(obs, -1, 0)
cm_hist = np.moveaxis(cm_hist, -1, 0)
cm_future = np.moveaxis(cm_future, -1, 0)
return obs, cm_hist, cm_future, dates
def read_in_debiased_testing_data(variable, isimip_data_path = "isimip3basd-master/data/"):
# Load in data
debiased_data = iris.load_cube(isimip_data_path+variable+"_sim-fut-basd_coarse_2065-2100.nc")
# Move time to first axis (our convention)
debiased_data = np.array(debiased_data.data)
debiased_data = np.moveaxis(debiased_data, -1, 0)
return debiased_data
Run the ibicus ISIMIP implementation:
variable = "tas"
obs, cm_hist, cm_future, dates = read_in_and_preprocess_isimip_testing_data_with_dates(variable)
debiaser = ISIMIP.from_variable(variable)
debiased_values = debiaser.apply(obs, cm_hist, cm_future, **dates)
INFO:root:----- Running debiasing for variable: Daily mean near-surface air temperature ----- INFO:root:obs is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_hist is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_future is a masked array, but contains no invalid data. It is converted to a normal numpy array. 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:20<00:00, 5.00s/it]
Compare:
debiased_values_isimip = read_in_debiased_testing_data(variable)
pct_agreement = np.sum(np.isclose(debiased_values,debiased_values_isimip))/debiased_values.size
print("Percentage agreement between the ibicus ISIMIP and the reference implementation is %s %%."%(pct_agreement*100))
Percentage agreement between the Ibicus ISIMIP and the reference implementation is 100.0 %.
We can plot the ibicus and reference implementation over time at one location:
min_time = 0
max_time = 100
lat = 1
lon = 1
time = np.arange(min_time, max_time)
plt.plot(time, debiased_values_isimip[min_time:max_time, lat, lon], label = "ISIMIP reference implementation")
plt.plot(time, debiased_values[min_time:max_time, lat, lon], label = "ibicus implementation", linestyle = ":")
plt.plot(time, cm_future[min_time:max_time, lat, lon], label = "Raw CM")
plt.legend()
plt.show()
We see that the dotted ibicus implementation line lies right above the ISIMIP implementation line.
We can also plot the values against eachother:
min_all_values = np.min([np.min(debiased_values), np.min(debiased_values_isimip)])
max_all_values = np.max([np.max(debiased_values), np.max(debiased_values_isimip)])
plt.scatter(debiased_values_isimip.flatten(), debiased_values.flatten())
plt.plot([min_all_values, max_all_values], [min_all_values, max_all_values])
plt.show()
A linear regression also shows that the reference implementation values' and the ibicus ones are consistent:
scipy.stats.linregress(debiased_values_isimip.flatten(), debiased_values.flatten())
LinregressResult(slope=0.9999998728613527, intercept=3.637642799958485e-05, rvalue=0.9999999999943137, pvalue=0.0, stderr=1.4705499476000933e-08, intercept_stderr=4.208917712180163e-06)
We conclude that tas is reproduced well by the ibicus implementation of ISIMIP.
Run the ibicus ISIMIP implementation:
variable = "pr"
obs, cm_hist, cm_future, dates = read_in_and_preprocess_isimip_testing_data_with_dates(variable)
debiaser = ISIMIP.from_variable(variable)
debiased_values = debiaser.apply(obs, cm_hist, cm_future, **dates)
INFO:root:----- Running debiasing for variable: Daily mean precipitation ----- INFO:root:obs is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_hist is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_future is a masked array, but contains no invalid data. It is converted to a normal numpy array. 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:16<00:00, 4.14s/it]
Compare:
debiased_values_isimip = read_in_debiased_testing_data(variable)
pct_agreement = np.sum(np.abs(debiased_values-debiased_values_isimip) < 1e-6)/debiased_values.size
print("%s %% of all values of the ibicus implementation are within 1e-6 distance of the reference implementation."%(pct_agreement*100))
99.96197140249467 % of all values of the Ibicus implementation are within 1e-6 distance of the reference implementation.
Plot over time:
min_time = 0
max_time = 100
lat = 1
lon = 1
time = np.arange(min_time, max_time)
plt.plot(time, debiased_values_isimip[min_time:max_time, lat, lon], label = "ISIMIP reference implementation")
plt.plot(time, debiased_values[min_time:max_time, lat, lon], label = "ibicus implementation", linestyle = ":")
plt.plot(time, cm_future[min_time:max_time, lat, lon], label = "Raw CM")
plt.legend()
plt.show()
Values of reference implementation against the ibicus one:
min_all_values = np.min([np.min(debiased_values), np.min(debiased_values_isimip)])
max_all_values = np.max([np.max(debiased_values), np.max(debiased_values_isimip)])
plt.scatter(debiased_values_isimip.flatten(), debiased_values.flatten())
plt.plot([min_all_values, max_all_values], [min_all_values, max_all_values])
plt.show()
It seems that pr is reproduced quite well and the values of the reference and ibicus implementation are in agreement. This can also be checked through a linear regression:
scipy.stats.linregress(debiased_values_isimip.flatten(), debiased_values.flatten())
LinregressResult(slope=1.0000577823976702, intercept=8.037901231378822e-10, rvalue=0.9999993853973769, pvalue=0.0, stderr=4.834882356821364e-06, intercept_stderr=2.2981006288726986e-10)
pr is reproduced well by ibicus. Some slight differences larger than floating point error exist. This is due to:
Randomization: pr includes some randomization between the defined bound and threshold. This can lead to differences.
The references implementation of nonparametric quantile mapping (preceding the parametric one in step 6), which uses linear interpolation, is inexact and differs from the ibicus implementation of nonparametric quantile mapping. This creates some differences. These decrease with the number of quantiles increasing, however they are slightly bigger than floating point error.
Accumulation of floating point errors in calculations. Especially floating point errors in the computation of quantiles can lead to slight numerical differences (larger than floating point) if those quantiles are mapped back to values. Similarly the distribution fits in step 6 are just slightly different (within floating point accuracy), meaning that the same values are mapped to slighty different ones (with difference potentially larger than floating point error) when transformed using an (inverse) CDF.
Run the ibicus ISIMIP implementation:
variable = "ps"
obs, cm_hist, cm_future, dates = read_in_and_preprocess_isimip_testing_data_with_dates(variable)
debiaser = ISIMIP.from_variable(variable)
debiased_values = debiaser.apply(obs, cm_hist, cm_future, **dates)
INFO:root:----- Running debiasing for variable: Daily mean sea-level pressure ----- INFO:root:obs is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_hist is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_future is a masked array, but contains no invalid data. It is converted to a normal numpy array. 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:21<00:00, 5.47s/it]
Compare:
debiased_values_isimip = read_in_debiased_testing_data("ps")
pct_agreement = np.sum(np.isclose(debiased_values,debiased_values_isimip))/debiased_values.size
print("Percentage agreement is %s %% between the ibicus and the reference implementation of ISIMIP."%(pct_agreement*100))
Percentage agreement is 100.0 % between the Ibicus and the reference implementation of ISIMIP.
Plot over time at location [1,1]:
min_time = 0
max_time = 100
lat = 1
lon = 1
time = np.arange(min_time, max_time)
plt.plot(time, debiased_values_isimip[min_time:max_time, lat, lon], label = "ISIMIP")
plt.plot(time, debiased_values[min_time:max_time, lat, lon], label = "ibicus implementation", linestyle = ":")
plt.plot(time, cm_future[min_time:max_time, lat, lon], label = "raw cm")
plt.legend()
plt.show()
Values of reference implementation against the ibicus one:
min_all_values = np.min([np.min(debiased_values), np.min(debiased_values_isimip)])
max_all_values = np.max([np.max(debiased_values), np.max(debiased_values_isimip)])
plt.scatter(debiased_values_isimip.flatten(), debiased_values.flatten())
plt.plot([min_all_values, max_all_values], [min_all_values, max_all_values])
plt.show()
A linear regression also shows that the reference implementation values' and the ibicus ones are consistent:
scipy.stats.linregress(debiased_values_isimip.flatten(), debiased_values.flatten())
LinregressResult(slope=0.9999999951431144, intercept=0.0004899313935311511, rvalue=0.9999999999776888, pvalue=0.0, stderr=2.9128922281955855e-08, intercept_stderr=0.0029384674434487085)
ps / psl is well reproduced by ibicus.
Run the ibicus ISIMIP implementation:
variable = "rlds"
obs, cm_hist, cm_future, dates = read_in_and_preprocess_isimip_testing_data_with_dates(variable)
debiaser = ISIMIP.from_variable(variable)
debiased_values = debiaser.apply(obs, cm_hist, cm_future, **dates)
INFO:root:----- Running debiasing for variable: Daily mean surface downwelling longwave radiation ----- INFO:root:obs is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_hist is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_future is a masked array, but contains no invalid data. It is converted to a normal numpy array. 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:20<00:00, 5.03s/it]
Compare:
debiased_values_isimip = read_in_debiased_testing_data(variable)
pct_agreement = np.sum(np.isclose(debiased_values,debiased_values_isimip))/debiased_values.size
print("Percentage agreement is %s %% between the ibicus and the reference implementation of ISIMIP."%(pct_agreement*100))
Percentage agreement is 100.0 % between the Ibicus and the reference implementation of ISIMIP.
We can plot the values over time at location [1,1]:
min_time = 0
max_time = 100
lat = 1
lon = 1
time = np.arange(min_time, max_time)
plt.plot(time, debiased_values_isimip[min_time:max_time, lat, lon], label = "ISIMIP")
plt.plot(time, debiased_values[min_time:max_time, lat, lon], label = "ibicus implementation", linestyle = ":")
plt.plot(time, cm_future[min_time:max_time, lat, lon], label = "raw cm")
plt.legend()
plt.show()
Values of reference implementation against the ibicus one:
min_all_values = np.min([np.min(debiased_values), np.min(debiased_values_isimip)])
max_all_values = np.max([np.max(debiased_values), np.max(debiased_values_isimip)])
plt.scatter(debiased_values_isimip.flatten(), debiased_values.flatten())
plt.plot([min_all_values, max_all_values], [min_all_values, max_all_values])
plt.show()
A linear regression also shows that the reference implementation values' and the ibicus ones are consistent:
scipy.stats.linregress(debiased_values_isimip.flatten(), debiased_values.flatten())
LinregressResult(slope=1.0000000224862213, intercept=-6.843278475798797e-05, rvalue=0.9999999999997392, pvalue=0.0, stderr=3.14927108323425e-09, intercept_stderr=1.043150503623679e-06)
rlds is well reproduced by ibicus.
Run the ibicus ISIMIP implementation:
variable = "sfcWind"
obs, cm_hist, cm_future, dates = read_in_and_preprocess_isimip_testing_data_with_dates(variable)
debiaser = ISIMIP.from_variable(variable)
debiased_values = debiaser.apply(obs, cm_hist, cm_future, **dates)
INFO:root:----- Running debiasing for variable: Daily mean near-surface wind speed ----- INFO:root:obs is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_hist is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_future is a masked array, but contains no invalid data. It is converted to a normal numpy array. 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:59<00:00, 14.92s/it]
Compare:
debiased_values_isimip = read_in_debiased_testing_data(variable)
pct_agreement = np.sum(np.abs(debiased_values -debiased_values_isimip) < 1e-3)/debiased_values.size
max_deviation = np.max(np.abs(debiased_values-debiased_values_isimip))
print(f"{pct_agreement*100}% of all values are within 1e-3. There is a maximum deviation of {max_deviation} which is {100*max_deviation/np.mean(debiased_values_isimip)}% of the average value")
98.03392150897476% of all values are within 1e-3. There is a maximum deviation of 0.0035228729248046875 which is 0.08669363071705068% of the average value
We can plot the values over time at location [1,1]:
min_time = 0
max_time = 100
lat = 1
lon = 1
time = np.arange(min_time, max_time)
plt.plot(time, debiased_values_isimip[min_time:max_time, lat, lon], label = "ISIMIP")
plt.plot(time, debiased_values[min_time:max_time, lat, lon], label = "ibicus implementation", linestyle = ":")
plt.plot(time, cm_future[min_time:max_time, lat, lon], label = "raw cm")
plt.legend()
plt.show()
Values of reference implementation against ibicus ones:
min_all_values = np.min([np.min(debiased_values), np.min(debiased_values_isimip)])
max_all_values = np.max([np.max(debiased_values), np.max(debiased_values_isimip)])
plt.scatter(debiased_values_isimip.flatten(), debiased_values.flatten())
plt.plot([min_all_values, max_all_values], [min_all_values, max_all_values])
plt.show()
A linear regression also shows that the reference implementation values' and the ibicus ones are consistent:
scipy.stats.linregress(debiased_values_isimip.flatten(), debiased_values.flatten())
LinregressResult(slope=1.0000021330821292, intercept=-8.191132122803424e-06, rvalue=0.9999999715602359, pvalue=0.0, stderr=1.039985640948038e-06, intercept_stderr=4.5461419994753345e-06)
sfcWind is reproduced well by ibicus. Some differences larger than floating point error exist. This is due to:
Randomization: sfcWind includes some randomization between lower bound and threshold. This can lead to differences.
The references implementation of nonparametric quantile mapping (preceding the parametric one in step 6), which uses linear interpolation, is inexact and differs from the ibicus ones. This creates some differences. These decrease with the number of quantiles increasing, however they are slightly bigger than floating point error.
Accumulation of floating point errors in calculations. Especially floating point errors in the computation of quantiles can lead to slight numerical differences (larger than floating point) if those quantiles are mapped back to values. Similarly the distribution fits in step 6 are just slightly different (within floating point accuracy), meaning that the same values are mapped to slighty different ones (with difference potentially larger than floating point error) if transformed with an (inverse) CDF.
Run the ibicus ISIMIP implementation:
variable = "tasrange"
obs, cm_hist, cm_future, dates = read_in_and_preprocess_isimip_testing_data_with_dates(variable)
debiaser = ISIMIP.from_variable(variable)
debiased_values = debiaser.apply(obs, cm_hist, cm_future, **dates)
INFO:root:----- Running debiasing for variable: Daily near-surface air temperature range (tasmax-tasmin) ----- INFO:root:obs is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_hist is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_future is a masked array, but contains no invalid data. It is converted to a normal numpy array. 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4/4 [01:13<00:00, 18.33s/it]
Compare:
debiased_values_isimip = read_in_debiased_testing_data(variable)
pct_agreement = np.sum(np.abs(debiased_values -debiased_values_isimip) < 1e-3)/debiased_values.size
max_deviation = np.max(np.abs(debiased_values-debiased_values_isimip))
print(f"{pct_agreement*100}% of all values are within 1e-3. There is a maximum deviation of {max_deviation} which is {100*max_deviation/np.mean(debiased_values_isimip)}% of the average value")
90.94348950410709% of all values are within 1e-3. There is a maximum deviation of 0.0056781768798828125 which is 0.0753957246673246% of the average value
We can plot the values over time at location [1,1]:
min_time = 0
max_time = 100
lat = 1
lon = 1
time = np.arange(min_time, max_time)
plt.plot(time, debiased_values_isimip[min_time:max_time, lat, lon], label = "ISIMIP")
plt.plot(time, debiased_values[min_time:max_time, lat, lon], label = "ibicus implementation", linestyle = ":")
plt.plot(time, cm_future[min_time:max_time, lat, lon], label = "raw cm")
plt.legend()
plt.show()
Values of reference implementation against the ibicus ones:
min_all_values = np.min([np.min(debiased_values), np.min(debiased_values_isimip)])
max_all_values = np.max([np.max(debiased_values), np.max(debiased_values_isimip)])
plt.scatter(debiased_values_isimip.flatten(), debiased_values.flatten())
plt.plot([min_all_values, max_all_values], [min_all_values, max_all_values])
plt.show()
A linear regression also shows that the reference implementation values' and the ibicus ones are consistent:
scipy.stats.linregress(debiased_values_isimip.flatten(), debiased_values.flatten())
LinregressResult(slope=0.9999986970893849, intercept=1.5534482123769067e-05, rvalue=0.9999999844855442, pvalue=0.0, stderr=7.681236296246608e-07, intercept_stderr=6.3467097659119465e-06)
ibicus can reproduce tasrange well. Some differences larger than floating point error exist. This is due to:
Randomization: tasrange includes some randomization between bound and threshold. This can lead to differences.
The references implementation of nonparametric quantile mapping (preceding the parametric one in step 6), which uses linear interpolation, is inexact and differs from the ibicus ones. This creates some differences. These decrease with the number of quantiles increasing, however they are slightly bigger than floating point error.
Accumulation of floating point errors in calculations. Especially floating point errors in the computation of quantiles can lead to slight numerical differences (larger than floating point) if those quantiles are mapped back to values. Similarly the distribution fits in step 6 are just slightly different (within floating point accuracy), meaning that the same values are mapped to slighty different ones (with difference potentially larger than floating point error) if transformed with an (inverse) CDF.
Run the ibicus ISIMIP implementation:
variable = "tasskew"
obs, cm_hist, cm_future, dates = read_in_and_preprocess_isimip_testing_data_with_dates(variable)
debiaser = ISIMIP.from_variable(variable)
debiased_values = debiaser.apply(obs, cm_hist, cm_future, **dates)
INFO:root:----- Running debiasing for variable: Daily near-surface air temperature skew (tas-tasmin)/tasrange ----- INFO:root:obs is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_hist is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_future is a masked array, but contains no invalid data. It is converted to a normal numpy array. 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:20<00:00, 5.15s/it]
Compare:
debiased_values_isimip = read_in_debiased_testing_data(variable)
pct_agreement = np.sum(np.abs(debiased_values -debiased_values_isimip) < 1e-3)/debiased_values.size
max_deviation = np.max(np.abs(debiased_values-debiased_values_isimip))
print(f"{pct_agreement*100}% of all values are within 1e-3. There is a maximum deviation of {max_deviation} which is {100*max_deviation/np.mean(debiased_values_isimip)}% of the average value")
99.9942957103742% of all values are within 1e-3. There is a maximum deviation of 0.0015919804573059082 which is 0.3169975187576745% of the average value
We can plot the values over time at location [1,1]:
min_time = 0
max_time = 100
lat = 1
lon = 1
time = np.arange(min_time, max_time)
plt.plot(time, debiased_values_isimip[min_time:max_time, lat, lon], label = "ISIMIP")
plt.plot(time, debiased_values[min_time:max_time, lat, lon], label = "ibicus implementation", linestyle = ":")
plt.plot(time, cm_future[min_time:max_time, lat, lon], label = "raw cm")
plt.legend()
plt.show()
Values of reference implementation against ibicus ones:
min_all_values = np.min([np.min(debiased_values), np.min(debiased_values_isimip)])
max_all_values = np.max([np.max(debiased_values), np.max(debiased_values_isimip)])
plt.scatter(debiased_values_isimip.flatten(), debiased_values.flatten())
plt.plot([min_all_values, max_all_values], [min_all_values, max_all_values])
plt.show()
Linear Regression:
scipy.stats.linregress(debiased_values_isimip.flatten(), debiased_values.flatten())
LinregressResult(slope=0.9999931206025992, intercept=3.6336887746513113e-06, rvalue=0.9999999241443093, pvalue=0.0, stderr=1.6984569092665668e-06, intercept_stderr=8.60771269891601e-07)
tasskew is reproduced well by ibicus. Some differences larger than floating point error exist. This is due to:
Randomization: tasrange includes some randomization between both upper and lower bound and threshold. This can lead to differences.
The references implementation of nonparametric quantile mapping using linear interpolation is inexact and differs from the ibicus ones. This creates some differences. These decrease with the number of quantiles increasing, however they are slightly bigger than floating point error.
Accumulation of floating point errors in calculations. Especially floating point errors in the computation of quantiles can lead to slight numerical differences (larger than floating point) if those quantiles are mapped back to values.
Run the ibicus ISIMIP implementation:
variable = "hurs"
obs, cm_hist, cm_future, dates = read_in_and_preprocess_isimip_testing_data_with_dates(variable)
debiaser = ISIMIP.from_variable(variable)
debiased_values = debiaser.apply(obs, cm_hist, cm_future, **dates)
INFO:root:----- Running debiasing for variable: Daily mean near-surface relative humidity ----- INFO:root:obs is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_hist is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_future is a masked array, but contains no invalid data. It is converted to a normal numpy array. 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:19<00:00, 4.99s/it]
Compare:
debiased_values_isimip = read_in_debiased_testing_data(variable)
pct_agreement = np.sum(np.abs(debiased_values -debiased_values_isimip) < 1e-2)/debiased_values.size
max_deviation = np.max(np.abs(debiased_values-debiased_values_isimip))
print(f"{pct_agreement*100}% of all values are within 1e-2 with a maximum deviation of {max_deviation} which is {np.round(100*max_deviation/np.mean(debiased_values_isimip), 4)}% of the average value")
99.57598113781565% of all values are within 1e-2 with a maximum deviation of 0.29793548583984375 which is 0.3875% of the average value
We can plot the values over time at location [1,1]:
min_time = 0
max_time = 100
lat = 1
lon = 1
time = np.arange(min_time, max_time)
plt.plot(time, debiased_values_isimip[min_time:max_time, lat, lon], label = "ISIMIP")
plt.plot(time, debiased_values[min_time:max_time, lat, lon], label = "ibicus implementation", linestyle = ":")
plt.plot(time, cm_future[min_time:max_time, lat, lon], label = "raw cm")
plt.legend()
plt.show()
Values of reference implementation against ibicus ones:
min_all_values = np.min([np.min(debiased_values), np.min(debiased_values_isimip)])
max_all_values = np.max([np.max(debiased_values), np.max(debiased_values_isimip)])
plt.scatter(debiased_values_isimip.flatten(), debiased_values.flatten())
plt.plot([min_all_values, max_all_values], [min_all_values, max_all_values])
plt.show()
Linear Regression:
scipy.stats.linregress(debiased_values_isimip.flatten(), debiased_values.flatten())
LinregressResult(slope=0.9999991927516926, intercept=6.969084543584358e-05, rvalue=0.9999999728811985, pvalue=0.0, stderr=1.0155429998726147e-06, intercept_stderr=7.894066571529341e-05)
hurs is reproduced well by ibicus. Some differences larger than floating point error exist. This is due to:
Randomization: as other variables hurs includes some randomization between both upper and lower bound and threshold. This can lead to differences.
The references implementation of nonparametric quantile mapping, which uses linear interpolation, is inexact and differs from the ibicus ones. This creates some differences. These decrease with the number of quantiles increasing, however they are slightly bigger than floating point error.
Accumulation of floating point errors in calculations. Especially floating point errors in the computation of quantiles can lead to slight numerical differences (larger than floating point) if those quantiles are mapped back to values.
Run the ibicus ISIMIP implementation:
variable = "rsds"
obs, cm_hist, cm_future, dates = read_in_and_preprocess_isimip_testing_data_with_dates(variable)
debiaser = ISIMIP.from_variable(variable)
debiased_values = debiaser.apply(obs, cm_hist, cm_future, **dates)
INFO:root:----- Running debiasing for variable: Daily mean surface downwelling shortwave radiation ----- INFO:root:obs is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_hist is a masked array, but contains no invalid data. It is converted to a normal numpy array. INFO:root:cm_future is a masked array, but contains no invalid data. It is converted to a normal numpy array. 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:19<00:00, 4.97s/it]
Compare:
debiased_values_isimip = read_in_debiased_testing_data(variable)
pct_agreement = np.sum(np.abs(debiased_values -debiased_values_isimip) < 0.1)/debiased_values.size
max_deviation = np.max(np.abs(debiased_values-debiased_values_isimip))
print(f"{pct_agreement*100}% of all values are within 0.1 with a maximum deviation of {max_deviation} which is {100*max_deviation/np.mean(debiased_values_isimip)}% of the average value")
99.77753270459385% of all values are within 0.1 with a maximum deviation of 1.1534500122070312 which is 0.9207561400338631% of the average value
We can plot the values over time at location [1,1]:
min_time = 0
max_time = 100
lat = 1
lon = 1
time = np.arange(min_time, max_time)
plt.plot(time, debiased_values_isimip[min_time:max_time, lat, lon], label = "ISIMIP")
plt.plot(time, debiased_values[min_time:max_time, lat, lon], label = "ibicus implementation", linestyle = ":")
plt.plot(time, cm_future[min_time:max_time, lat, lon], label = "raw cm")
plt.legend()
plt.show()
Values of reference implementation against ibicus ones:
min_all_values = np.min([np.min(debiased_values), np.min(debiased_values_isimip)])
max_all_values = np.max([np.max(debiased_values), np.max(debiased_values_isimip)])
plt.scatter(debiased_values_isimip.flatten(), debiased_values.flatten())
plt.plot([min_all_values, max_all_values], [min_all_values, max_all_values])
plt.show()
Linear Regression:
scipy.stats.linregress(debiased_values_isimip.flatten(), debiased_values.flatten())
LinregressResult(slope=1.0000023217484482, intercept=-8.485653290790651e-05, rvalue=0.999999989836804, pvalue=0.0, stderr=6.216981478242796e-07, intercept_stderr=9.553003222892304e-05)
ibicus reproduces rsds well. Some differences larger than floating point error exist. This is due to:
Randomization: rsds includes some randomization between both upper and lower bound and threshold. This can lead to differences.
The references implementation of nonparametric quantile mapping, which uses linear interpolation, is inexact and differs from the ibicus ones. This creates some differences. These decrease with the number of quantiles increasing, however they are slightly bigger than floating point error.
Accumulation of floating point errors in calculations. Especially floating point errors in the computation of quantiles can lead to slight numerical differences (larger than floating point) if those quantiles are mapped back to values.
Run the ibicus ISIMIP implementation:
variable = "prsnratio"
obs, cm_hist, cm_future, dates = read_in_and_preprocess_isimip_testing_data_with_dates(variable)
debiaser = ISIMIP.from_variable(variable)
debiased_values = debiaser.apply(obs, cm_hist, cm_future, **dates)
INFO:root:----- Running debiasing for variable: Daily mean snowfall flux / Daily mean precipitation ----- WARNING:root:obs is a masked array and contains cells with invalid data. Not all debiasers support invalid/missing values and their presence might lead to infs or nans inside the debiased values. Consider infilling them. For computation the masked values here are filled in by nan-values. WARNING:root:cm_hist is a masked array and contains cells with invalid data. Not all debiasers support invalid/missing values and their presence might lead to infs or nans inside the debiased values. Consider infilling them. For computation the masked values here are filled in by nan-values. WARNING:root:cm_future is a masked array and contains cells with invalid data. Not all debiasers support invalid/missing values and their presence might lead to infs or nans inside the debiased values. Consider infilling them. For computation the masked values here are filled in by nan-values. 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:08<00:00, 2.07s/it]
Compare:
debiased_values_isimip = read_in_debiased_testing_data(variable)
pct_agreement = np.sum(np.isclose(debiased_values,debiased_values_isimip))/debiased_values.size
max_deviation = np.max(np.abs(debiased_values-debiased_values_isimip))
print("Percentage agreement is %s with max deviation %s"%(pct_agreement*100, max_deviation))
Percentage agreement is 92.59963492546395 with max deviation 1.0
We can plot the values over time at location [1,1]:
min_time = 0
max_time = 100
lat = 1
lon = 1
time = np.arange(min_time, max_time)
plt.plot(time, debiased_values_isimip[min_time:max_time, lat, lon], label = "ISIMIP")
plt.plot(time, debiased_values[min_time:max_time, lat, lon], label = "ibicus implementation", linestyle = ":")
plt.plot(time, cm_future[min_time:max_time, lat, lon], label = "raw cm")
plt.legend()
plt.show()
Values of reference implementation against ibicus ones:
min_all_values = np.min([np.min(debiased_values), np.min(debiased_values_isimip)])
max_all_values = np.max([np.max(debiased_values), np.max(debiased_values_isimip)])
plt.scatter(debiased_values_isimip.flatten(), debiased_values.flatten())
plt.plot([min_all_values, max_all_values], [min_all_values, max_all_values])
plt.show()
Linear Regression:
scipy.stats.linregress(debiased_values_isimip.flatten(), debiased_values.flatten())
LinregressResult(slope=0.9903892593000797, intercept=0.0008416249599787345, rvalue=0.9910176848082406, pvalue=0.0, stderr=0.0005827799956240896, intercept_stderr=0.0001465072709228235)
prsnratio includes a large amount of randomization: the metric is not defined on any day without rain (given that it is the ratio of prsn over pr) and in step 2 of the ISIMIP algorithm these values where the metric is not defined are filled with draws from all other available values in a given window (in all three: obs, cm_hist, cm_future). This randomization and random draws are not exactly reproducible using ibicus and they influence all other computed values: the transfer of the climate change signal in step 5, the computation which values are set to bounds and the nonparametric quantile mapping in step 6. However the values of ibicus and the reference implementation nevertheless seem in decent agreement and special care was taken to ensure that the outputs of all steps not including randomzation do agree.
The results of the ibicus and reference implementation of ISIMIP agree. Some smaller differences do exist -- at points slightly larger than floating point error, in particular for variables that involve randomization. However, these differences are within the small numerical differences that are to be expected in two different implementations of an algorithm involving a large variety of steps.
The results above also only serve as last "sanity check" and demonstration for the ibicus implementation of ISIMIP. The authors of the software package took special care that the results do also agree in all intermediate steps. Differences seen in the outputs here were therefore shown to be due to the differences named above (in the variables).