In [1]:
import numpy as np
import xarray as xr
import pandas as pd
from salishsea_tools import viz_tools, places, visualisations
from matplotlib import pyplot as plt, dates
from datetime import datetime, timedelta
from calendar import month_name
from scipy.io import loadmat
from tqdm.notebook import tqdm
from salishsea_tools import nc_tools
from dask.diagnostics import ProgressBar
import cmocean
from scipy.stats import sem
import scipy.stats as stats
%matplotlib inline
In [2]:
plt.rcParams.update({'font.size': 12, 'axes.titlesize': 'medium'})
Load files from monthly averages¶
Plot Seasonal Cycles for Individual Years¶
In [3]:
#years, months, data
monthly_array_PAR_slice = np.zeros([14,12,50,50])
# Load monthly averages
mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc')
slc = {'y': slice(450,500), 'x': slice(250,300)}
e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')]
years, variables = range(2007, 2021), ['PAR']
# PARorary list dict
data = {}
# Permanent aggregate dict
aggregates = {var: {} for var in variables}
monthlydat = {var: {} for var in variables}
# Loop through years
for year in [2007,2008,2009,2010,2011,2012,2016,2017,2018,2019,2020]:
# Initialize lists
for var in variables: data[var] = []
# Load monthly averages
for month in range(1, 13):
datestr = f'{year}{month:02d}'
prefix = f'/data/sallen/results/MEOPAR/v201905r/SalishSea_1m_{datestr}_{datestr}'
# Load grazing variables
with xr.open_dataset(prefix + '_carp_T.nc') as ds:
q = ds.PAR.isel(deptht=0, **slc).values
q2 = q[0,:,:]
monthly_array_PAR_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension
for var in ['PAR']:
data[var].append(ds.PAR.isel(deptht=0, **slc).values)
# Concatenate months
for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)
# Loop through years for wrap files
for year in [2013,2014,2015]:
# Initialize lists
for var in variables: data[var] = []
# Load monthly averages
for month in range(1, 13):
datestr = f'{year}{month:02d}'
prefix = f'/data/sallen/results/MEOPAR/v201905r_wrap/SalishSea_1m_{datestr}_{datestr}'
# Load grazing variables
with xr.open_dataset(prefix + '_carp_T.nc') as ds:
q = ds.PAR.isel(deptht=0, **slc).values
q2 = q[0,:,:]
monthly_array_PAR_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension
for var in ['PAR']:
data[var].append(ds.PAR.isel(deptht=0, **slc).values)
# Concatenate months
for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)
In [4]:
monthly_array_PAR_slice[monthly_array_PAR_slice == 0 ] = np.nan
monthly_array_PAR_slicemean = \
np.nanmean(np.nanmean(monthly_array_PAR_slice, axis = 2),axis = 2)
print(np.shape(monthly_array_PAR_slicemean))
(14, 12)
Plot Seasonal Cycles for Individual Years¶
In [5]:
#plot monthly means for all years
fig, ax = plt.subplots(figsize=(15, 6))
bbox = {'boxstyle': 'round', 'facecolor': 'w', 'alpha': 0.9}
cmap = plt.get_cmap('tab10')
palette = [cmap(0), cmap(0.2), 'k', cmap(0.1), cmap(0.3)]
for i in range(0,6):
ax.plot(np.arange(1,13), monthly_array_PAR_slicemean[i,:],label=2007+i)
ax.set_title('Central SoG PAR Seasonal Cycle',fontsize=18)
ax.legend(frameon=False)
ax.set_ylim(40,375)
ax.set_ylabel('\u03bcmol N')
for i in range(7,14):
ax.plot(np.arange(1,13), monthly_array_PAR_slicemean[i,:],linestyle='--',label=2007+i)
ax.set_title('Central SoG Surface PAR Seasonal Cycle',fontsize=18)
ax.legend(frameon=False,bbox_to_anchor=(1, 1))
ax.set_ylim(0,150)
ax.set_ylabel('m$^{-2}$')
Select 4 warmest and 4 coldest years; leave NPGO "neutral" years out¶
In [6]:
#2008, 2010, 2011, 2012
NPGO_C=(((+monthly_array_PAR_slicemean[1,:]+\
monthly_array_PAR_slicemean[3,:]+\
monthly_array_PAR_slicemean[4,:]+monthly_array_PAR_slicemean[5,:])/4))
In [7]:
#2015, 2018, 2019, 2020
NPGO_W=(((monthly_array_PAR_slicemean[8,:]+\
monthly_array_PAR_slicemean[11,:]+monthly_array_PAR_slicemean[12,:]+\
monthly_array_PAR_slicemean[13,:])/4))
In [8]:
## Plot the coldest and warmest years only; Supp Fig. S3
fig, ax = plt.subplots(figsize=(14, 2))
bbox = {'boxstyle': 'round', 'facecolor': 'w', 'alpha': 0.9}
cmap = plt.get_cmap('tab10')
palette = [cmap(0), cmap(0.2), 'k', cmap(0.1), cmap(0.3)]
xticks=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov',"Dec"]
ax.plot(xticks, monthly_array_PAR_slicemean[1,:],color='b',linestyle='-',label='2008')
ax.plot(xticks, monthly_array_PAR_slicemean[3,:],color='b',linestyle='--',label='2010')
ax.plot(xticks, monthly_array_PAR_slicemean[4,:],color='b',linestyle='-.',label='2011')
ax.plot(xticks, monthly_array_PAR_slicemean[5,:],color='b',linestyle=':',label='2012')
ax.plot(xticks, monthly_array_PAR_slicemean[8,:],color='r',linestyle='-',label='2015')
ax.plot(xticks, monthly_array_PAR_slicemean[11,:],color='r',linestyle='--',label='2018')
ax.plot(xticks, monthly_array_PAR_slicemean[12,:],color='r',linestyle='-.',label='2019')
ax.plot(xticks, monthly_array_PAR_slicemean[13,:],color='r',linestyle=':',label='2020')
ax.set_title('Surface PAR',fontsize=14)
ax.legend(frameon=False,loc='center left', bbox_to_anchor=(1, 0.5),fontsize=10)
ax.set_ylim(0,150)
ax.set_ylabel('uE/m$^{2}$/s')
ax.xaxis.set_tick_params(labelsize=12)
ax.yaxis.set_tick_params(labelsize=12)
ax.set_xticklabels([])
Out[8]:
[Text(0, 0, ''), Text(1, 0, ''), Text(2, 0, ''), Text(3, 0, ''), Text(4, 0, ''), Text(5, 0, ''), Text(6, 0, ''), Text(7, 0, ''), Text(8, 0, ''), Text(9, 0, ''), Text(10, 0, ''), Text(11, 0, '')]
Data for Figure 4: Calculate monthly standard error values for cold and warm years¶
In [9]:
NPGO_W_years=[monthly_array_PAR_slicemean[8,:],monthly_array_PAR_slicemean[11,:],monthly_array_PAR_slicemean[12,:],monthly_array_PAR_slicemean[13,:]]
In [10]:
sem(NPGO_W_years)
Out[10]:
array([1.47123463, 0.57532029, 3.05257932, 2.50705727, 4.39118173,
5.75917672, 2.30639185, 2.69061531, 2.1721261 , 0.98056129,
1.51775549, 0.60773262])
In [11]:
NPGO_W_SEM=[1.47123463, 0.57532029, 3.05257932, 2.50705727, 4.39118173,
5.75917672, 2.30639185, 2.69061531, 2.1721261 , 0.98056129,
1.51775549, 0.60773262]
In [12]:
NPGO_C_years=[monthly_array_PAR_slicemean[1,:],
monthly_array_PAR_slicemean[3,:],
monthly_array_PAR_slicemean[4,:],monthly_array_PAR_slicemean[5,:]]
In [13]:
sem(NPGO_C_years)
Out[13]:
array([0.98927229, 0.80121518, 2.21281148, 2.65062751, 4.36595316,
6.5682177 , 3.54313251, 2.28383506, 2.1619708 , 1.18555839,
0.89343873, 1.10119006])
In [14]:
NPGO_C_SEM=[0.98927229, 0.80121518, 2.21281148, 2.65062751, 4.36595316,
6.5682177 , 3.54313251, 2.28383506, 2.1619708 , 1.18555839,
0.89343873, 1.10119006]
In [15]:
NPGO_C
Out[15]:
array([12.04400661, 23.23528714, 36.51879487, 54.5113959 , 76.02100952,
87.81722026, 99.97068933, 86.76086059, 61.21366767, 31.11335206,
15.38539501, 10.86040996])
In [17]:
## Preliminary Figure 4d
fig, ax = plt.subplots(figsize=(15, 3))
bbox = {'boxstyle': 'round', 'facecolor': 'w', 'alpha': 0.9}
cmap = plt.get_cmap('tab10')
palette = [cmap(0), cmap(0.2), 'k', cmap(0.1), cmap(0.3)]
xticks=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov',"Dec"]
ax.errorbar(xticks, NPGO_C, yerr= NPGO_C_SEM, capsize=3,label='NPGO+ coldest',linewidth=2)
ax.errorbar(xticks, NPGO_W,yerr= NPGO_W_SEM, capsize=3,linestyle='--',label='NPGO- warmest',color='r',linewidth=2)
ax.set_title('Surface PAR',fontsize=18)
ax.legend((),frameon=False)
ax.set_ylim(0,150)
ax.set_ylabel('uE/m$^{2}$/s',fontsize=14)
ax.xaxis.set_tick_params(labelsize=14)
ax.yaxis.set_tick_params(labelsize=14)
ax.set_xticklabels([])
a=2
b=5
c=5
d=8
ax.text(-.4, 130, '(d)', fontsize=15, color='k')
#plt.fill([a, a, b, b], [0, 25, 25, 0], color = 'lightgreen', alpha = 0.1)
#plt.fill([c, c, d, d], [0, 25, 25, 0], color = 'wheat', alpha = 0.2)
#plt.savefig('Fig3d_PAR.png', bbox_inches='tight',dpi=1000,transparent=False)
Out[17]:
Text(-0.4, 130, '(d)')
Data for calculating seasonal mean values for box plots¶
In [18]:
Spring_C=(((monthly_array_PAR_slicemean[(1,3,4,5),2]+
monthly_array_PAR_slicemean[(1,3,4,5),3]+monthly_array_PAR_slicemean[(1,3,4,5),4]))/3) #
In [19]:
Spring_W=(((monthly_array_PAR_slicemean[(8,11,12,13),2]+
monthly_array_PAR_slicemean[(8,11,12,13),3]+monthly_array_PAR_slicemean[(8,11,12,13),4]))/3) #
In [20]:
Summer_C=(((monthly_array_PAR_slicemean[(1,3,4,5),5]+
monthly_array_PAR_slicemean[(1,3,4,5),6]+monthly_array_PAR_slicemean[(1,3,4,5),7]))/3) #
In [21]:
Summer_C
Out[21]:
array([91.57734924, 97.41059098, 94.13498441, 82.94210228])
In [22]:
Summer_W=(((monthly_array_PAR_slicemean[(8,11,12,13),5]+
monthly_array_PAR_slicemean[(8,11,12,13),6]+monthly_array_PAR_slicemean[(8,11,12,13),7]))/3) #
In [23]:
Summer_W
Out[23]:
array([94.68697684, 93.73494871, 92.82362585, 85.96508285])
In [24]:
## Preliminary figure 4d
def color_boxplot(data, color, pos=[0], ax=None):
ax = ax or plt.gca()
bp = ax.boxplot(data, patch_artist=True, showmeans=False, positions=pos,widths=0.4)
for item in ['boxes']:
plt.setp(bp[item], color=color)
for item in ['whiskers', 'fliers', 'medians', 'caps']:
plt.setp(bp[item], color='k')
data1 = [Spring_C]
data2 = [Spring_W]
data3 = [Summer_C]
data4 = [Summer_W]
fig, ax = plt.subplots(figsize=(3,3))
bp1 = color_boxplot(data1, 'royalblue', [1])
bp2 = color_boxplot(data2, 'r', [1.5])
bp3 = color_boxplot(data3, 'royalblue', [2.5])
bp4 = color_boxplot(data4, 'r', [3])
#ax.autoscale()
ax.set(xticks=[1.25,2.75], xticklabels=['Spring','Summer'])
ax.set_ylim(25,125)
ax.set_ylabel('uE/m$^{2}$/s')
#ax.legend([bp1["boxes"], bp2["boxes"], ['A', 'B'], loc='upper right')
plt.show()
In [25]:
Summer_C
Out[25]:
array([91.57734924, 97.41059098, 94.13498441, 82.94210228])
In [26]:
Spring_C.mean()
Out[26]:
55.683733429515605
In [28]:
Spring_W.mean()
Out[28]:
66.93513386529416
t tests for differences between spring and summer of cold and warm years¶
In [29]:
stats.ttest_ind(a=Spring_C, b=Spring_W, equal_var=True)
Out[29]:
Ttest_indResult(statistic=-4.6381287268844655, pvalue=0.003547412512743722)
In [30]:
stats.ttest_ind(a=Summer_C, b=Summer_W, equal_var=True)
Out[30]:
Ttest_indResult(statistic=-0.07787819822628278, pvalue=0.9404571289640673)
In [ ]: