#!/usr/bin/env python # Wenchang Yang (wenchang@princeton.edu) # Mon Feb 1 16:59:42 EST 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 #import matplotlib.pyplot as plt #more imports import xfilter from shared import spines_off, emulate_hu from misc.cim import cim # if __name__ == '__main__': tt.check('end import') # #start from here def wyplot(ax=None, legend_on=True, legend_ncol=3): if ax is None: fig, ax = plt.subplots(figsize=(8,3)) n_window = 40 padtype = 'even' cases = ['fullForcing', 'GHG', 'VOLC_GRA', 'LULC_HurttPongratz', 'SSI_VSK_L', 'ORBITAL', 'OZONE_AER']#, '850forcing']#, '0850cntl'] #case info: https://www.cesm.ucar.edu/experiments/cesm1.1/LM/ for case in cases: color = 'k' if case == 'fullForcing' else None lw = 2 if case == 'fullForcing' else 1.5 alpha = 1 if case == 'fullForcing' else 0.5 ifile = os.path.join(os.path.dirname(__file__), 'data', f'LME_TC.{case}.*fullForcingRef.nc') ds = xr.open_mfdataset(ifile, combine='by_coords') if 'en' in ds.dims: #da = ds.HU.mean('en').filter.lowpass(1./n_window, dim='year', padtype=padtype) da = emulate_hu(ds, 'MDR', 'TROP', yearsRef=slice(1982,2005), lmeCase=True) \ .mean('en').filter.lowpass(1./n_window, dim='year', padtype=padtype) n_ens = ds.en.size else: #da = ds.HU.filter.lowpass(1./n_window, dim='year', padtype=padtype) da = emulate_hu(ds, 'MDR', 'TROP', yearsRef=slice(1982,2005), lmeCase=True) \ .filter.lowpass(1./n_window, dim='year', padtype=padtype) n_ens = None if n_ens: label = f'{case}({n_ens})' else: label = case da.plot(ax = ax, color=color, label=label, alpha=alpha, lw=lw) if False and case == 'fullForcing': #da_cim = ds.HU.filter.lowpass(1./n_window, dim='year', padtype=padtype).pipe(cim, confidence=0.95) #da_cim = emulate_hu(ds, 'MDR', 'TROP', yearsRef=slice(1982,2005), lmeCase=True) \ # .filter.lowpass(1./n_window, dim='year', padtype=padtype).pipe(cim, confidence=0.95) #ax.fill_between(da.year, da-da_cim, da+da_cim, color='gray', alpha=0.3) da_ = emulate_hu(ds, 'MDR', 'TROP', yearsRef=slice(1982,2005), lmeCase=True) rng = np.random.default_rng(0) zz = rng.poisson(da_, size=(10000,)+da_.shape) da_mc = xr.DataArray(zz, dims=('sim',)+da_.dims).assign_coords(year=da_.year.values).stack(s=('sim','en')) #se = da_mc.filter.lowpass(1./n_window, dim='year', padtype=padtype).std('s') #ci = da - se, da + se ci = da_mc.filter.lowpass(1./n_window, dim='year', padtype=padtype).quantile([0.025, 0.975], dim='s') #da_mc = xr.DataArray(zz, dims=('sim',)+da_.dims).assign_coords(year=da_.year.values).mean('en') #ci = da_mc.filter.lowpass(1./n_window, dim='year', padtype=padtype).quantile([0.025, 0.975], dim='sim') ax.fill_between(da.year, ci[0], ci[1], color='gray', alpha=0.2) ifile1 = os.path.join(os.path.dirname(__file__), 'data', f'LME_TC.0850cntl.fullForcingRef.nc') ifile2 = os.path.join(os.path.dirname(__file__), 'data', f'LME_TC.850forcing.fullForcingRef.nc') #da_ctl = ( xr.open_dataset(ifile1).HU.mean('year') + xr.open_dataset(ifile2).HU.mean('year') )/2 # average between two ctl runs #da_ctl_std = ( xr.open_dataset(ifile1).HU.filter.lowpass(1./n_window, dim='year', padtype=padtype).std('year') + xr.open_dataset(ifile2).HU.filter.lowpass(1./n_window, dim='year', padtype=padtype).std('year') )/2 # average std between two ctl runs #da_ctl = emulate_hu(xr.open_dataset(ifile1), 'MDR', 'TROP', yearsRef=slice(1982,2005), lmeCase=True) \ # .mean('year') #da_ctl_std = emulate_hu(xr.open_dataset(ifile1), 'MDR', 'TROP', yearsRef=slice(1982,2005), lmeCase=True) \ # .filter.lowpass(1./n_window, dim='year', padtype=padtype).std('year') #lower = da_ctl - da_ctl_std #upper = da_ctl + da_ctl_std da_ = emulate_hu(xr.open_dataset(ifile1), 'MDR', 'TROP', yearsRef=slice(1982,2005), lmeCase=True) rng = np.random.default_rng(0) zz = rng.poisson(da_, size=(10000,)+da_.shape) da_mc = xr.DataArray(zz, dims=('mcrun',)+da_.dims).assign_coords(year=da_.year.values) \ .filter.lowpass(1./n_window, dim='year', padtype=padtype) da_ctl = da_mc.mean(['mcrun', 'year']) lower = da_mc.quantile(0.025, dim=['mcrun', 'year']) upper = da_mc.quantile(0.975, dim=['mcrun', 'year']) #ax.axhline(da_ctl, color='gray', ls='--', label='0850cntl mean and mean$\pm$STD') ax.axhline(da_ctl, color='gray', ls='--', label='0850cntl mean and 95% CI') ax.axhline(lower, color='gray', ls='--') ax.axhline(upper, color='gray', ls='--') #plt.axhspan(da_ctl-da_ctl_std, da_ctl+da_ctl_std, color='gray', alpha=0.3) if legend_on: ax.legend(ncol=legend_ncol, frameon=False) ax.set_xticks(np.arange(2000, 850, -100)) ax.set_xlim(850, 2000) ax.set_ylim(5, 11) spines_off(ax) #ax.set_ylabel('Hurricane #') ax.set_ylabel('LME HU count') if __name__ == '__main__': from wyconfig import * #my plot settings wyplot() """ plt.figure(figsize=(6,6*8/16)) n_window = 40 padtype = 'even' cases = ['fullForcing', 'GHG', 'VOLC_GRA', 'LULC_HurttPongratz', 'SSI_VSK_L', 'ORBITAL', 'OZONE_AER']#, '850forcing']#, '0850cntl'] #case info: https://www.cesm.ucar.edu/experiments/cesm1.1/LM/ for case in cases: color = 'k' if case == 'fullForcing' else None lw = 2 if case == 'fullForcing' else 1.5 alpha = 1 if case == 'fullForcing' else 0.5 ifile = os.path.join(os.path.dirname(__file__), 'data', f'LME_TC.{case}.*fullForcingRef.nc') ds = xr.open_mfdataset(ifile, combine='by_coords') print('[loaded]:', ifile) if 'en' in ds.dims: #da = ds.HU.mean('en').filter.lowpass(1./n_window, dim='year', padtype=padtype) da = emulate_hu(ds, 'MDR', 'TROP', yearsRef=slice(1982,2005), lmeCase=True) \ .mean('en').filter.lowpass(1./n_window, dim='year', padtype=padtype) else: #da = ds.HU.filter.lowpass(1./n_window, dim='year', padtype=padtype) da = emulate_hu(ds, 'MDR', 'TROP', yearsRef=slice(1982,2005), lmeCase=True) \ .filter.lowpass(1./n_window, dim='year', padtype=padtype) da.plot(color=color, label=case, alpha=alpha, lw=lw) if case == 'fullForcing': #da_cim = ds.HU.filter.lowpass(1./n_window, dim='year', padtype=padtype).pipe(cim, confidence=0.95) da_cim = emulate_hu(ds, 'MDR', 'TROP', yearsRef=slice(1982,2005), lmeCase=True) \ .filter.lowpass(1./n_window, dim='year', padtype=padtype).pipe(cim, confidence=0.95) plt.fill_between(da.year, da-da_cim, da+da_cim, color='gray', alpha=0.3) ifile1 = os.path.join(os.path.dirname(__file__), 'data', f'LME_TC.0850cntl.fullForcingRef.nc') ifile2 = os.path.join(os.path.dirname(__file__), 'data', f'LME_TC.850forcing.fullForcingRef.nc') #da_ctl = ( xr.open_dataset(ifile1).HU.mean('year') + xr.open_dataset(ifile2).HU.mean('year') )/2 # average between two ctl runs #da_ctl_std = ( xr.open_dataset(ifile1).HU.filter.lowpass(1./n_window, dim='year', padtype=padtype).std('year') + xr.open_dataset(ifile2).HU.filter.lowpass(1./n_window, dim='year', padtype=padtype).std('year') )/2 # average std between two ctl runs #da_ctl = xr.open_dataset(ifile1).HU.mean('year') #da_ctl_std = xr.open_dataset(ifile1).HU.filter.lowpass(1./n_window, dim='year', padtype=padtype).std('year') da_ctl = emulate_hu(xr.open_dataset(ifile1), 'MDR', 'TROP', yearsRef=slice(1982,2005), lmeCase=True) \ .mean('year') da_ctl_std = emulate_hu(xr.open_dataset(ifile1), 'MDR', 'TROP', yearsRef=slice(1982,2005), lmeCase=True) \ .filter.lowpass(1./n_window, dim='year', padtype=padtype).std('year') plt.axhline(da_ctl, color='gray', ls='--', label='0850cntl mean and mean$\pm$STD') plt.axhline(da_ctl - da_ctl_std, color='gray', ls='--') plt.axhline(da_ctl + da_ctl_std, color='gray', ls='--') #plt.axhspan(da_ctl-da_ctl_std, da_ctl+da_ctl_std, color='gray', alpha=0.3) plt.legend(ncol=2, frameon=False) plt.xticks(np.arange(2000, 850, -100)) plt.xlim(850, 2000) plt.ylim(None, 11) spines_off(plt.gca()) plt.ylabel('Hurricane #') """ figname = __file__.replace('.py', f'.png') if len(sys.argv) > 1 and sys.argv[1] == 'savefig': wysavefig(figname) else: print(f'figure not saved. to save figure: python {__file__} savefig') tt.check(f'**Done**') plt.show()