/
ag6_calc_yearly_cmip6_bias_trends.py
65 lines (55 loc) · 3.36 KB
/
ag6_calc_yearly_cmip6_bias_trends.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
import xarray as xr
import xesmf as xe
import numpy as np
import matplotlib.pyplot as plt
import scipy
import statsmodels.api as sm
import cartopy
import cartopy.util
import cartopy.crs as ccrs
import glob
import sys, os
import pickle, gzip
import datetime
dirProj = '/home/edcoffel/drive/MAX-Filer/Research/Climate-01/Personal-F20/edcoffel-F20/research/2020-ag-cmip6'
region = 'global'
model = sys.argv[1]
print('loading yearly bias for %s'%model)
with open('cmip6_output/bias/yearly-cmip6-era5-tasmax-max-bias-%s-%s.dat'%(region, model), 'rb') as f:
yearly_tasmax_max_bias = pickle.load(f)
with open('cmip6_output/bias/yearly-cmip6-era5-tasmax-monthly-max-bias-%s-%s.dat'%(region, model), 'rb') as f:
yearly_tasmax_monthly_max_bias = pickle.load(f)
with open('cmip6_output/bias/yearly-cmip6-era5-tasmax-mean-bias-%s-%s.dat'%(region, model), 'rb') as f:
yearly_tasmax_mean_bias = pickle.load(f)
yearly_tasmax_max_bias_trend = np.full([yearly_tasmax_max_bias.shape[1], yearly_tasmax_max_bias.shape[2]], np.nan)
yearly_tasmax_max_bias_trend_sig = np.full([yearly_tasmax_max_bias.shape[1], yearly_tasmax_max_bias.shape[2]], np.nan)
yearly_tasmax_monthly_max_bias_trend = np.full([12, yearly_tasmax_max_bias.shape[1], yearly_tasmax_max_bias.shape[2]], np.nan)
yearly_tasmax_monthly_max_bias_trend_sig = np.full([12, yearly_tasmax_max_bias.shape[1], yearly_tasmax_max_bias.shape[2]], np.nan)
yearly_tasmax_mean_bias_trend = np.full([yearly_tasmax_mean_bias.shape[1], yearly_tasmax_mean_bias.shape[2]], np.nan)
yearly_tasmax_mean_bias_trend_sig = np.full([yearly_tasmax_mean_bias.shape[1], yearly_tasmax_mean_bias.shape[2]], np.nan)
print('processing trends in cmip6 bias %s...'%model)
for xlat in range(yearly_tasmax_max_bias_trend.shape[0]):
for ylon in range(yearly_tasmax_max_bias_trend.shape[1]):
curBias = np.squeeze(yearly_tasmax_max_bias[:, xlat, ylon])
X = sm.add_constant(range(1981, 2015))
mdl = sm.RLM(curBias, X).fit()
yearly_tasmax_max_bias_trend[xlat, ylon] = mdl.params[1]*10
yearly_tasmax_max_bias_trend_sig[xlat, ylon] = mdl.pvalues[1]
curBias = np.squeeze(yearly_tasmax_mean_bias[:, xlat, ylon])
X = sm.add_constant(range(1981, 2015))
mdl = sm.RLM(curBias, X).fit()
yearly_tasmax_mean_bias_trend[xlat, ylon] = mdl.params[1]*10
yearly_tasmax_mean_bias_trend_sig[xlat, ylon] = mdl.pvalues[1]
for month in range(12):
curBias = np.squeeze(yearly_tasmax_monthly_max_bias[month, :, xlat, ylon])
X = sm.add_constant(range(1981, 2015))
mdl = sm.RLM(curBias, X).fit()
yearly_tasmax_monthly_max_bias_trend[month, xlat, ylon] = mdl.params[1]*10
yearly_tasmax_monthly_max_bias_trend_sig[month, xlat, ylon] = mdl.pvalues[1]
with open('cmip6_output/bias/yearly-cmip6-era5-tasmax-bias-trend-%s-%s.dat'%(region, model), 'wb') as f:
pickle.dump({'yearly_tasmax_max_bias_trend': yearly_tasmax_max_bias_trend,
'yearly_tasmax_max_bias_trend_sig':yearly_tasmax_max_bias_trend_sig,
'yearly_tasmax_monthly_max_bias_trend': yearly_tasmax_monthly_max_bias_trend,
'yearly_tasmax_monthly_max_bias_trend_sig':yearly_tasmax_monthly_max_bias_trend_sig,
'yearly_tasmax_mean_bias_trend':yearly_tasmax_mean_bias_trend,
'yearly_tasmax_mean_bias_trend_sig':yearly_tasmax_mean_bias_trend_sig}, f)