Created by Ivan Lima on Mon Nov 1 2021 10:34:26 -0400
import xarray as xr
import numpy as np
import os, datetime, warnings
from cartopy import crs as ccrs
from cartopy import feature as cfeature
from cesm_utils import fixtime
from pencil_cases import *
datadir = '/home/ivan/Data/Postproc/1DPOP_CESM2'
case = 'cesm2.B1850.f09_g17.1dpop.008'
print('{}: {}'.format(case, case_info[case]))
cesm2.B1850.f09_g17.1dpop.008: 1D + G-terms (start year) + SSS restoring (tau = 30 days) + EBM off
import hvplot.xarray
import holoviews as hv
proj = ccrs.Robinson(central_longitude=210)
plist = hv.Layout()
for vname in ['HMXL', 'HBLT', 'IFRAC']:
with xr.open_dataset(os.path.join(datadir, '{}.{}.annual.nc'.format(case,vname)), decode_times=False) as ds:
ds = fixtime(ds, 1851)
ds['time'] = ds.time.to_index().to_datetimeindex()
ds['time'] = ds.time.dt.year
poly = ds[vname].polyfit('time', 1) # fit polynomial of degree 1
slope = poly.polyfit_coefficients[0,...] * 100 # per year -> per century
ds['slope'] = slope
# vmax = np.abs([slope.min().item(), slope.max().item()]).max()
vmax = np.abs([slope.quantile(0.0025,dim=(['nlat','nlon'])).item(),
slope.quantile(0.9975,dim=(['nlat','nlon'])).item()]).max()
clabel = '{} slope ({} per century)'.format(vname, ds[vname].units)
p = ds.slope.hvplot.quadmesh(x='TLONG', y='TLAT', geo=True, coastline=True, projection=proj, project=True,
cmap='coolwarm', clim=(-vmax, vmax), rasterize=True, global_extent=True,
clabel=clabel, title=vname)
plist += p
plist.cols(1)
plist = hv.Layout()
for vname in ['TEMP', 'SALT', 'RHO']:
with xr.open_dataset(os.path.join(datadir, '{}.{}.annual.nc'.format(case,vname)), decode_times=False) as ds:
ds = ds.isel(z_t=0) # select suface values
ds = fixtime(ds, 1851)
ds['time'] = ds.time.to_index().to_datetimeindex()
ds['time'] = ds.time.dt.year
poly = ds[vname].polyfit('time', 1) # fit polynomial of degree 1
slope = poly.polyfit_coefficients[0,...] * 100 # per year -> per century
ds['slope'] = slope
# vmax = np.abs([slope.min().item(), slope.max().item()]).max()
vmax = np.abs([slope.quantile(0.0005,dim=(['nlat','nlon'])).item(),
slope.quantile(0.9995,dim=(['nlat','nlon'])).item()]).max()
clabel = '{} slope ({} per century)'.format(vname, ds[vname].units)
p = ds.slope.hvplot.quadmesh(x='TLONG', y='TLAT', geo=True, coastline=True, projection=proj, project=True,
cmap='coolwarm', clim=(-vmax, vmax), rasterize=True, global_extent=True,
clabel=clabel, title='Surface {}'.format(vname))
plist += p
plist.cols(1)
plist = hv.Layout()
for vname in ['TEMP', 'SALT', 'RHO']:
with xr.open_dataset(os.path.join(datadir, '{}.{}.annual.nc'.format(case,vname)), decode_times=False) as ds:
ds = ds.isel(z_t=slice(0,20)) # upper 200 m
ds = fixtime(ds, 1851)
ds['time'] = ds.time.to_index().to_datetimeindex()
ds['time'] = ds.time.dt.year
poly = ds[vname].weighted(ds.dz).mean(dim='z_t').polyfit('time', 1) # fit polynomial of degree 1
slope = poly.polyfit_coefficients[0,...] * 100 # per year -> per century
ds['slope'] = slope
# vmax = np.abs([slope.min().item(), slope.max().item()]).max()
vmax = np.abs([slope.quantile(0.0005,dim=(['nlat','nlon'])).item(),
slope.quantile(0.9995,dim=(['nlat','nlon'])).item()]).max()
clabel = '{} slope ({} per century)'.format(vname, ds[vname].units)
p = ds.slope.hvplot.quadmesh(x='TLONG', y='TLAT', geo=True, coastline=True, projection=proj, project=True,
cmap='coolwarm', clim=(-vmax, vmax), rasterize=True, global_extent=True,
clabel=clabel, title='{} (upper 200 m)'.format(vname))
plist += p
plist.cols(1)
from my_bokeh_themes import *
# coordinates [348:351,37:38]
# list(zip(ds.z_t.values, ds.SALT.isel(time=0, nlat=349, nlon=37).values))
vname = 'SALT'
ds = xr.open_dataset(os.path.join(datadir, '{}.{}.annual.nc'.format(case,vname)), decode_times=False)
ds = fixtime(ds, 1851)
ds['time'] = ds.time.to_index().to_datetimeindex()
pt1 = ds[vname].isel(z_t=slice(0,39), nlat=349, nlon=37)
p1 = pt1.hvplot.heatmap(x='time', y='z_t', C=vname, cmap='bmy', title='{} (annual means)'.format(vname), clabel=ds[vname].units,
frame_height=400, aspect=1.1, shared_axes=False)
p1.opts(invert_yaxis=True)
pt2 = ds[vname].isel(z_t=0, nlat=349, nlon=37)
p2 = pt2.hvplot.line(x='time', y=vname, title='Surface {} (annual means)'.format(vname), frame_height=400,
aspect=1.1, shared_axes=False, grid=True)
(p1.opts(invert_yaxis=True) + p2).cols(2)
vname = 'TEMP'
ds = xr.open_dataset(os.path.join(datadir, '{}.{}.annual.nc'.format(case,vname)), decode_times=False)
ds = fixtime(ds, 1851)
ds['time'] = ds.time.to_index().to_datetimeindex()
pt1 = ds[vname].isel(z_t=slice(0,39), nlat=349, nlon=37)
p1 = pt1.hvplot.heatmap(x='time', y='z_t', C=vname, cmap='bmy', title='{} (annual means)'.format(vname), clabel=ds[vname].units,
frame_height=400, aspect=1.1, shared_axes=False)
p1.opts(invert_yaxis=True)
pt2 = ds[vname].isel(z_t=0, nlat=349, nlon=37)
p2 = pt2.hvplot.line(x='time', y=vname, title='Surface {} (annual means)'.format(vname), frame_height=400,
aspect=1.1, shared_axes=False, grid=True)
(p1.opts(invert_yaxis=True) + p2).cols(2)
vname = 'RHO'
ds = xr.open_dataset(os.path.join(datadir, '{}.{}.annual.nc'.format(case,vname)), decode_times=False)
ds = fixtime(ds, 1851)
ds['time'] = ds.time.to_index().to_datetimeindex()
pt1 = ds[vname].isel(z_t=slice(0,39), nlat=349, nlon=37)
p1 = pt1.hvplot.heatmap(x='time', y='z_t', C=vname, cmap='bmy', title='{} (annual means)'.format(vname), clabel=ds[vname].units,
frame_height=400, aspect=1.1, shared_axes=False)
p1.opts(invert_yaxis=True)
pt2 = ds[vname].isel(z_t=0, nlat=349, nlon=37)
p2 = pt2.hvplot.line(x='time', y=vname, title='Surface {} (annual means)'.format(vname), frame_height=400,
aspect=1.1, shared_axes=False, grid=True)
(p1.opts(invert_yaxis=True) + p2).cols(2)