-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathec_earth_timeseries_ensemble.py
73 lines (57 loc) · 2.96 KB
/
ec_earth_timeseries_ensemble.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
66
67
68
69
70
71
72
73
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 29 22:28:28 2021
@author: morenodu
"""
def open_regularize(address_file, reference_file):
DS_ec = xr.open_dataset(address_file,decode_times=True)
if list(DS_ec.keys())[1] == 'tasmax':
da_ec = DS_ec[list(DS_ec.keys())[1]] - 273.15
elif list(DS_ec.keys())[1] == 'pr':
da_ec = DS_ec[list(DS_ec.keys())[1]] * 1000
else:
da_ec = DS_ec[list(DS_ec.keys())[1]]
DS_ec = da_ec.to_dataset()
DS_ec_crop = DS_ec.where(reference_file.mean('time') > -300 )
return DS_ec_crop
DS_test2 = xr.open_dataset("EC_earth_2C/pr_m_ECEarth_2C_ensemble_2062-4062.nc",decode_times=True)
def open_regularize_3(address_file, reference_file):
DS_ec = xr.open_dataset(address_file,decode_times=True)
DS_ec['lat']=DS_test2.lat
if list(DS_ec.keys())[1] == 'tasmax':
da_ec = DS_ec[list(DS_ec.keys())[1]] - 273.15
elif list(DS_ec.keys())[1] == 'pr':
da_ec = DS_ec[list(DS_ec.keys())[1]] * 1000
else:
da_ec = DS_ec[list(DS_ec.keys())[1]]
DS_ec = da_ec.to_dataset()
DS_ec_crop = DS_ec.where(reference_file.mean('time') > -300 )
return DS_ec_crop
DS_tmx_ec_un = open_regularize("EC_earth_PD/tasmax_m_ECEarth_PD_ensemble_2035-4035.nc", DS_pre_cru_us['pre'])
DS_tmx_ec = open_regularize("EC_earth_PD/tasmax_m_ECEarth_PD_ensemble_2035-4035.nc", DS_pre_cru_us['pre']).groupby('time.year').mean('time')
DS_tmx_ec_2C = open_regularize("EC_earth_2C/tasmax_m_ECEarth_2C_ensemble_2062-4062.nc", DS_pre_cru_us['pre']).groupby('time.year').mean('time')
DS_tmx_ec_3C = open_regularize_3("EC_earth_3C/tasmax_d_ECEarth_3C_ensemble_2082-4082.nc", DS_pre_cru_us['pre']).groupby('time.year').mean('time')
DS_tmx_ec_mean = DS_tmx_ec.to_dataframe().groupby(['year']).mean()
DS_tmx_ec_mean.index = range((2000))
DS_tmx_ec_2C_mean = DS_tmx_ec_2C.to_dataframe().groupby(['year']).mean()
DS_tmx_ec_2C_mean.index = range((2000))
DS_tmx_ec_3C_mean = DS_tmx_ec_3C.to_dataframe().groupby(['year']).mean()
DS_tmx_ec_3C_mean.index = range((2000))
DS_tmx_ec_mean.plot()
fig1 = plt.figure(figsize=(10,5))
plt.hist(DS_tmx_ec_mean, bins=50, label = 'Present Day')
plt.hist(DS_tmx_ec_2C_mean,bins=50, label = '2C')
plt.hist(DS_tmx_ec_3C_mean,bins=50, label = '3C')
plt.legend(loc="upper left")
plt.title('Histogram of EC-Earth ensembles')
plt.ylabel('Frequency')
plt.xlabel('Mean yearly temperature (°C)')
fig1.savefig('paper_figures/ec_earth_moving_average.png', format='png', dpi=250)
plt.show()
plt.figure(figsize=(20,10)) #plot clusters
ax=plt.axes(projection=ccrs.Mercator())
DS_tmx_ec_un['tasmax'].mean('time').plot(x='lon', y='lat',transform=ccrs.PlateCarree(), robust=True)
ax.add_geometries(us1_shapes, ccrs.PlateCarree(),edgecolor='black', facecolor=(0,1,0,0.0))
# ax.set_extent([-105,-25,-50,50], ccrs.PlateCarree())
ax.set_title('Spatial variability of bias between CRU and EC-earth')
plt.show()