#!/usr/bin/env python # Wenchang Yang (wenchang@princeton.edu) # Mon Nov 18 11:56:55 EST 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 #ifile = '/projects/w/wenchang/MODEL_OUT/HIRAM/amipLMR2019SST0850ic_tigercpu_intelmpi_18_540PE/en01/analysis_wy/TCIplev/19990101.atmos_month.VI.nc' dsname = get_kws_from_argv('dsname', 'VI') #dsname could be VI, PI, eta, omega500, Vshear, To, z300 if dsname == 'PI': daname = 'vmax' else: daname = dsname model = 'AM2.5C360' expname = 'amipHadISSTrcp45_tigercpu_intelmpi_18_1080PE' #ens = range(16,16+1) #ens = range(17,17+1) ens = range(18,18+1) n_ens = len(ens) years = range(1871,2100+1) ofile = f'{dsname}_{model}_{expname}_{years[0]:04d}-{years[-1]:04d}.nc' if n_ens == 1: ofile = ofile.replace(expname, f'{expname}_ens{ens[0]:02d}') else: ofile = ofile.replace(expname, f'{expname}_{n_ens}ens') if os.path.exists(ofile): print('[exists]:', ofile); sy.exit() #concat over ens dim dss_ens = [] for en in ens: #concat over the year dim dss_time = [] for year in years: #ifile = f'/tigress/wenchang/MODEL_OUT/{model}/{expname}/en{en:02d}/analysis_wy/TCIplev/{year:04d}0101.atmos_month.{dsname}.nc' ifile = f'/scratch/gpfs/GEOCLIM/wenchang/analysis/TC/{model}/{expname}/en{en:02d}/analysis_wy/TCIplev/{year:04d}0101.atmos_month.{dsname}.nc' if dsname in ('To',): ifile = ifile.replace(f'{dsname}.nc', 'PI.nc') #omega500 is under different directory (plev instead of TCIplev) if dsname in ('omega500',): ifile = ifile.replace('TCIplev', 'plev').replace('atmos_month', 'atmos_month.plev') if dsname in ('z300',): ifile = ifile.replace('TCIplev', 'plev').replace('atmos_month', 'atmos_month.plev').replace(f'{dsname}', 'z_full') print(f'ens {en:02d} year {year:04d} {dsname}: {ifile}') ds = xr.open_dataset(ifile) if dsname in ('PI', 'To'): ds = ds[[daname]] #only select the specified variable: either vmax or pmin if dsname in ('omega500',): ds = ds[['omega']] #cause error if keep all variables when loading if dsname in ('z300',): ds = ds[['z_full']].sel(level=300) ds.load() dss_time.append(ds) print('concat over time...') ds = xr.concat(dss_time, dim='time') dss_ens.append(ds) print('concat over ens...') ds = xr.concat(dss_ens, dim=pd.Index(ens, name='ens')) #saveing print('saving...') ds = ds.rename(grid_xt='lon', grid_yt='lat').load() ds.to_netcdf(ofile) print('[saved]:', ofile) if __name__ == '__main__': #from wyconfig import * #my plot settings #savefig if 'savefig' in sys.argv or 's' in sys.argv: figname = __file__.replace('.py', f'.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()