#!/usr/bin/env python # Wenchang Yang (wenchang@princeton.edu) # Wed Dec 14 09:50:06 EST 2022 if __name__ == '__main__': import sys from misc.timer import Timer tt = Timer('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 import geoxarray from misc.seasons import sel_season from misc import get_kws_from_argv # if __name__ == '__main__': tt.check('end import') # #start from here region = get_kws_from_argv('region', 'NorthPP') ifile = glob.glob(f'precip_*_{region}_masked.nc')[0] ofile = ifile.replace('masked.nc', 'index.nc') #https://docs.google.com/document/d/1NvC7xYOEd5ejLpntthWvddubywUZk_mCEDhdAZHIulQ/edit?pli=1#heading=h.5gmk8vl67plq #Working event definition(s): #3-day accumulated precipitation June-Sep (more related to the monsoons) maximum in the Philippines, for all land areas above 12.5 deg N #2-day accumulated precipitation June-September maximum in Taiwan #2-day accumulated precipitation June-September maximum in Hunan province, China season = 'JJAS' N = 3 if region == 'NorthPP' else 2 #rolling window daname = f'rx{N}day' units = 'mm' long_name = f'Rx{N}day precip' func = lambda da: da.geo.fldmean() \ .rolling(time=N).sum() \ .pipe(sel_season, season) \ .groupby('time.year').max('time') \ .assign_attrs(units=units, long_name=long_name) if os.path.exists(ofile): print('[exists]:', ofile) da = xr.open_dataarray(ofile) else: #make index print('loading...') da = xr.open_dataarray(ifile).load() #daname = da.name print('making index...') #rxNday da = da.pipe(func) #print(da); sys.exit() #save print('saving...') da.to_dataset(name=daname).to_netcdf(ofile) print('[saved]:', ofile) if __name__ == '__main__': from wyconfig import * #my plot settings #da.mean('ens', keep_attrs=True).plot() da.plot(hue='ens') #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) tt.check(f'**Done**') print() if 'notshowfig' in sys.argv: pass else: plt.show()