Notebook

In [7]:

import numpy as np
import matplotlib.pyplot as plt
from salishsea_tools import evaltools as et, viz_tools
import cmocean as cmo
import datetime as dt
import netCDF4 as nc
import matplotlib.dates as mdates
import matplotlib as mpl
import matplotlib.gridspec as gridspec
from matplotlib.colors import LogNorm
import cmocean
import pandas as pd

%matplotlib inline

definitions for matching model output¶

In [2]:

# path to model files:
PATH= '/results/SalishSea/nowcast-green.201812/'

# start and end dates for analysis:
start_date = dt.datetime(2015,1,1)
end_date = dt.datetime(2016,1,1)

# number of days per model file:
flen=1

# dictionary mapping desired model variables to the file types where they are found
filemap={'nitrate':'ptrc_T','silicon':'ptrc_T','ammonium':'ptrc_T','diatoms':'ptrc_T','ciliates':'ptrc_T',
         'flagellates':'ptrc_T','vosaline':'grid_T','votemper':'grid_T'}

# dictionary mapping model file types to their time resolution in hours (1 is hourly files, 24 is daily)
fdict={'ptrc_T':1,'grid_T':1}

# results format
# -- nowcast: files like 01jan15/SalishSea_1h_20150101_20150101_ptrc_T.nc
# -- long: files like SalishSea_1h_20150206_20150804_ptrc_T_20150427-20150506.nc, all in one directory
namfmt='nowcast'

load PSF data and create csv file¶

# load PSF bottle data (returns pandas dataframe) # NO23 is nitrate+nitrate, Si is Silicate; Chl is chlorophyll; Z is depth (m); dtUTC is datetime in UTC # SA is reference salinity (g/kg), CT is Conservative Temperature (deg C) df1=et.loadPSF(loadCTD=True) df1.drop(columns=['pLat','pLon','tdiffH']) print(len(df1)) df1.head()df1.to_csv('/data/eolson/results/MEOPAR/oldDBs/PSFBotChl.csv',index=False)

load csv file¶

In [18]:

df2=pd.read_csv('/data/eolson/results/MEOPAR/oldDBs/PSFBotChl.csv')
# convert dtUTC back to datetime
df2['dtUTC']=[dt.datetime.strptime(ii,'%Y-%m-%d %H:%M:%S') for ii in df2['dtUTC']]
print(len(df2))
df2.head()

Out[18]:

	Lat	Lon	NO23	PO4	Si	Station	Z	dtUTC	Chl	Depth_m	Phaeo	SA	CT	pLat	pLon	tdiffH
0	49.608333	-124.866667	18.65	1.459	68.43	BS1	2.0	2015-02-18 23:07:00	NaN	NaN	NaN	24.3455	8.986387	49.603803	-124.867129	0.059306
1	49.608333	-124.866667	26.62	2.175	55.73	BS1	20.0	2015-02-18 23:07:00	NaN	NaN	NaN	28.0035	8.810689	49.603803	-124.867129	0.059306
2	49.608333	-124.866667	16.25	1.809	39.14	BS1	2.0	2015-03-21 22:55:00	NaN	NaN	NaN	26.9730	9.675332	49.608083	-124.865769	0.060139
3	49.608333	-124.866667	16.83	2.069	39.11	BS1	20.0	2015-03-21 22:55:00	NaN	NaN	NaN	28.3310	9.113702	49.608083	-124.865769	0.060139
4	49.608333	-124.866667	24.27	0.551	49.90	BS1	20.0	2015-04-05 23:10:00	NaN	NaN	NaN	28.6415	9.296229	49.607844	-124.867736	0.029306

In [28]:

# match model output to observations and return both in a dataframe
# the model variables will have their original names prefixed by mod_, eg mod_vosaline
# the observation file names are unchanged. 
data=et.matchData(data=df2,filemap=filemap, fdict=fdict, mod_start=start_date, mod_end=end_date, 
                  mod_nam_fmt=namfmt, mod_basedir=PATH, mod_flen=flen)

In [29]:

#Chl:N ratio used later in plots: grab from appropriate namelist
mod_chl_N=et.getChlNRatio(basedir=PATH,nam_fmt=namfmt)

In [30]:

# plot matched data locations
fig, ax = plt.subplots(figsize = (6,6))
viz_tools.set_aspect(ax, coords = 'map')
ax.plot(data['Lon'], data['Lat'], 'ro',label='matched obs')
grid = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
viz_tools.plot_coastline(ax, grid, coords = 'map')
ax.set_ylim(48, 50.5)
ax.legend()
ax.set_xlim(-125.7, -122.5);

/data/eolson/results/MEOPAR/tools/SalishSeaTools/salishsea_tools/viz_tools.py:123: UserWarning: No contour levels were found within the data range.
  contour_lines = axes.contour(

In [31]:

# create dictionary of dataframe views by year
datyr=dict()
yy=np.array([ii.year for ii in data['dtUTC']])
yys=np.unique(yy)
for yr in yys:
    datyr[yr]=data.loc[yy==yr]

In [32]:

# plot matched data sampling times
clist=('gold','aqua','plum','c','m','r','g','b','brown','gray')
fig,axL=plt.subplots(1,1,figsize=(12,1))
for ii, yr in enumerate(yys):
    dshift=dt.datetime(yys[0],1,1)-dt.datetime(yr,1,1)
    axL.plot(datyr[yr].dtUTC+dshift,np.zeros(np.shape(datyr[yr].dtUTC))+.1*ii,'.',
             color=clist[ii],markersize=4,label=str(yr))
axL.set_yticks(());
yearsFmt = mdates.DateFormatter('%d %b')
axL.xaxis.set_major_formatter(yearsFmt)
axL.xaxis.set_ticks([dt.datetime(int(yys[0]),1,1), dt.datetime(int(yys[0]),3,1),dt.datetime(int(yys[0]),5,1),dt.datetime(int(yys[0]),7,1),
                     dt.datetime(int(yys[0]),9,1),dt.datetime(int(yys[0]),11,1),dt.datetime(int(yys[0])+1,1,1)])
for tick in axL.get_xticklabels():
        tick.set_rotation(90)
        tick.set_horizontalalignment('center')
axL.set_ylim(-.1,.1*(len(datyr.keys())+1))
axL.set_xlim(dt.datetime(yys[0],1,1),dt.datetime(yys[0],12,31))
axL.legend()
axL.set_frame_on(False)
xmin, xmax = axL.get_xaxis().get_view_interval()
ymin, ymax = axL.get_yaxis().get_view_interval()
axL.add_artist(mpl.lines.Line2D((xmin, xmax), (ymin, ymin), color='black', linewidth=2))

Out[32]:

<matplotlib.lines.Line2D at 0x7f9abe60beb0>

Display stats¶

In [34]:

print('NO3 depth<15:') 
et.printstats(data.loc[data.Z<15,:],'NO23','mod_nitrate')
print('dSi depth<15:')
et.printstats(data.loc[data.Z<15,:],'Si','mod_silicon')

NO3 depth<15:
  N: 391
  bias: -0.7725768276013598
  RMSE: 4.703735413902654
  WSS: 0.9085849267808243
dSi depth<15:
  N: 391
  bias: 6.810807028312716
  RMSE: 16.32044125034323
  WSS: 0.6644405124413911

In [35]:

# chlorophyll calculations
data['l10_obsChl']=np.log10(data['Chl']+0.001)
data['l10_modChl']=np.log10(mod_chl_N*(data['mod_diatoms']+data['mod_ciliates']+data['mod_flagellates'])+0.001)
data['mod_Chl']=mod_chl_N*(data['mod_diatoms']+data['mod_ciliates']+data['mod_flagellates'])

mod vs. obs plots¶

In [37]:

cm1=cmocean.cm.thermal
bounds = np.array([0,5,10,35,100,200,250,300,450])
norm = mpl.colors.BoundaryNorm(boundaries=bounds, ncolors=256)
#pcm = ax[1].pcolormesh(X, Y, Z1, norm=norm, cmap='RdBu_r')
fig,ax=plt.subplots(1,3,figsize=(12,3))
args={'marker':'.','s':8,'norm':norm}
ps=et.varvarScatter(ax[0],data,'NO23','mod_nitrate','Z',cm=cm1,args=args)
ps=et.varvarScatter(ax[1],data,'Si','mod_silicon','Z',cm=cm1,args=args)
ps=et.varvarScatter(ax[2],data,'l10_obsChl','l10_modChl','Z',cm=cm1,args=args)

cb=fig.colorbar(ps,ax=ax[2],boundaries=np.linspace(0,450,46))
cb.set_label('Depth (m)')
for iax in ax:
    iax.set_xlabel('Model')
    iax.set_ylabel('Obs')
ax[0].set_title('Nitrate ($\mu$M)')
ax[1].set_title('dSi ($\mu$M)')
ax[2].set_title('log10[Chl (mg/m3)]')

Out[37]:

Text(0.5, 1.0, 'log10[Chl (mg/m3)]')

In [ ]: