Climatological global means¶

Created by Ivan Lima on Thu Aug 25 2022 09:39:37 -0400

In [1]:
import xarray as xr
import pandas as pd
import numpy as np
import holoviews as hv
import os, datetime, my_bokeh_themes, warnings
import hvplot.pandas
import hvplot.xarray
from cartopy import crs as ccrs
from cartopy import feature as cfeature
from cesm_utils import spd
print('Last updated on {}'.format(datetime.datetime.now().ctime()))

Last updated on Fri Sep  2 10:42:00 2022
In [2]:
from bokeh.plotting import show
hv.extension('bokeh')
warnings.filterwarnings('ignore')
In [3]:
datadir = '/home/ivan/Data/Postproc/1DPOP_CESM2'

def fixtime(ds_in):
    ds_out = ds_in.assign(time = ds_in.time.values - 1)
    ds_out['time'] = ds_out.time.assign_attrs(**ds_in.time.attrs)
    return xr.decode_cf(ds_out)

def get_weights(ds, vname):
    """Get weights & dimensions to compute weighted means"""
    if 'z_t' in ds[vname].dims:
        weights = ds.TAREA * ds.dz
        dims = ['nlat','nlon','z_t']
    else:
        weights = ds.TAREA
        dims = ['nlat','nlon']

    weights = weights.where(ds[vname].isel(time=0).notnull())
    return weights, dims

def global_mean(ds, vname, upper200m=False, surf=False):
    if upper200m:
        ds = ds.isel(z_t=slice(0,20)) # upper 200 m
    elif surf:
        ds = ds.isel(z_t=0) # surface
    weights, dims = get_weights(ds, vname)
    da = ds[vname].weighted(weights.fillna(0)).mean(dim=dims)
    da.attrs = ds[vname].attrs
    da.time.attrs = ds.time.attrs
    return da

PI control run¶

In [4]:
ds_picontrol = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.ts_global.nc'), decode_times=False)
ds_picontrol = fixtime(ds_picontrol)
df_picontrol = ds_picontrol[['SHF','TEMP']].groupby(ds_picontrol.time.dt.month).mean().to_dataframe()

ds_sst_picontrol = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.TEMP.surf.nc'),
                                   decode_times=False, chunks={'time':12})
ds_sst_picontrol = fixtime(ds_sst_picontrol)
da_sst = global_mean(ds_sst_picontrol, 'TEMP', surf=True)

units = dict(SHF = ds_picontrol.SHF.units, TEMP = ds_picontrol.TEMP.units, SST = ds_sst_picontrol.TEMP.units)

df_picontrol['SST'] = da_sst.groupby(da_sst.time.dt.month).mean().to_series()
df_picontrol[['SHF', 'SST', 'TEMP']]
Out[4]:
SHF SST TEMP
month
1 9.908335 18.470272 4.224045
2 11.123707 18.568558 4.225834
3 8.283001 18.553994 4.227527
4 0.567431 18.469868 4.228333
5 -10.351532 18.399188 4.227492
6 -21.351665 18.354504 4.224697
7 -21.680179 18.302807 4.220732
8 -10.186318 18.270044 4.217817
9 1.832370 18.202444 4.217120
10 10.142300 18.131908 4.218244
11 12.324347 18.154001 4.220344
12 10.543674 18.287396 4.222405
In [5]:
# plist = df_picontrol[['SHF','SST']].hvplot(aspect=1.1, frame_height=350, grid=True, subplots=True).opts(shared_axes=False)
# for p in plist.items():
#     p[1].opts(ylabel=units[p[0]])
# plist

1DPOP.010 run¶

In [6]:
ds_rest = xr.open_dataset(os.path.join(datadir, 'cesm2.B1850.f09_g17.1dpop.010.ts_global.nc'), decode_times=False)
ds_rest = fixtime(ds_rest)
ds_rest = ds_rest.drop('z_t')
df_rest = ds_rest[['SHF', 'SST', 'TEMP']].groupby(ds_rest.time.dt.month).mean().to_dataframe()
df_rest
Out[6]:
SHF SST TEMP
month
1 3.899234 18.468780 4.223508
2 4.034639 18.561812 4.225195
3 1.805531 18.525048 4.226832
4 -5.765069 18.410479 4.227540
5 -16.902548 18.326183 4.226625
6 -28.366549 18.279786 4.223813
7 -29.181880 18.241320 4.220029
8 -18.655397 18.217971 4.217273
9 -7.310988 18.128286 4.216533
10 1.976833 18.048827 4.217599
11 5.836022 18.095667 4.219670
12 4.915738 18.273954 4.221781
In [7]:
# plist = df_rest.hvplot(aspect=1.1, frame_height=350, grid=True, subplots=True).opts(shared_axes=False)
# for p in plist.items():
#     p[1].opts(ylabel=units[p[0]])
# plist

