<jupyter_start><jupyter_text>Notebook based on Elise's version<jupyter_code>import datetime as dt
import importlib
import matplotlib as mpl
from matplotlib.colors import LogNorm
import matplotlib.dates as mdates
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from salishsea_tools import evaltools as et
from salishsea_tools import viz_tools

%matplotlib inline
importlib.reload(et)
# path to model files:
PATH = '/results2/SalishSea/nowcast-green.201905/'

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

# number of days per model file:
flen = 1

# dictionary mapping desired model variables to the file types where they are found
filemap = {'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 DFO CTD data (returns pandas dataframe)
# SA is Absolute (actually reference) Salinity, and CT is Conservative Temperature
df1 = et.loadDFOCTD(datelims=(start_date, end_date))
print(len(df1))
df1.head()
from sqlalchemy import create_engine, case
from sqlalchemy.orm import create_session
from sqlalchemy.ext.automap import automap_base
from sqlalchemy.sql import and_, or_, not_, func

basedir='/ocean/shared/SalishSeaCastData/DFO/CTD/'
dbname='DFO_CTD.sqlite'
engine = create_engine('sqlite:///' + basedir + dbname, echo = False)
Base = automap_base()

Base.prepare(engine, reflect=True)
# mapped classes have been created
# existing tables:
StationTBL=Base.classes.StationTBL
ObsTBL=Base.classes.ObsTBL
CalcsTBL=Base.classes.CalcsTBL
session = create_session(bind = engine, autocommit = False, autoflush = True)
SA=case([(CalcsTBL.Salinity_T0_C0_SA!=None, CalcsTBL.Salinity_T0_C0_SA)], else_=
         case([(CalcsTBL.Salinity_T1_C1_SA!=None, CalcsTBL.Salinity_T1_C1_SA)], else_=
         case([(CalcsTBL.Salinity_SA!=None, CalcsTBL.Salinity_SA)], else_= None)))
CT=case([(CalcsTBL.Temperature_Primary_CT!=None, CalcsTBL.Temperature_Primary_CT)], else_=
         case([(CalcsTBL.Temperature_Secondary_CT!=None, CalcsTBL.Temperature_Secondary_CT)], else_=CalcsTBL.Temperature_CT))
ZD=case([(ObsTBL.Depth!=None,ObsTBL.Depth)], else_= CalcsTBL.Z)
FL=case([(ObsTBL.Fluorescence_URU_Seapoint!=None,ObsTBL.Fluorescence_URU_Seapoint)], else_= ObsTBL.Fluorescence_URU_Wetlabs)

datelims = (start_date, end_date)
start_y=datelims[0].year
start_m=datelims[0].month
start_d=datelims[0].day
end_y=datelims[1].year
end_m=datelims[1].month
end_d=datelims[1].day
qry=session.query(StationTBL.StartYear.label('Year'),StationTBL.StartMonth.label('Month'),
              StationTBL.StartDay.label('Day'),StationTBL.StartHour.label('Hour'),
              StationTBL.Lat,StationTBL.Lon,ZD.label('Z'),SA.label('SA'),CT.label('CT'),FL.label('Fluor')).\
        select_from(StationTBL).join(ObsTBL,ObsTBL.StationTBLID==StationTBL.ID).\
        join(CalcsTBL,CalcsTBL.ObsTBLID==ObsTBL.ID).filter(and_(or_(StationTBL.StartYear>start_y,
                                                                 and_(StationTBL.StartYear==start_y, StationTBL.StartMonth>start_m),
                                                                 and_(StationTBL.StartYear==start_y, StationTBL.StartMonth==start_m, StationTBL.StartDay>=start_d)),
                                                                or_(StationTBL.StartYear<end_y,
                                                                 and_(StationTBL.StartYear==end_y,StationTBL.StartMonth<end_m),
                                                                 and_(StationTBL.StartYear==end_y,StationTBL.StartMonth==end_m, StationTBL.StartDay<end_d)),
                                                            StationTBL.Lat>47-3/2.5*(StationTBL.Lon+123.5),
                                                            StationTBL.Lat<47-3/2.5*(StationTBL.Lon+121),
                                                            StationTBL.Include==True,ObsTBL.Include==True,CalcsTBL.Include==True))
df1=pd.DataFrame(qry.all())

df1
# 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=df1, filemap=filemap, fdict=fdict, mod_start=start_date, mod_end=end_date, 
                  mod_nam_fmt=namfmt, mod_basedir=PATH, mod_flen=flen)
data
data['year'] = data['Year'].astype(int)
data['month'] = data['Month'].astype(int)
data['day'] = data['Day'].astype(int)
pd.to_datetime(data[['year', 'month', 'day']])
year = data['Year'].astype(int)
year
# 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')
bathy = '/home/sallen/MEOPAR/grid/bathymetry_201702.nc'
viz_tools.plot_land_mask(ax, bathy, coords='map', color='burlywood')
ax.set_ylim(48, 50.5)
ax.legend()
ax.set_xlim(-125.7, -122.5);
# create dictionary of dataframe views by year
datyr = {}
yy = np.array([ii.year for ii in data['dtUTC']])
yys = np.unique(yy)
for yr in np.unique(yy):
    datyr[yr] = data.loc[yy]
# 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(datyr[yr].dtUTC.shape) + 0.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(-0.1, 0.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));<jupyter_output><empty_output><jupyter_text>Display stats<jupyter_code>print('CT:')
et.printstats(data.loc[:], 'CT', 'mod_votemper')
print('SA:')
et.printstats(data.loc[:], 'SA', 'mod_vosaline')<jupyter_output>CT:
  N: 90946
  bias: -0.08334274459062385
  RMSE: 0.38189412450767063
  WSS: 0.9638351717153325
SA:
  N: 90946
  bias: 0.15792155556050957
  RMSE: 0.4877281764689603
  WSS: 0.9690926188489902<jupyter_text>mod vs. obs plots<jupyter_code>fig = plt.figure(figsize = (7.5, 3.25))
gs1 = gridspec.GridSpec(1, 4, left=.08, right=0.91, bottom=0.1, top=0.97,
                      wspace=0.2, width_ratios=[1, 0.1, 1, 0.1])
axT = fig.add_subplot(gs1[0, 0])
axS = fig.add_subplot(gs1[0, 2])
gscb = gridspec.GridSpecFromSubplotSpec(7, 1, subplot_spec=gs1[:, 3])
axcb = fig.add_subplot(gscb[1:-1])

axS.plot((0,36), (0,36), 'k-', alpha=0.2)
axT.plot((0,25), (0,25), 'k-', alpha=0.2)

iiT = (~np.isnan(data['CT'])) & (~np.isnan(data['mod_votemper']))
iiS = (~np.isnan(data['SA'])) & (~np.isnan(data['mod_vosaline']))
counts, xedges, yedges, m1 = axT.hist2d(data.loc[iiT,['CT']].values.flatten(),
                                      data.loc[iiT,['mod_votemper']].values.flatten(), bins=25*3, norm=LogNorm())
counts, xedges, yedges, m2 = axS.hist2d(data.loc[iiS,['SA']].values.flatten(),
                                      data.loc[iiS,['mod_vosaline']].values.flatten(), bins=36*3, norm=LogNorm())

cb = fig.colorbar(m1, cax=axcb)
cb.set_label('Count')

Smax = 35
ntick = np.arange(0, Smax, 10)
ntickl = [str(i) for i in ntick]
axS.set_xlim((0, Smax))
axS.set_ylim((0, Smax))
axS.set_xticks(ntick)
axS.set_xticklabels(ntickl)
axS.set_yticks(ntick)
axS.set_yticklabels(ntickl)
    
Tmax = 22
Tmin = 4
stick=np.arange(Tmin, Tmax, 10)
stickl=[str(i) for i in stick]
axT.set_xlim((Tmin, Tmax))
axT.set_ylim((Tmin, Tmax))
axT.set_xticks(stick)
axT.set_xticklabels(stickl)
axT.set_yticks(stick)
axT.set_yticklabels(stickl)
    
for ax in (axT,axS):
    ax.set_aspect(1, adjustable='box')
    
axT.set_ylabel('Modeled', fontsize=12)
axS.set_ylabel('Modeled', fontsize=12)
axT.set_xlabel('Observed', fontsize=12)
axS.set_xlabel('Observed', fontsize=12)

axS.set_title('S$_A$ (g/kg)', fontsize=12)
axT.set_title('$\Theta$ ($^{\circ}$C)', fontsize=12);
plt.savefig('firsthalf2016_201905.png')
## where is that weird point from?
fig, ax = plt.subplots(figsize = (5,5))
viz_tools.set_aspect(ax, coords = 'map')
ax.plot(data.loc[((data.mod_votemper-data.CT)<-3),['Lon']], data.loc[((data.mod_votemper-data.CT)<-3),['Lat']], 'ro',label='high 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.loc[(data.mod_votemper-data.CT)<-3]<jupyter_output>/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(<jupyter_text>Alternatively, color by depth (Z), or other variable:<jupyter_code>fig, ax = plt.subplots(figsize = (8,8))
ps=et.varvarPlot(ax,data,'CT','mod_votemper','Z',(15,50,100,200),'z','m',('darkorange','gold','mediumseagreen','darkturquoise','navy'))
ax.legend(handles=ps)
ax.set_xlabel('Obs')
ax.set_ylabel('Model')
ax.set_title('$\Theta$ ($^{\circ}$C)')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16149870801033592, 'fraction_numerical': 0.044126416219439475, 'mean_word_length': 4.26032409827496, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '178824', 'license': 'Apache-2.0', 'n_tokens_mistral': 4211, 'n_tokens_neox': 3890, 'n_words': 661, 'path': 'notebooks/PhysTuning/evalTSExample_DFOCTD.ipynb', 'repo_name': 'SalishSeaCast/analysis-susan', 'size': '10062'}