#!/usr/bin/env python # Wenchang Yang (wenchang@princeton.edu) # Sat Dec 18 21:22:10 EST 2021 #noIO: exclude Indian Ocean (60-90E) if __name__ == '__main__': import sys from misc.timer import Timer s = ' ' tt = Timer(f'start {s.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 from misc.landmask import flagland import xlinregress # if __name__ == '__main__': tt.check('end import') # #start from here #cal the weights used to merge NH and SH to global da = xr.open_dataarray('/tigress/wenchang/data/cmip6/variables/ts/hist-GHG/wy_regrid_all_members/ts.hist-GHG.ACCESS-CM2.r1i1p1f1.1850-2020.nc') is_ocean = flagland(da) < 0.1 is_ocean = is_ocean.where((is_ocean.lon>90)|(is_ocean.lon<60)) #exclude 60-90E NH = is_ocean.sel(lat=slice(10,30), lon=slice(40,360-20)) #NH = is_ocean.sel(lat=slice(10,30), lon=slice(90,360-20)) #noIO NH = NH.weighted(np.cos(np.deg2rad(NH.lat))).sum(['lat', 'lon']) SH = is_ocean.sel(lat=slice(-30,-10), lon=slice(30,360-150)) #SH = is_ocean.sel(lat=slice(-30,-10), lon=slice(90,360-150)) #noIO SH = SH.weighted(np.cos(np.deg2rad(SH.lat))).sum(['lat', 'lon']) wgtSH = 1/(NH/SH + 1) wgtNH = 1 - wgtSH wgts = (wgtNH.item(), wgtSH.item()) print('weights of NH and SH: ', wgts) #savedata = True if len(sys.argv)>1 and 'savedata' in sys.argv[1:] else False def get_data(dsname): #dsname = 'PI' print(dsname) if dsname == 'PI': daname = dsname long_name = 'Potential Intensity' units = 'm/s' elif dsname == 'SST': daname = 'ts' long_name = dsname units = 'K' elif dsname == 'RH': daname = 'hur' long_name = 'Relative Humidity' units = '%' elif dsname == 'Vshear': daname = dsname long_name = 'Vertical Wind Shear' units = 'm/s' trend_units = f'{units}/year' len_season = 3 if len_season == 3: ifile_hist_NH = [f'{os.path.join(dsname, f)}' for f in os.listdir(dsname) if f.startswith(f'{daname}_historical_NH_ASO_') and f.endswith('_noIO.nc')][0] ifile_hist_SH = [f'{os.path.join(dsname, f)}' for f in os.listdir(dsname) if f.startswith(f'{daname}_historical_SH_FMA_') and f.endswith('_noIO.nc')][0] ifile_nat_NH = [f'{os.path.join(dsname, f)}' for f in os.listdir(dsname) if f.startswith(f'{daname}_hist-nat_NH_ASO_') and f.endswith('_noIO.nc')][0] ifile_nat_SH = [f'{os.path.join(dsname, f)}' for f in os.listdir(dsname) if f.startswith(f'{daname}_hist-nat_SH_FMA_') and f.endswith('_noIO.nc')][0] ifile_ghg_NH = [f'{os.path.join(dsname, f)}' for f in os.listdir(dsname) if f.startswith(f'{daname}_hist-GHG_NH_ASO_') and f.endswith('_noIO.nc')][0] ifile_ghg_SH = [f'{os.path.join(dsname, f)}' for f in os.listdir(dsname) if f.startswith(f'{daname}_hist-GHG_SH_FMA_') and f.endswith('_noIO.nc')][0] #elif len_season == 5: # ifile_nat_NH = [f for f in os.listdir() if f.startswith(f'{daname}_hist-nat_NH_JJASO_') and f.endswith('.nc')][0] # ifile_nat_SH = [f for f in os.listdir() if f.startswith(f'{daname}_hist-nat_SH_DJFMA_') and f.endswith('.nc')][0] # ifile_ghg_NH = [f for f in os.listdir() if f.startswith(f'{daname}_hist-GHG_NH_JJASO_') and f.endswith('.nc')][0] # ifile_ghg_SH = [f for f in os.listdir() if f.startswith(f'{daname}_hist-GHG_SH_DJFMA_') and f.endswith('.nc')][0] print(f'{ifile_hist_NH = }') print(f'{ifile_hist_SH = }') print(f'{ifile_nat_NH = }') print(f'{ifile_nat_SH = }') print(f'{ifile_ghg_NH = }') print(f'{ifile_ghg_SH = }') #historical da_hist = xr.open_dataarray(ifile_hist_NH)*wgts[0] + xr.open_dataarray(ifile_hist_SH)*wgts[1] da_hist = da_hist.sel(year=slice(1982, 2014)) da_hist = da_hist.assign_attrs(units=units, long_name=long_name) lintrend_hist = da_hist.linregress.on(da_hist.year) lintrend_hist.slope.attrs['units'] = trend_units #nat da_nat = xr.open_dataarray(ifile_nat_NH)*wgts[0] + xr.open_dataarray(ifile_nat_SH)*wgts[1] da_nat = da_nat.sel(year=slice(1982, 2014)) da_nat = da_nat.assign_attrs(units=units, long_name=long_name) lintrend_nat = da_nat.linregress.on(da_nat.year) lintrend_nat.slope.attrs['units'] = trend_units #ghg da_ghg = xr.open_dataarray(ifile_ghg_NH)*wgts[0] + xr.open_dataarray(ifile_ghg_SH)*wgts[1] da_ghg = da_ghg.sel(year=slice(1982, 2014)) da_ghg = da_ghg.assign_attrs(units=units, long_name=long_name) lintrend_ghg = da_ghg.linregress.on(da_ghg.year) lintrend_ghg.slope.attrs['units'] = trend_units #savedata #if savedata: # da_nat.drop(['model', 'members']).to_dataframe().to_csv(__file__.replace('.py', '.hist-nat.timeseries.csv')) # lintrend_nat.slope.to_dataframe().to_csv(__file__.replace('.py', '.hist-nat.lintrend.csv')) # da_ghg.drop(['model', 'members']).to_dataframe().to_csv(__file__.replace('.py', '.hist-GHG.timeseries.csv')) # lintrend_ghg.slope.to_dataframe().to_csv(__file__.replace('.py', '.hist-GHG.lintrend.csv')) # print('**csv data saved**') ds_hist = lintrend_hist ds_hist['timeseries'] = da_hist ds_nat = lintrend_nat ds_nat['timeseries'] = da_nat ds_ghg = lintrend_ghg ds_ghg['timeseries'] = da_ghg #save data #historical ofile_hist = f'{dsname}.global.tcseason{len_season}.historical.noIO.nc' if os.path.exists(ofile_hist): print('[exists]:', ofile_hist) else: ds_hist.to_netcdf(ofile_hist) print('[saved]:', ofile_hist) #nat ofile_nat = f'{dsname}.global.tcseason{len_season}.hist-nat.noIO.nc' if os.path.exists(ofile_nat): print('[exists]:', ofile_nat) else: ds_nat.to_netcdf(ofile_nat) print('[saved]:', ofile_nat) #ghg ofile_ghg = f'{dsname}.global.tcseason{len_season}.hist-GHG.noIO.nc' if os.path.exists(ofile_ghg): print('[exists]:', ofile_ghg) else: ds_ghg.to_netcdf(ofile_ghg) print('[saved]:', ofile_ghg) if __name__ == '__main__': #from wyconfig import * #my plot settings dsnames = ('PI', 'SST', 'RH', 'Vshear') for dsname in dsnames: get_data(dsname) #savefig if len(sys.argv)>1 and 'savefig' in sys.argv[1:]: figname = __file__.replace('.py', f'_{daname}_season{len_season}mon.png') wysavefig(figname) tt.check(f'**Done**') #plt.show()