#!/usr/bin/env python # Wenchang Yang (wenchang@princeton.edu) # Wed Jul 12 13:23:03 EDT 2023 if __name__ == '__main__': import sys,os from misc.timer import Timer tt = Timer(f'[{os.getcwd()}] start ' + ' '.join(sys.argv)) import sys, os.path, os, glob, datetime import xarray as xr, numpy as np, pandas as pd, matplotlib.pyplot as plt #more imports wython = '/tigress/wenchang/wython' if wython not in sys.path: sys.path.append(wython); print('added to python path:', wython) #from misc import get_kws_from_argv from misc.seasons import sel_season import xfilter from wyplot_lines_LMR import das as das_LMR # if __name__ == '__main__': tt.check('end import') # #start from here model = 'AM2.1' lmr = """ 2018 2019 """.split()[-1] season = 'annual' tspan = [ None, slice('1551', '1650'), ][-1] def get_data(region, season=season, tspan=tspan): if region == 'North China': lonlatbox = '106.120.33.40' elif region == 'Northeast China': lonlatbox = '122.130.40.46' elif region == 'Mongolia': lonlatbox = '107.122.40.47' ifile = f'precip_AM2.1_amipLMR{lmr}SST0850_ens_tigercpu_intelmpi_18_30PE_5ens_0850-1999_fldmean.{lonlatbox}.nc' print(ifile) ### if tspan is None: da = xr.open_dataarray(ifile).load() #only the 5 ens of the last millennium else: #concatenate two exps along the ens dim da = xr.open_dataarray(ifile).sel(time=tspan) da_850 = da ifile = ifile.replace('0850-1999', '1551-1650').replace('0850', '1551') print(ifile) da = xr.open_dataarray(ifile).sel(time=tspan) da = da.assign_coords(ens=range(6,11)) da_1551 = da da = xr.concat([da_850, da_1551], dim='ens') daa = da.groupby('time.month') - da.groupby('time.month').mean('time').mean('ens') daa = daa.pipe(sel_season, season).groupby('time.year').mean('time') return daa das = dict() region = 'North China' das[region] = get_data(region=region) region = 'Northeast China' das[region] = get_data(region=region) region = 'Mongolia' das[region] = get_data(region=region) dass = dict() dass['AM2.1'] = das dass['LMR'] = das_LMR def wyplot(ax=None, region='North China', source='AM2.1', **kws): if ax is None: fig,ax = plt.subplots() #lp = lambda x: x.rolling(year=9, center=True, min_periods=1).mean() n_window = 99 if tspan is None else 9 lp = lambda x: x.filter.lowpass(1/n_window, dim='year', padtype='even') da = dass[source][region] da_lp = da.pipe(lp) mm = da_lp.mean('ens') err = da_lp.std('ens')/np.sqrt(da_lp.ens.size) ax.fill_between(mm.year, mm-err, mm+err, alpha=0.1, **kws) mm.plot(ax=ax, label=region+', '+source, **kws) ax.axhline(0, color='gray', ls='--') ax.axvspan(1625, 1644, color='gray', alpha=0.1) ax.set_ylabel(f'mm/day') ax.set_title(f'{season} precip anom, {n_window}LP') if __name__ == '__main__': from wyconfig import * #my plot settings region = 'North China' fig, ax = plt.subplots() wyplot(region='North China', source='LMR', ax=ax, color='k') wyplot(region='North China', source='AM2.1', ax=ax, color='C0') ax.legend() #savefig if 'savefig' in sys.argv or 's' in sys.argv: figname = __file__.replace('.py', f'__{region.replace(" ", "")}.png') if tspan is None: figname = figname.replace('.png', '_LM.png') #last millennium results, only 5 ens if 'overwritefig' in sys.argv or 'o' in sys.argv: wysavefig(figname, overwritefig=True) else: wysavefig(figname) region = 'Northeast China' fig, ax = plt.subplots() wyplot(region=region, source='LMR', ax=ax, color='k') wyplot(region=region, source='AM2.1', ax=ax, color='C1') ax.legend() #savefig if 'savefig' in sys.argv or 's' in sys.argv: figname = __file__.replace('.py', f'__{region.replace(" ", "")}.png') if tspan is None: figname = figname.replace('.png', '_LM.png') #last millennium results, only 5 ens if 'overwritefig' in sys.argv or 'o' in sys.argv: wysavefig(figname, overwritefig=True) else: wysavefig(figname) region = 'Mongolia' fig, ax = plt.subplots() wyplot(region=region, source='LMR', ax=ax, color='k') wyplot(region=region, source='AM2.1', ax=ax, color='C2') ax.legend() #savefig if 'savefig' in sys.argv or 's' in sys.argv: figname = __file__.replace('.py', f'__{region.replace(" ", "")}.png') if tspan is None: figname = figname.replace('.png', '_LM.png') #last millennium results, only 5 ens if 'overwritefig' in sys.argv or 'o' in sys.argv: wysavefig(figname, overwritefig=True) else: wysavefig(figname) tt.check(f'**Done**') print() if 'notshowfig' in sys.argv or 'n' in sys.argv: pass else: if 'plt' in globals(): plt.show()