#!/usr/bin/env python # Wenchang Yang (wenchang@princeton.edu) # Thu May 21 11:24:22 EDT 2020 if __name__ == '__main__': from misc.timer import Timer tt = Timer(f'start {__file__}') import sys, os.path, os, glob import xarray as xr, numpy as np, pandas as pd #import matplotlib.pyplot as plt #more imports from xtc import tc_count from xtc.mask import get_landflag from misc.cim import cim, sem #confidence interval of the mean # if __name__ == '__main__': tt.check('end import') #start from here thisdir = os.path.dirname(__file__) source = 'AM2.5_amipHadISST' basin = 'NA' def get_nseed(basin='NA'): #ifile = '/tigress/wenchang/analysis/TC/HIRAM/amipHadISST_tigercpu_intelmpi_18_540PE/netcdf/HIRAM.amipHadISST_tigercpu_intelmpi_18_540PE.tc_tracks.allstorms.1971-2018.nc' #ifile = '/tigress/wenchang/analysis/TC/AM2.5/amipHadISST_tigercpu_intelmpi_18_540PE/netcdf/tc_tracks.allstorms.1971-2018.nc' ifile = '/tigress/wenchang/analysis/rainyseed/AM2.5_amipHadISST_tigercpu_intelmpi_18_540PE_en01-05_yr1971to2018_mo1to12_rainPerc99.5_distThresh2.9_sizeThresh3.1_latLim30.nc' #ofile = os.path.join(thisdir, f'data_nseed_rainy_{basin}.nc') ofile = __file__.replace('.py', f'_{basin}.nc') if os.path.exists(ofile): print('[exists]:', ofile) return xr.open_dataset(ofile)['N_SEED'] seed_kinds = [6, 12, 18, 24] # least hours of life """ ds = (xr.open_dataset(ifile) # genesis locations within the tropics .pipe(lambda x: x.where(x.lat.isel(stage=0)>-30).where(x.lat.isel(stage=0)<30)) ) """ ds = xr.open_dataset(ifile).rename(track='storm', lifetime='stage') #ds = ds.where(ds.vort.max('stage')>4e-4) ds = ds.where(ds.vort.max('stage')>6e-4) n_ocean = get_landflag(ds).pipe(lambda x: x<0.5).sum('stage') # calculate # of seed with different threshold of life hours n_seed = [ds.where(n_ocean>(seed_life//6)).pipe(tc_count, basin=basin) for seed_life in seed_kinds] #if 'vmax' in ds and 'windmax' not in ds: # ds = ds.rename(vmax='windmax') n_seed = xr.concat(n_seed, dim=pd.Index(seed_kinds, name='life')) n_seed = n_seed.rename('N_SEED') #save to nc n_seed.to_dataset().to_netcdf(ofile) print('[saved]:', ofile) return n_seed if __name__ == '__main__': from wyconfig import * #my plot settings fig, ax = plt.subplots() #figname = f'data_nseed_rainy_{basin}_{tt.today()}.png' figname = __file__.replace('.py', f'_{basin}_{tt.today()}.png') da = get_nseed(basin=basin) da.groupby('time.year').sum('time').mean('en').plot(hue='life') plt.title(f'{basin} N_SEED(rainy): {source}') plt.savefig(figname) print('[saved]:', figname) tt.check(f'**Done**') plt.show()