
import xarray as xr
import pandas as pd
import numpy as np
import hvplot.pandas
import hvplot.xarray
import os, datetime, warnings
import cartopy.crs as ccrs
from cesm_utils import fixtime, spd
from my_bokeh_themes import *
warnings.filterwarnings('ignore')
print('Last updated on {}'.format(datetime.datetime.now().ctime()))

Last updated on Wed Sep 15 13:14:43 2021


from pencil_cases import *
case = 'cesm2.B1850.f09_g17.1dpop.006'
datadir = '/home/ivan/Data/Postproc/1DPOP_CESM2/'
# datadir = '/mnt/chromeos/removable/WD/Data/Postproc/1DPOP_CESM2/'
infile = datadir + '{}.abs.max_2D.nc'.format(case)
ds_max_2D = xr.open_dataset(infile, decode_times=False, chunks='auto',
                         drop_variables=['HMXL','TEMP','UVEL','VVEL','ULAT','ULONG','z_w','dz','dzw'])
ds_max_2D = fixtime(ds_max_2D, 1851)
print('{}: {}'.format(case, case_info[case]))

stacked = ds_max_2D.SALT.stack(z=('nlat','nlon'))
maxind = stacked.argmax(axis=1)
maxpos = stacked['z'][maxind]
ntime = ds_max_2D.dims['time']
imax = [maxpos.values[t][0] for t in range(ntime)]
jmax = [maxpos.values[t][1] for t in range(ntime)]

cesm2.B1850.f09_g17.1dpop.006: 1D + G-terms + SSS restoring (tau = 365 days) + EBM off


smax = ds_max_2D.SALT.max(dim='time')
data = []
tlat, tlon = ds_max_2D.TLAT.values, ds_max_2D.TLONG.values

for i, j in sorted(set([p for p in zip(imax,jmax)])):
    data.append([smax.values[i,j], tlat[i,j], tlon[i,j], i, j])
    
df = pd.DataFrame(data, columns=['SALT', 'lat', 'lon', 'i', 'j'])
df


df.hvplot.points('lon', 'lat', geo=True, color='red', marker='x', coastline=True, global_extent=True, frame_width=400, tiles='EsriTerrain')
# df.hvplot.points('lon', 'lat', geo=True, color='red', marker='x', coastline=True, frame_width=400, tiles='EsriTerrain')


# df.hvplot.points?
# from bokeh.sampledata.sea_surface_temperature import sea_surface_temperature as sst
# sst.hvplot()


ds_max = xr.open_dataset(os.path.join(datadir, '{}.abs.max.nc'.format(case)), decode_times=False)
ds_max = fixtime(ds_max, 1851)
ds_max['time'] = ds_max.time.to_index().to_datetimeindex()

data_dict = {vname: ds_max[vname].to_series() for vname in ['HMXL','TEMP','SALT','UVEL','VVEL']}
df_series = pd.DataFrame(data_dict)#.reset_index()

df_series.head()


plist = df_series.hvplot.line(aspect=1, frame_height=350, subplots=True, shared_axes=False, grid=True).cols(3)
for p in plist.items():
#     print(p[0], ds_max[p[0]].units)
    p[1].opts(ylabel=ds_max[p[0]].units)

plist


infile = datadir + '{}.weddell_sea.nc'.format(case)
ds_tseries = xr.open_dataset(infile, decode_times=False)
ds_tseries = fixtime(ds_tseries, 1851)
ds_tseries['time'] = ds_tseries.time.to_index().to_datetimeindex()
ds_tseries = ds_tseries.squeeze(drop=True)


subset = ds_tseries.SALT.sel(time=slice('1930-01-01','1950-12-31'), z_t=slice(500,6.3886262e+04))
df_salt = subset.to_dataframe().reset_index()
# df_salt = ds_tseries.SALT.to_dataframe().reset_index()
df_salt['z_t'] = df_salt.z_t / 100 # cm -> m
df_salt.hvplot.heatmap(x='time', y='z_t', C='SALT', cmap='bgy', height=400, title='Salinity - Weddell Sea',
                       clabel='SALT ({})'.format(ds_tseries.SALT.units)).opts(invert_yaxis=True)


ds_tseries['FRAZIL'] = - ds_tseries.QFLUX/ds_tseries.latent_heat_fusion_mks
ds_tseries['FRAZIL'].attrs = {'long_name':'frazil ice formation', 'units':'kg/m^2/s'}

# - [Kg of salt / m^2 / s] x (1/0.0347) [ Kg of sea water / Kg of salt] x (1/1027) [m^3/kg of sea water] x 1000 [kg of freshwater / m^3]
ds_tseries['SALT_F'] = - ds_tseries.SALT_F * 1/0.0347 * 1/1027 * 1000
ds_tseries['SALT_F'].attrs = {'long_name':'fresh water flux', 'units':'kg/m^2/s'}

varlist = ['SFWF','EVAP_F','PREC_F','ROFF_F','MELT_F','IOFF_F','SALT_F','FRAZIL','IFRAC']
data_dict = {vname: ds_tseries[vname][:].values for vname in varlist}
df_fluxes = pd.DataFrame(data_dict, index=ds_tseries.time)
df_fluxes.index.name = 'time'

df_fluxes.head()


flist = df_fluxes.hvplot.line(aspect=1, frame_height=350, subplots=True, shared_axes=False, grid=True).cols(3)
for p in flist.items():
    p[1].opts(ylabel=ds_tseries[p[0]].units)

flist


# ds_tseries.SALT.hvplot.quadmesh(x='time', y='z_t', rasterize=True)#.opts(invert_yaxis=True)
# foo = ds_tseries.SALT.sel(time=slice('1871-01-01','1890-12-31'), z_t=slice(500,10500))
# foo.hvplot.quadmesh(x='time', y='z_t', rasterize=True)
# foo


# ds_tseries.load()
# ds_tseries.sel(time=ds_tseries.get_index('time').duplicated())

	HMXL	TEMP	SALT	UVEL	VVEL
time
1851-01-01	52526.035156	32.270226	45.964340	136.165497	90.170372
1851-01-02	78700.289062	32.123447	45.939728	95.371529	92.089088
1851-01-03	107081.210938	32.055328	45.932510	82.988930	102.024765
1851-01-04	106187.531250	32.008568	45.922672	74.433929	91.994499
1851-01-05	105998.296875	32.028236	45.913082	71.651810	82.429817

	SFWF	EVAP_F	PREC_F	ROFF_F	MELT_F	IOFF_F	SALT_F	FRAZIL	IFRAC
time
1851-01-01	0.000100	-0.000002	1.553492e-10	7.373395e-09	0.000115	0.000000e+00	-0.000012	-0.0	0.676629
1851-01-02	0.000333	-0.000004	1.564804e-07	4.623667e-09	0.000377	8.600123e-10	-0.000039	-0.0	0.672917
1851-01-03	0.000278	-0.000003	2.343692e-09	7.673421e-09	0.000314	0.000000e+00	-0.000032	-0.0	0.640486
1851-01-04	0.000246	-0.000002	3.036637e-08	7.299467e-09	0.000278	0.000000e+00	-0.000029	-0.0	0.604148
1851-01-05	0.000224	-0.000004	2.818818e-08	5.713488e-09	0.000255	1.100518e-11	-0.000026	-0.0	0.583724

Find location of SALT max value¶

Find position of salinity maximum absolute values¶

Time series of absolute max values¶

Time series at grid point in Weddell Sea¶