1DPOP.015 run¶

In [8]:
ds_rest2 = xr.open_dataset(os.path.join(datadir, 'cesm2.B1850.f09_g17.1dpop.015.ts_global.nc'), decode_times=False)
ds_rest2 = fixtime(ds_rest2)
ds_rest2 = ds_rest2.drop('z_t')
df_rest2 = ds_rest2[['SHF', 'SST', 'TEMP']].groupby(ds_rest2.time.dt.month).mean().to_dataframe()
df_rest2
Out[8]:
SHF SST TEMP
month
1 3.529463 18.441522 4.222950
2 2.584329 18.536024 4.224578
3 1.346435 18.487939 4.226164
4 -6.096360 18.364811 4.226965
5 -17.146361 18.271832 4.226070
6 -29.101312 18.214796 4.223156
7 -29.762764 18.172787 4.219276
8 -19.760546 18.145294 4.216412
9 -7.900819 18.049781 4.215638
10 1.814320 17.965044 4.216770
11 6.832998 18.032124 4.219086
12 4.137487 18.227906 4.221193

Year 1161 of PI control run¶

In [10]:
ds_1161_shf = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.SHF.yr_1161.nc'), decode_times=False)
ds_1161_shf = fixtime(ds_1161_shf)
da_1161_shf = global_mean(ds_1161_shf, 'SHF')
da_1161_shf['time'] = range(1,13)
da_1161_shf = da_1161_shf.rename({'time':'month'})

ds_1161_temp = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.TEMP.yr_1161.nc'), decode_times=False)
ds_1161_temp = fixtime(ds_1161_temp)
da_1161_sst = global_mean(ds_1161_temp, 'TEMP', surf=True)
da_1161_temp = global_mean(ds_1161_temp, 'TEMP')
da_1161_sst['time'] = range(1,13)
da_1161_sst = da_1161_sst.rename({'time':'month'})
da_1161_temp['time'] = range(1,13)
da_1161_temp = da_1161_temp.rename({'time':'month'})

df_1161 = pd.DataFrame({'SHF': da_1161_shf.to_series(),
                        'SST': da_1161_sst.to_series(),
                        'TEMP': da_1161_temp.to_series()})
df_1161
Out[10]:
SHF SST TEMP
month
1 10.912917 18.499498 4.222409
2 9.633392 18.571053 4.224317
3 10.012563 18.560223 4.226013
4 0.360650 18.469493 4.226733
5 -11.141526 18.404369 4.225935
6 -20.389392 18.393965 4.223153
7 -19.051755 18.352014 4.219430
8 -13.911678 18.300126 4.216522
9 3.127374 18.220129 4.215540
10 9.083420 18.150645 4.216613
11 11.101028 18.177402 4.218770
12 12.858325 18.295630 4.220761

Mean annual cycle¶

In [11]:
shf_dict = {'PI control': df_picontrol.SHF, '1dpop.010': df_rest.SHF, '1dpop.015': df_rest2.SHF, 'yr_1161': df_1161.SHF}
sst_dict = {'PI control': df_picontrol.SST, '1dpop.010': df_rest.SST, '1dpop.015': df_rest2.SST, 'yr_1161': df_1161.SST}
temp_dict = {'PI control': df_picontrol.TEMP, '1dpop.010': df_rest.TEMP, '1dpop.015': df_rest2.TEMP, 'yr_1161': df_1161.TEMP}

df_shf = pd.DataFrame(shf_dict)
df_sst = pd.DataFrame(sst_dict)
df_temp = pd.DataFrame(temp_dict)

p_shf = df_shf.hvplot(aspect=1.1, frame_height=350, grid=True, title='SHF', ylabel=units['SHF'])
p_sst = df_sst.hvplot(aspect=1.1, frame_height=350, grid=True, title='SST', ylabel=units['SST'])
p_temp = df_temp.hvplot(aspect=1.1, frame_height=350, grid=True, title='TEMP', ylabel=units['TEMP'])
show(hv.render((p_shf + p_sst + p_temp).opts(shared_axes=False).cols(2)))

Temporal mean maps¶

SHF¶

In [12]:
ds_pi_shf = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.SHF.mean.nc'), decode_times=False)
ds_pi_shf = fixtime(ds_pi_shf)
pi_shf = ds_pi_shf.SHF.squeeze(drop=True)

ds_010_shf = xr.open_dataset(os.path.join(datadir, 'cesm2.B1850.f09_g17.1dpop.010.SHF.annual.nc'), decode_times=False)
ds_010_shf = fixtime(ds_010_shf)
da_010_shf = ds_010_shf.SHF.mean(dim='time')

