Monthly standard deviations of the regional means

Created by Ivan Lima on Tue Apr 20 2021 12:42:40 -0400

%matplotlib inline
import xarray as xr
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import os, datetime, warnings
from cesm_utils import fixtime
print('Last updated on {}'.format(datetime.datetime.now().ctime()))

Last updated on Wed Apr 21 09:08:28 2021

sns.set_theme()
sns.set_context('paper')
warnings.filterwarnings('ignore')
plt.rcParams['figure.dpi'] = 100

def region_mean(ds, vname, regnum):
    """Compute time series of area weighted means for given region"""
    data = ds[vname].where(ds.REGION_MASK==regnum)
    weights = ds.TAREA.where(ds.REGION_MASK==regnum)
    return ((data * weights).sum(dim=['nlat','nlon']) / weights.sum()).to_series()

def socean_mean(ds, vname):
    """Compute time series of area weighted means for Southern Ocean"""
    data = ds[vname].where(ds.TLAT<=-60)
    weights = ds.TAREA.where((ds.TLAT<=-60) & (ds.REGION_MASK!=0))
    return ((data * weights).sum(dim=['nlat','nlon']) / weights.sum()).to_series()

datadir = '/home/ivan/Data/Postproc/1DPOP_CESM2'
start_yr = 1851 # use "real" year to keep pandas datetime functionality 

ds_gterm = xr.open_dataset(os.path.join(datadir, 'cesm2.B1850.f09_g17.1dpop.002.0005.monthly.nc'), decode_times=False, chunks='auto')
ds_gterm = fixtime(ds_gterm, start_yr)

Labrador Sea

stats_labsea = {}
varlist = ['TEMP','SALT','UVEL','VVEL']
for vname in varlist:
    ds = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.{}.upper_200m.nc'.format(vname)),
                         decode_times=False, chunks='auto')
    ds = fixtime(ds, start_yr-449)
    ds['TAREA'] = ds_gterm['TAREA']
    ds['REGION_MASK'] = ds_gterm['REGION_MASK']
    df = pd.DataFrame({'month':ds.time.dt.month, 'year':ds.time.dt.year, vname:region_mean(ds,vname,8)})
    df_stats = pd.DataFrame({'mean': df.groupby('month')[vname].mean(), 'std': df.groupby('month')[vname].std()})
    stats_labsea[vname] = df_stats

df_labsea = pd.concat([stats_labsea['TEMP'], stats_labsea['SALT'], stats_labsea['UVEL'], stats_labsea['VVEL']],
                      keys=varlist, axis=1)

df_std = pd.DataFrame({vname: stats_labsea[vname]['std'] for vname in ['TEMP','SALT','UVEL','VVEL']})
axs = df_std.plot(subplots=True, figsize=(8,6), layout=(2,2), title='STD')

df_labsea

	TEMP		SALT		UVEL		VVEL
	mean	std	mean	std	mean	std	mean	std
month
1	1.983920	0.304538	33.925045	0.078882	0.703006	0.120689	-2.132192	0.223754
2	1.967708	0.290566	34.011390	0.078924	0.811217	0.139318	-2.120554	0.249357
3	1.986412	0.263498	34.065488	0.073579	0.817975	0.122949	-1.950958	0.215034
4	2.034051	0.254256	34.082101	0.069899	0.742525	0.102929	-1.738567	0.170556
5	2.187921	0.255091	34.054258	0.066570	0.671951	0.084705	-1.573552	0.134365
6	2.489094	0.259417	33.984632	0.064439	0.647401	0.076301	-1.591031	0.119944
7	2.858791	0.264076	33.897623	0.063009	0.665884	0.071384	-1.641827	0.114551
8	3.199733	0.267308	33.836706	0.061682	0.680718	0.069572	-1.701105	0.122657
9	3.305709	0.262877	33.808679	0.060531	0.681730	0.073449	-1.804396	0.138790
10	3.110800	0.257337	33.804798	0.059831	0.661153	0.077866	-1.869563	0.169176
11	2.732701	0.253671	33.821815	0.059861	0.634859	0.086675	-1.942890	0.201490
12	2.273833	0.265502	33.854710	0.063514	0.621267	0.093924	-2.048271	0.209548

Arctic

stats_arctic = {}
for vname in varlist:
    ds = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.{}.upper_200m.nc'.format(vname)),
                         decode_times=False, chunks='auto')
    ds = fixtime(ds, start_yr-449)
    ds['TAREA'] = ds_gterm['TAREA']
    ds['REGION_MASK'] = ds_gterm['REGION_MASK']
    df = pd.DataFrame({'month':ds.time.dt.month, 'year':ds.time.dt.year, vname:region_mean(ds,vname,10)})
    df_stats = pd.DataFrame({'mean': df.groupby('month')[vname].mean(), 'std': df.groupby('month')[vname].std()})
    stats_arctic[vname] = df_stats

