Plotting a histogram based on maximum temperature frequencies

histogram.py

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]

data = xarray.open_dataset('/g/data3/w97/dc8106/AMZ_def_EXPs/test/tasmax_sc_001GPsc_E0_test.nc', chunks={'time':1000})
tasmax = data.tasmax
#tasmin = data.tasmin

lat = data.lat
lon = data.lon
lons,lats = np.meshgrid(lon,lat)

ind_label = indices

#x= np.arange(-4,4,0.001)                                                                                                                                                                                   \
#plt.suptitle(ind_label +' in 121GPsc_E0', fontsize=16)                                                                                                                                                      
#plt.savefig('/g/data3/w97/dc8106/images/'+ind_label+'_ensmean_121GPsc_E0', format='png')                                                                                                                   \

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"]), 10, normed=1, facecolor='green', alpha=0.75)

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--', linewidth=1)

plt.xlabel(indices)
# 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()