Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import xarray
- from scipy.stats import norm
- import numpy as np
- import matplotlib.mlab as mlab
- import matplotlib.pyplot as plt
- list_indices = ['tasmax']
- indices = list_indices[0]
- exp = ['CTL_E0', '121GPsc_E0']
- exp1 = '121GPsc_E0'
- for i,ind in enumerate(exp):
- #print ind
- #files = sorted(glob.glob('/g/data3/w97/dc8106/AMZ_def_EXPs/'+exp+'/tasmax_sc-only_1978-2011_'+exp+'.nc', chunks={'time':1000})
- #analysis/ensmean/'+indices+'/'+indices+'**_*diff_to_CTL.nc'))
- data = xarray.open_dataset('/g/data3/w97/dc8106/AMZ_def_EXPs/'+exp[i]+'/tasmax_sc-only_1978-2011_'+exp[i]+'.nc', chunks={'time':1000})
- #data1 = xarray.open_dataset('/g/data3/w97/dc8106/AMZ_def_EXPs/'+exp1+'/tasmax_sc-only_1978-2011_'+exp1+'.nc', chunks={'time':1000})
- tasmax = data.tasmax - 272.15
- #tasmax1 = data1.tasmax - 272.15
- #tasmin = data.tasmin
- lat = data.lat
- lon = data.lon
- lons,lats = np.meshgrid(lon,lat)
- ind_label = indices
- print(tasmax)
- print("tasmax")
- print(tasmax.stack(dim=["lat","lon","time"]))
- mu, sigma = tasmax.mean().values, tasmax.std().values
- # Print the values of mu and sigma which forces them to be evaluated so I can see how long it takes to do this, then I can tune the time chunking
- print(mu,sigma)
- # the histogram of the data \
- n, bins, patches = plt.hist(tasmax.stack(dim=["lat","lon","time"]), bins = 1000, normed=1, facecolor='green', alpha=0.75)
- plt.xticks(np.arange(20, 50, 2.0))
- print(n)
- print(bins)
- print(patches)
- # add a 'best fit' line \
- y = mlab.normpdf( bins, mu, sigma)
- print(y)
- l = plt.plot(bins, y, 'r--', label=ind, linewidth=10
- l_legend = plt.legend(handles=[l], loc=1)
- l1 = plt.plot(bins, y, 'b--', label=ind, linewidth=1)
- l1_legend = plt.legend(handles=[l1], loc=1
- units = 'Celsius'
- plt.axis([20, 50, 0, 0.15])
- plt.xlabel(indices+' in '+units)
- plt.suptitle(ind_label+ ' in ' +ind, fontsize=16)
- #plt.savefig('/g/data3/w97/dc8106/images/'+ind_label+'_'+exp1, format='png')
- # plt.ylabel('Probability')
- # plt.title(r'$\mathrm{Histogram of '+indices+'\ \mu='+mu+',\ \sigma='+sigma+')')
- #plt.axis([25, 45, 0, 0.03]) \
- #plt.grid(True) \
- plt.show()
Add Comment
Please, Sign In to add comment