shf_diff = pi_shf - da_010_shf

proj = ccrs.Robinson(central_longitude=210)

vmax = np.abs([pi_shf.min().item(), pi_shf.max().item()]).max() * 0.8
p1 = pi_shf.hvplot.quadmesh(x='TLONG', y='TLAT', geo=True, coastline=True, projection=proj, project=True,
                            cmap='bwr', clim=(-vmax, vmax),rasterize=True, global_extent=True,
                            clabel=pi_shf.units, title='SHF PI control temporal mean')

vmax = np.abs([da_010_shf.min().item(), da_010_shf.max().item()]).max() * 0.8
p2 = da_010_shf.hvplot.quadmesh(x='TLONG', y='TLAT', geo=True, coastline=True, projection=proj, project=True,
                            cmap='bwr', clim=(-vmax, vmax),rasterize=True, global_extent=True,
                            clabel=ds_010_shf.SHF.units, title='SHF 1dpop.010 temporal mean')

vmax = np.abs([shf_diff.min().item(), shf_diff.max().item()]).max() * 0.8
p3 = shf_diff.hvplot.quadmesh(x='TLONG', y='TLAT', geo=True, coastline=True, projection=proj, project=True,
                            cmap='bwr', clim=(-vmax, vmax),rasterize=True, global_extent=True,
                            clabel=pi_shf.units, title='SHF PI control - 1dpop.010')

(p1 + p2 + p3).cols(2)

WARNING:param.main: Calling the .opts method with options broken down by options group (i.e. separate plot, style and norm groups) is deprecated. Use the .options method converting to the simplified format instead or use hv.opts.apply_groups for backward compatibility.
WARNING:param.main: Calling the .opts method with options broken down by options group (i.e. separate plot, style and norm groups) is deprecated. Use the .options method converting to the simplified format instead or use hv.opts.apply_groups for backward compatibility.
WARNING:param.main: Calling the .opts method with options broken down by options group (i.e. separate plot, style and norm groups) is deprecated. Use the .options method converting to the simplified format instead or use hv.opts.apply_groups for backward compatibility.
Out[12]:

SST¶

In [13]:
ds_pi_sst = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.TEMP.surf.mean.nc'), decode_times=False)
ds_pi_sst = fixtime(ds_pi_sst)
pi_sst = ds_pi_sst.TEMP.squeeze(drop=True)

ds_010_sst = xr.open_dataset(os.path.join(datadir, 'cesm2.B1850.f09_g17.1dpop.010.TEMP.annual.nc'), decode_times=False)
ds_010_sst = fixtime(ds_010_sst)
da_010_sst = ds_010_sst.TEMP.isel(z_t=0).mean(dim='time')

sst_diff = pi_sst - da_010_sst
In [14]:
vmin, vmax = pi_sst.min().item(), pi_sst.max().item()
p1 = pi_sst.hvplot.quadmesh(x='TLONG', y='TLAT', geo=True, coastline=True, projection=proj, project=True,
                            cmap='kbc', clim=(vmin, vmax), rasterize=True, global_extent=True,
                            clabel=pi_sst.units, title='SST PI control temporal mean')

vmin, vmax = da_010_sst.min().item(), da_010_sst.max().item()
p2 = da_010_sst.hvplot.quadmesh(x='TLONG', y='TLAT', geo=True, coastline=True, projection=proj, project=True,
                            cmap='kbc', clim=(vmin, vmax), rasterize=True, global_extent=True,
                            clabel=ds_010_sst.TEMP.units, title='SST 1dpop.010 temporal mean')

vmax = np.abs([sst_diff.min().item(), sst_diff.max().item()]).max()
p3 = sst_diff.hvplot.quadmesh(x='TLONG', y='TLAT', geo=True, coastline=True, projection=proj, project=True,
                            cmap='bwr', clim=(-vmax, vmax), rasterize=True, global_extent=True,
                            clabel=pi_sst.units, title='SST PI control - 1dpop.010')

(p1 + p2 + p3).cols(2)

WARNING:param.main: Calling the .opts method with options broken down by options group (i.e. separate plot, style and norm groups) is deprecated. Use the .options method converting to the simplified format instead or use hv.opts.apply_groups for backward compatibility.
WARNING:param.main: Calling the .opts method with options broken down by options group (i.e. separate plot, style and norm groups) is deprecated. Use the .options method converting to the simplified format instead or use hv.opts.apply_groups for backward compatibility.
WARNING:param.main: Calling the .opts method with options broken down by options group (i.e. separate plot, style and norm groups) is deprecated. Use the .options method converting to the simplified format instead or use hv.opts.apply_groups for backward compatibility.
Out[14]: