#!/usr/bin/env python # Wenchang Yang (wenchang@princeton.edu) # Wed Mar 1 11:08:25 EST 2023 if __name__ == '__main__': import sys,os from misc.timer import Timer tt = Timer(f'[{os.getcwd()}] 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 from misc import get_kws_from_argv import salem from misc.wysalem import get_world_shape from misc import get_kws_from_argv # if __name__ == '__main__': tt.check('end import') # #start from here region = get_kws_from_argv('region', 'jamaica') shapefile = f'../shapefiles/sf_{region}.shp' #shdf = salem.read_shapefile(shapefile) shdf_raw = salem.read_shapefile(shapefile) #the original shapefile causes error with salem. so extract its geometry values to assign to another shapefile that works shdf = get_world_shape('Jamaica').iloc[0:1, :] shdf['geometry'] = shdf_raw['geometry'].values #ifiles = glob.glob('precip_*_hurricaneHelene.nc') #ifiles.sort() #nfiles = len(ifiles) #if nfiles == 3: #for AM2.5C360 # ofile = ifiles[0].replace('ens06', '3ens').replace('.nc', f'_{region}_masked.nc') #elif nfiles == 10: #for FLOR # ofile = ifiles[0].replace('ens01', '10ens').replace('.nc', f'_{region}_masked.nc') ifile = glob.glob('precip_*_Melissa.nc')[0] ofile = ifile.replace('.nc', f'_{region}_masked.nc') if os.path.exists(ofile): print('[exists]:', ofile) sys.exit() print('mask data...') #ds = xr.open_mfdataset(ifiles).load() \ ds = xr.open_dataset(ifile).load() \ .pipe(lambda x: x.assign_coords(lon=x.lon.values-360)) \ .salem.subset(shape=shdf, margin=1) \ .salem.roi(shape=shdf) print('saving data...') encoding = {vname:{'zlib': True, 'complevel': 1, 'dtype': 'float32'} for vname in list(ds.data_vars)} ds.to_netcdf(ofile, encoding=encoding, unlimited_dims='time') 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) tt.check(f'**Done**') print() if 'notshowfig' in sys.argv or 'n' in sys.argv: pass else: if 'plt' in globals(): plt.show()