#!/usr/bin/env python # Wenchang Yang (wenchang@princeton.edu) # Mon Jul 15 13:55:05 EDT 2024 if __name__ == '__main__': import sys,os try: from misc.timer import Timer tt = Timer(f'[{os.getcwd()}] start ' + ' '.join(sys.argv)) except: pass 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 # if __name__ == '__main__': try: tt.check('end import') except: pass # #start from here daname = 'so4' #ref year year0 = get_kws_from_argv('year0', '1970') year0 = int(year0) #pertub year year = get_kws_from_argv('year', '1990') year = int(year) ifile = '/tigress/wenchang/CM2.5_tiger2_wy/CM2.5/input_wy/input_for_FLOR/aerosol.climatology.nc' ds = xr.open_dataset(ifile) da = ds[daname].sel(time=slice(f'{year0-5}', f'{year0+5}')).sum('sigma_full', keep_attrs=True).mean('time', keep_attrs=True) da0 = da da = ds[daname].sel(time=slice(f'{year-5}', f'{year+5}')).sum('sigma_full', keep_attrs=True).mean('time', keep_attrs=True) with xr.set_options(keep_attrs=True): daa = da - da0 if __name__ == '__main__': from wyconfig import * #my plot settings from geoplots import mapplot fig,ax = plt.subplots() da0.plot.contourf(robust=True, levels=21, cmap='Reds') mapplot() ax.set_title(f'FLOR {daname} annual mean: {year0}') #savefig if 'savefig' in sys.argv or 's' in sys.argv: figname = __file__.replace('.py', f'__{daname}_{year0}.png') if 'overwritefig' in sys.argv or 'o' in sys.argv: wysavefig(figname, overwritefig=True) else: wysavefig(figname) fig,ax = plt.subplots() daa.plot.contourf(robust=True, levels=21) mapplot() ax.set_title(f'FLOR {daname} annual mean anom: {year} minus {year0}') #savefig if 'savefig' in sys.argv or 's' in sys.argv: figname = __file__.replace('.py', f'__{daname}_{year}minus{year0}.png') if 'overwritefig' in sys.argv or 'o' in sys.argv: wysavefig(figname, overwritefig=True) else: wysavefig(figname) try: tt.check(f'**Done**') except: pass print() if 'notshowfig' in sys.argv or 'n' in sys.argv: pass else: if 'plt' in globals(): plt.show()