Commit a825a2f4 authored by Carl Schreck's avatar Carl Schreck

Working wave_phase

parent 1f1bbcf1
...@@ -13,7 +13,7 @@ import cjs # '~carl/lib-py' ...@@ -13,7 +13,7 @@ import cjs # '~carl/lib-py'
cjs.tstamp('Here we go!') cjs.tstamp('Here we go!')
# These are some parameters that could be useful to have up top # These are some parameters that could be useful to have up top
wave_name = 'kelvin' wave_name = 'mjo'
var_name = 'olr' var_name = 'olr'
min_lon = 100 min_lon = 100
max_lon = 100 max_lon = 100
...@@ -26,21 +26,30 @@ if var_name == 'olr': ...@@ -26,21 +26,30 @@ if var_name == 'olr':
in_path = os.environ['NASA_YMC'] + '/fcst_verif/olr.k09.nc' in_path = os.environ['NASA_YMC'] + '/fcst_verif/olr.k09.nc'
elif var_name == 'trmm': elif var_name == 'trmm':
in_path = os.environ['NASA_YMC'] + '/fcst_verif/rain.k09.nc' in_path = os.environ['NASA_YMC'] + '/fcst_verif/rain.k09.nc'
ds = xr.open_dataset(in_path) ds_in = xr.open_dataset(in_path)
cjs.tstamp('lon')
if min_lon == max_lon: if min_lon == max_lon:
wave = ds[wave_name].sel(lon=min_lon, method='nearest') wave = ds_in[wave_name].sel(lon=min_lon, method='nearest')
out_path = out_path + '_' + cjs.lon_string(min_lon)
else: else:
wave = ds[wave_name].sel(lon=slice(min_lon,max_lon)).mean(dim='lon') wave = ds_in[wave_name].sel(lon=slice(min_lon,max_lon)).mean(dim='lon')
out_path = (out_path + '_' + cjs.lon_string(min_lon) + '-'
+ cjs.lon_string(max_lon))
cjs.tstamp('lat')
if min_lat == max_lat: if min_lat == max_lat:
wave = wave.sel(lat=min_lat, method='nearest') wave = wave.sel(lat=min_lat, method='nearest')
out_path = out_path + '_' + cjs.lat_string(min_lat)
else: else:
wave = wave.sel(lat=slice(min_lat,max_lat)).mean(dim='lat') wave = wave.sel(lat=slice(min_lat,max_lat)).mean(dim='lat')
out_path = (out_path + '_' + cjs.lat_string(min_lat) + '-'
+ cjs.lat_string(max_lat))
cjs.tstamp('Calculate tendency') cjs.tstamp('Calculate tendency')
tend = wave.differentiate('time') tend = wave.differentiate('time')
tend.name = 'tend'
cjs.tstamp('Standardize both') cjs.tstamp('Standardize both')
...@@ -49,13 +58,17 @@ tend = (tend - tend.mean()) / tend.std() ...@@ -49,13 +58,17 @@ tend = (tend - tend.mean()) / tend.std()
cjs.tstamp('Calculate amplitude and phase') cjs.tstamp('Calculate amplitude and phase')
amp = (wave**2 + tend**2)**0.5 amp = (wave**2 + tend**2)**0.5
amp.name = 'amp'
phase = xr.where((wave >= 0) & (abs(wave) >= abs(tend)), 'H', '') phase = xr.where((wave >= 0) & (abs(wave) >= abs(tend)), 'H', '')
phase = xr.where((wave < 0) & (abs(wave) >= abs(tend)), 'L', phase) phase = xr.where((wave < 0) & (abs(wave) >= abs(tend)), 'L', phase)
phase = xr.where((tend >= 0) & (abs(tend) >= abs(wave)), 'I', phase) phase = xr.where((tend >= 0) & (abs(tend) >= abs(wave)), 'I', phase)
phase = xr.where((tend < 0) & (abs(tend) >= abs(wave)), 'D', phase) phase = xr.where((tend < 0) & (abs(tend) >= abs(wave)), 'D', phase)
phase.name = 'phase'
cjs.tstamp('Write result') cjs.tstamp('Write result')
ds_out = xr.merge([wave, tend, amp, phase])
ds_out.to_netcdf(out_path + '.nc')
df = pd.DataFrame({'time':wave.time, 'wave':wave, 'tend':tend, 'amp':amp, df = pd.DataFrame({'time':wave.time, 'wave':wave, 'tend':tend, 'amp':amp,
'phase':phase}) 'phase':phase})
df.to_csv(out_path + '.csv', index=False, float_format='%6.3f') df.to_csv(out_path + '.csv', index=False, float_format='%6.3f')
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment