#!/usr/bin/env python # Wenchang Yang (wenchang@princeton.edu) # Fri May 28 11:01:11 EDT 2021 if __name__ == '__main__': from misc.timer import Timer tt = Timer(f'start {__file__}') import sys, os.path, os, glob, datetime import xarray as xr, numpy as np, pandas as pd, matplotlib.pyplot as plt #more imports # if __name__ == '__main__': tt.check('end import') # #start from here ofile = __file__.replace('.py', '.nc').replace('wyget_', '') if os.path.exists(ofile): ds = xr.open_dataset(ofile) print('[loaded]:', ofile) else: #recorded HU ifile = 'recorded_hurricanes.nc' da = xr.open_dataset(ifile)['HU'].rename(YRT='year') da = da.assign_attrs(long_name='Recorded HU #') HU_raw = da #HU correction ifile = 'hu_bootstrap_correction.nc' da = xr.open_dataset(ifile)['BOOTMISS40'].rename(YAYA='year', XAX1_10000='sample') da = da.assign_attrs(long_name='HU # correction') dHU = da #HU adjusted HU = HU_raw + dHU HU = HU.assign_attrs(long_name='Adjusted HU #') #landfall HU ifile = 'landfall.nc' da = xr.open_dataset(ifile)['USHU'].rename(TYEAR='year') da = da.assign_attrs(long_name='Landfall HU #') HUlandfall = da #recorded MH ifile = 'obs_major.nc' da = xr.open_dataset(ifile)['OBS_MAJOR_ON_Y'].rename(YOO='year') da = da.assign_attrs(long_name='Recorded MH #') MH_raw = da #MH correction ifile = 'mh_bootstrap_correction.nc' da = xr.open_dataset(ifile)['BOOTMISS'].rename(YAYA='year', XAX1_10000='sample') da = da.assign_attrs(long_name='MH # correction') dMH = da #MH adjusted MH = MH_raw + dMH MH = MH.assign_attrs(long_name='Adjusted MH #') #landfall MH ifile = 'landfall.nc' da = xr.open_dataset(ifile)['USMA'].rename(TYEAR='year') da = da.assign_attrs(long_name='Landfall MH #') MHlandfall = da #dataset ds = xr.Dataset(dict(HU=HU, HU_raw=HU_raw, dHU=dHU, HUlandfall=HUlandfall, MH=MH, MH_raw=MH_raw, dMH=dMH, MHlandfall=MHlandfall)) #saving ds.to_netcdf(ofile) print('[saved]:', ofile) if __name__ == '__main__': """ from wyconfig import * #my plot settings import xfilter plt.close() lowpass = lambda x: x.filter.lowpass(1/15, dim='year', padtype='constant') fig,ax = plt.subplots(figsize=(8,4)) #raw recorded #ds.HU_raw.plot(color='gray') ds.HU_raw.pipe(lowpass).plot(color='gray', label='HU recorded 15-year lowpass') #adjusted q = ds.HU.quantile([0.025, 0.975], dim='sample') ax.fill_between(q.year, q[0], q[1], alpha=0.5) ds.HU.mean('sample').plot(label='HU adjusted') #lowpass q = ds.HU.pipe(lowpass).quantile([0.025, 0.975], dim='sample') ax.fill_between(q.year, q[0], q[1], alpha=0.5) ds.HU.pipe(lowpass).mean('sample').plot(label='HU adjusted 15-year lowpass') ax.set_xlim(1850, 2020) ax.set_xticks(range(1850,2021,10)) ax.legend(frameon=True, loc='upper center') ax.grid('on') """ #savefig if 'savefig' in sys.argv: figname = __file__.replace('.py', f'.png') wysavefig(figname) tt.check(f'**Done**') plt.show()