df_arctic = pd.concat([stats_arctic['TEMP'], stats_arctic['SALT'], stats_arctic['UVEL'], stats_arctic['VVEL']],
                      keys=varlist, axis=1)

df_std = pd.DataFrame({vname: stats_arctic[vname]['std'] for vname in ['TEMP','SALT','UVEL','VVEL']})
axs = df_std.plot(subplots=True, figsize=(8,6), layout=(2,2), title='STD')

df_arctic

	TEMP		SALT		UVEL		VVEL
	mean	std	mean	std	mean	std	mean	std
month
1	-1.110802	0.037254	32.438846	0.064759	0.099071	0.164736	0.352153	0.071961
2	-1.141129	0.032685	32.533054	0.064267	0.098720	0.165973	0.376181	0.069712
3	-1.162925	0.030919	32.614050	0.063905	0.083311	0.153204	0.406938	0.056535
4	-1.173221	0.031133	32.675568	0.063725	0.007304	0.149489	0.432869	0.048573
5	-1.142362	0.037308	32.681625	0.063562	-0.159665	0.127622	0.464744	0.044131
6	-1.003857	0.059281	32.560425	0.065171	-0.225588	0.137554	0.462692	0.044439
7	-0.766733	0.099968	32.347998	0.066670	-0.145637	0.151565	0.428501	0.044019
8	-0.527092	0.135728	32.183151	0.066629	-0.091351	0.176412	0.406737	0.052305
9	-0.464945	0.143437	32.117046	0.066727	0.014117	0.216542	0.365333	0.067767
10	-0.673018	0.115512	32.129476	0.066021	0.052415	0.232807	0.344996	0.075070
11	-0.929636	0.072583	32.215740	0.065348	0.050964	0.218556	0.332399	0.085132
12	-1.058506	0.048934	32.328450	0.065263	0.074568	0.182161	0.334063	0.080590

Southern Ocean

stats_socean = {}
for vname in varlist:
    ds = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.{}.upper_200m.nc'.format(vname)),
                         decode_times=False, chunks='auto')
    ds = fixtime(ds, start_yr-449)
    ds['TAREA'] = ds_gterm['TAREA']
    ds['REGION_MASK'] = ds_gterm['REGION_MASK']
    df = pd.DataFrame({'month':ds.time.dt.month, 'year':ds.time.dt.year, vname:socean_mean(ds,vname)})
    df_stats = pd.DataFrame({'mean': df.groupby('month')[vname].mean(), 'std': df.groupby('month')[vname].std()})
    stats_socean[vname] = df_stats

df_socean = pd.concat([stats_socean['TEMP'], stats_socean['SALT'], stats_socean['UVEL'], stats_socean['VVEL']],
                      keys=varlist, axis=1)

df_std = pd.DataFrame({vname: stats_socean[vname]['std'] for vname in ['TEMP','SALT','UVEL','VVEL']})
axs = df_std.plot(subplots=True, figsize=(8,6), layout=(2,2), title='STD')

df_socean

	TEMP		SALT		UVEL		VVEL
	mean	std	mean	std	mean	std	mean	std
month
1	-0.152994	0.071203	33.865804	0.017691	0.870104	0.119599	-0.034023	0.072204
2	0.055251	0.077155	33.840427	0.018070	0.826463	0.136213	0.023124	0.082879
3	0.095116	0.078524	33.840053	0.018390	0.805919	0.150832	0.116732	0.088003
4	-0.021623	0.075793	33.858584	0.018677	0.729375	0.169668	0.137611	0.095067
5	-0.188140	0.070075	33.884705	0.018831	0.681011	0.179336	0.111559	0.098773
6	-0.337619	0.063946	33.908998	0.018847	0.660044	0.184264	0.086959	0.103727
7	-0.465940	0.058742	33.929652	0.018787	0.663166	0.173973	0.099172	0.105707
8	-0.575693	0.053405	33.947444	0.018621	0.656694	0.175447	0.114246	0.110391
9	-0.651218	0.050351	33.961274	0.018399	0.668899	0.170968	0.146276	0.110242
10	-0.667294	0.050888	33.968111	0.018123	0.777465	0.164039	0.174405	0.105914
11	-0.587276	0.054997	33.958711	0.017945	0.912381	0.156843	0.131430	0.099496
12	-0.405327	0.062110	33.919118	0.017862	0.920342	0.134059	0.012986	0.084957

Save stats to CSV file

for vname in varlist:
    stats_arctic[vname].to_csv('stats/arctic_{}.csv'.format(vname))
    stats_labsea[vname].to_csv('stats/labsea_{}.csv'.format(vname))
    stats_socean[vname].to_csv('stats/socean_{}.csv'.format(vname))