#!/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 could be: VI, omega500, eta(absolute vorticity) dsname = get_kws_from_argv('dsname', 'VI') if dsname == 'PI': daname = 'vmax' elif dsname == 'omega500': daname = 'w' else: daname = dsname #model = 'AM2.5C360' expname = 'ERA5' #ens = range(16,16+1) #n_ens = len(ens) years = range(1979,2024+1) ofile = f'{dsname}_{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') #concat over ens dim dss_ens = [] #for en in ens: if True: #concat over the year dim dss_time = [] for year in years: #ifile = f'/scratch/gpfs/GEOCLIM/wenchang/analysis/TC/{model}/{expname}/en{en:02d}/analysis_wy/TCIplev/{year:04d}0101.atmos_month.{dsname}.nc' #ifile = f'/tigress/wenchang/MODEL_OUT/{model}/{expname}/en{en:02d}/analysis_wy/TCIplev/{year:04d}0101.atmos_month.{dsname}.nc' ifile = f'/projects/w/wenchang/data/era5/analysis_wy/TCI_v202201/era5.monthly.{year:04d}.{dsname}.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 ('omega500',): ifile = f'/projects/w/wenchang/data/era5/raw/monthly/plevels/vertical_velocity/era5.vertical_velocity.monthly.{year:04d}.nc' #print(f'ens {en:02d} year {year:04d} {dsname}: {ifile}') print(f'{expname} year {year:04d} {dsname}: {ifile}') ds = xr.open_dataset(ifile) #if dsname == 'PI': 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 ('omega500',): ds = ds[['w']] #cause error if keep all variables when loading #ERA5 omega has name w ds = ds[[daname]] #only keep the data variable if 'expver' in ds.coords: ds = ds.drop('expver') #drop the coords that cause problems if 'valid_time' in ds.dims: ds = ds.rename(valid_time='time') print('[renamed]: valid_time -> time') if 'pressure_level' in ds.dims: levels_new = ds.pressure_level.values.astype(int) ds = ds.rename(pressure_level='level').assign_coords(level=levels_new) print('[converted]: pressure_level -> level') if dsname in ('omega500',): ds = ds.sel(level=500) #only select the 500hPa level 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 = ds.rename(longitude='lon', latitude='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()