#!/usr/bin/env python
# Wenchang Yang (wenchang@princeton.edu)
# Wed May 26 21:21:01 EDT 2021
if __name__ == '__main__':
    from misc.timer import Timer
    tt = Timer(f'start {__file__}')
import sys, os.path, os, glob, datetime
import xarray as xr, numpy as np, pandas as pd, matplotlib.pyplot as plt
#more imports
import xfilter
#
if __name__ == '__main__':
    tt.check('end import')
#
#start from here
ifile = 'ssttc.past1000histrcp85.0850-2100.ens06.nc'
lowpass = lambda x: x.filter.lowpass(1/40, dim='year', padtype='constant')
ds = xr.open_dataset(ifile)
da = ds.HU
rng = np.random.default_rng(0)
da_mc = xr.DataArray(rng.poisson(da, size=(10000,)+da.shape), dims=('sample',)+da.dims)
meanv = da.pipe(lowpass)
ci = da_mc.pipe(lowpass).quantile([0.025, 0.975], dim='sample')

 
 
if __name__ == '__main__':
    from wyconfig import * #my plot settings
    g = meanv.plot(col='model', col_wrap=2, figsize=(11,6))
    for ax,model in zip(g.axes.flat,ci.model.values):
        ax.fill_between(meanv.year, ci.sel(model=model, quantile=0.025), ci.sel(model=model, quantile=0.975), alpha=1/4, label='95% CI')
        ax.grid('on')
        ax.axvspan(1850,2005, color='gray', alpha=1/4)
    ax.set_xlim(850,2100)
    ax.set_xticks(range(2100,850,-100))
    ax = g.axes[0,0]
    ax.text(0, 0.98, 'CMIP5 past1000 + historical(gray shaded) + RCP8.5',
        transform=ax.transAxes, ha='left', va='top') 
    ax.legend(loc='lower left')
    #ax.set_yscale('log')
    
    #savefig
    if 'savefig' in sys.argv:
        figname = __file__.replace('.py', f'.png')
        wysavefig(figname)
    tt.check(f'**Done**')
    plt.show()