DFO Nutrient Comparison¶

In [1]:

import sqlalchemy
from sqlalchemy import (create_engine, Column, String, Integer, Float, MetaData, 
                        Table, type_coerce, ForeignKey, case)
from sqlalchemy.orm import mapper, create_session, relationship, aliased, Session
from sqlalchemy.ext.declarative import declarative_base
from sqlalchemy import case
import numpy as np
from sqlalchemy.ext.automap import automap_base
import matplotlib.pyplot as plt
import sqlalchemy.types as types
from sqlalchemy.sql import and_, or_, not_, func
from sqlalchemy.sql import select
import os
from os.path import isfile
import pandas as pd
import netCDF4 as nc
import datetime as dt
from salishsea_tools import evaltools as et, viz_tools
import datetime
import glob
import gsw
import matplotlib.gridspec as gridspec
import matplotlib as mpl
import matplotlib.dates as mdates

mpl.rc('xtick', labelsize=16)
mpl.rc('ytick', labelsize=16)
mpl.rc('legend', fontsize=12)
mpl.rc('axes', titlesize=16)
mpl.rc('axes', labelsize=16)
mpl.rc('figure', titlesize=16)
mpl.rc('font', size=16)

%matplotlib inline

In [ ]:

PATH= '/results2/SalishSea/hindcast/'
start_date = datetime.datetime(2015,1,1)
end_date = datetime.datetime(2018,1,1)
flen=1
namfmt='nowcast'
#varmap={'N':'nitrate','Si':'silicon','Ammonium':'ammonium'}
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'}
#gridmap={'nitrate':'tmask','silicon':'tmask','ammonium':'tmask'}
fdict={'ptrc_T':1,'grid_T':1}

df1=et.loadDFO()
df1.head()

Out[ ]:

	Year	Month	Day	Hour	Lat	Lon	Pressure	Depth	Ammonium	Ammonium_units	Chlorophyll_Extracted	Chlorophyll_Extracted_units	N	Si	Silicate_units	AbsSal	T	T_units	Z	dtUTC
0	2015.0	2.0	11.0	11.068611	48.300833	-124.000333	1.9	NaN	None	None	NaN	mg/m^3	15.31	32.14	umol/L	29.227507	9.7647	'deg_C_(ITS90)'	1.883998	2015-02-11 11:04:07
1	2015.0	2.0	11.0	11.068611	48.300833	-124.000333	6.6	NaN	None	None	2.57	mg/m^3	17.13	33.90	umol/L	29.484341	9.6880	'deg_C_(ITS90)'	6.544340	2015-02-11 11:04:07
2	2015.0	2.0	11.0	11.068611	48.300833	-124.000333	6.7	NaN	None	None	NaN	mg/m^3	NaN	NaN	umol/L	29.484839	9.6828	'deg_C_(ITS90)'	6.643495	2015-02-11 11:04:07
3	2015.0	2.0	11.0	11.068611	48.300833	-124.000333	11.0	NaN	None	None	NaN	mg/m^3	NaN	NaN	umol/L	30.144549	9.3646	'deg_C_(ITS90)'	10.907117	2015-02-11 11:04:07
4	2015.0	2.0	11.0	11.068611	48.300833	-124.000333	11.0	NaN	None	None	NaN	mg/m^3	20.62	37.65	umol/L	30.157913	9.3586	'deg_C_(ITS90)'	10.907117	2015-02-11 11:04:07

In [ ]:

data=et.matchData(df1,filemap, fdict, start_date, end_date, namfmt, PATH, flen)

In [ ]:

fig, ax = plt.subplots(figsize = (6,6))
viz_tools.set_aspect(ax, coords = 'map')
ax.plot(data['Lon'], data['Lat'], 'ro',label='data')
ax.plot(data.loc[(data.Lon < -123.5) & (data.Lat < 48.6),['Lon']], 
        data.loc[(data.Lon < -123.5) & (data.Lat < 48.6),['Lat']], 
        'bo', label = 'Juan de Fuca')

ax.plot(data.loc[data.Si>75,['Lon']],data.loc[data.Si>75,['Lat']],'*',color='y',label='high Si')
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);

2015¶

In [ ]:

# 2015
data2015=data.loc[(data.dtUTC>=dt.datetime(2015,1,1))&(data.dtUTC<dt.datetime(2016,1,1))]
idata=data2015
yy='2015'
fig = plt.figure(figsize = (12,4))
gs1 = gridspec.GridSpec(2,3,left=.08,right=.93,hspace=.8,wspace=.7,bottom=.08,top=.9,height_ratios=[5.8,.5],width_ratios=[4,4,3])
ax1 = fig.add_subplot(gs1[0,0])
ax2 = fig.add_subplot(gs1[0,1])
ax3 = fig.add_subplot(gs1[0,2])
axb = fig.add_subplot(gs1[1, :])

ps=et.varvarPlot(ax1,idata,'N','mod_nitrate','Z',(15,22),'z','m',('mediumseagreen','darkturquoise','navy'))
for el in ps:
    el.set_alpha(.3)
#ax1.legend(handles=ps)
ax1.set_xlabel('Observed')
ax1.set_ylabel('Modelled')
ax1.set_title(yy+' NO$_3$ ($\mu$M)')
ax1.set_xlim((0,45))
ax1.set_ylim((0,45))
ntick=np.arange(0,56,15)
ntickl=[str(i) for i in ntick]
ax1.set_xticks(ntick)
ax1.set_xticklabels(ntickl)
ax1.set_yticks(ntick)
ax1.set_yticklabels(ntickl)

ps=et.varvarPlot(ax2,idata,'Si','mod_silicon','Z',(15,22),'z','m',('mediumseagreen','darkturquoise','navy'))
for el in ps:
    el.set_alpha(.3)
ax2.legend(handles=ps,bbox_to_anchor=[.65,.45])
ax2.set_xlabel('Observed')
ax2.set_ylabel('Modelled')
ax2.set_title(yy+' dSi ($\mu$M)')
ax2.set_xlim((0,80))
ax2.set_ylim((0,80))
ntick=np.arange(0,81,20)
ntickl=[str(i) for i in ntick]
ax2.set_xticks(ntick)
ax2.set_xticklabels(ntickl)
ax2.set_yticks(ntick)
ax2.set_yticklabels(ntickl)

viz_tools.set_aspect(ax3, coords = 'map')
datanut=idata.loc[(idata.N>=0)|(idata.Si>=0)]
ax3.plot(datanut['Lon'], datanut['Lat'], '.',color='darkorange',label='data',markersize=6)
grid = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
viz_tools.plot_coastline(ax3, grid, coords = 'map')
ax3.set_ylim(48, 50.7)
#ax3.legend()
ax3.set_xlim(-125.8, -122.4);
ax3.axis('off')
test=ax3.get_position()
ax3.set_position(mpl.transforms.Bbox(points=[[.7,.2],[.97,.99]]))

axb.plot(datanut.dtUTC,np.ones(np.shape(datanut.dtUTC)),'.',color='darkorange',markersize=6)
axb.set_yticks(());
yearsFmt = mdates.DateFormatter('%d %b %Y')
axb.xaxis.set_major_formatter(yearsFmt)
axb.xaxis.set_ticks([dt.datetime(int(yy),1,1), dt.datetime(int(yy),4,1),dt.datetime(int(yy),7,1),dt.datetime(int(yy),10,1),dt.datetime(int(yy)+1,1,1)])
#axb.autofmt_xdate(bottom=0.2, rotation=15, ha='right')
#axb.set_title('Sampling Dates')
#plt.setp(plt.xticks()[1], rotation=10, ha='right')
axb.set_title('Sample Dates')
#fig.savefig('/home/eolson/pyCode/notebooks/figs/EvalDFO2015.png',dpi=200,transparent=True)

2016¶

In [ ]:

# 2016
data2016=data.loc[(data.dtUTC>=dt.datetime(2016,1,1))&(data.dtUTC<dt.datetime(2017,1,1))]
idata=data2016
yy='2016'
fig = plt.figure(figsize = (12,4))
gs1 = gridspec.GridSpec(2,3,left=.08,right=.93,hspace=.8,wspace=.7,bottom=.08,top=.9,height_ratios=[5.8,.5],width_ratios=[4,4,3])
ax1 = fig.add_subplot(gs1[0,0])
ax2 = fig.add_subplot(gs1[0,1])
ax3 = fig.add_subplot(gs1[0,2])
axb = fig.add_subplot(gs1[1, :])

ps=et.varvarPlot(ax1,idata,'N','mod_nitrate','Z',(15,22),'z','m',('mediumseagreen','darkturquoise','navy'))
for el in ps:
    el.set_alpha(.3)
#ax1.legend(handles=ps)
ax1.set_xlabel('Observed')
ax1.set_ylabel('Modelled')
ax1.set_title(yy+' NO$_3$ ($\mu$M)')
ax1.set_xlim((0,45))
ax1.set_ylim((0,45))
ntick=np.arange(0,56,15)
ntickl=[str(i) for i in ntick]
ax1.set_xticks(ntick)
ax1.set_xticklabels(ntickl)
ax1.set_yticks(ntick)
ax1.set_yticklabels(ntickl)

ps=et.varvarPlot(ax2,idata,'Si','mod_silicon','Z',(15,22),'z','m',('mediumseagreen','darkturquoise','navy'))
for el in ps:
    el.set_alpha(.3)
#ax2.legend(handles=ps,bbox_to_anchor=[.65,.4])
ax2.set_xlabel('Observed')
ax2.set_ylabel('Modelled')
ax2.set_title(yy+' dSi ($\mu$M)')
ax2.set_xlim((0,80))
ax2.set_ylim((0,80))
ntick=np.arange(0,81,20)
ntickl=[str(i) for i in ntick]
ax2.set_xticks(ntick)
ax2.set_xticklabels(ntickl)
ax2.set_yticks(ntick)
ax2.set_yticklabels(ntickl)

viz_tools.set_aspect(ax3, coords = 'map')
datanut=idata.loc[(idata.N>=0)|(idata.Si>=0)]
ax3.plot(datanut['Lon'], datanut['Lat'], '.',color='darkorange',label='data',markersize=6)
grid = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
viz_tools.plot_coastline(ax3, grid, coords = 'map')
ax3.set_ylim(48, 50.7)
#ax3.legend()
ax3.set_xlim(-125.8, -122.4);
ax3.axis('off')
test=ax3.get_position()
ax3.set_position(mpl.transforms.Bbox(points=[[.7,.2],[.97,.99]]))

axb.plot(datanut.dtUTC,np.ones(np.shape(datanut.dtUTC)),'.',color='darkorange',markersize=6)
axb.set_yticks(());
yearsFmt = mdates.DateFormatter('%d %b %Y')
axb.xaxis.set_major_formatter(yearsFmt)
axb.xaxis.set_ticks([dt.datetime(int(yy),1,1), dt.datetime(int(yy),4,1),dt.datetime(int(yy),7,1),dt.datetime(int(yy),10,1),dt.datetime(int(yy)+1,1,1)])
#axb.autofmt_xdate(bottom=0.2, rotation=15, ha='right')
#axb.set_title('Sampling Dates')
#plt.setp(plt.xticks()[1], rotation=10, ha='right')
axb.set_title('Sample Dates')
#fig.savefig('/home/eolson/pyCode/notebooks/figs/EvalDFO2016.png',dpi=200,transparent=True)

2017¶

In [ ]:

# 2017
data2016=data.loc[(data.dtUTC>=dt.datetime(2016,1,1))&(data.dtUTC<dt.datetime(2017,1,1))]
idata=data2016
yy='2017'
fig = plt.figure(figsize = (12,4))
gs1 = gridspec.GridSpec(2,3,left=.08,right=.93,hspace=.8,wspace=.7,bottom=.08,top=.9,height_ratios=[5.8,.5],width_ratios=[4,4,3])
ax1 = fig.add_subplot(gs1[0,0])
ax2 = fig.add_subplot(gs1[0,1])
ax3 = fig.add_subplot(gs1[0,2])
axb = fig.add_subplot(gs1[1, :])

ps=et.varvarPlot(ax1,idata,'N','mod_nitrate','Z',(15,22),'z','m',('mediumseagreen','darkturquoise','navy'))
for el in ps:
    el.set_alpha(.3)
#ax1.legend(handles=ps)
ax1.set_xlabel('Observed')
ax1.set_ylabel('Modelled')
ax1.set_title(yy+' NO$_3$ ($\mu$M)')
ax1.set_xlim((0,45))
ax1.set_ylim((0,45))
ntick=np.arange(0,56,15)
ntickl=[str(i) for i in ntick]
ax1.set_xticks(ntick)
ax1.set_xticklabels(ntickl)
ax1.set_yticks(ntick)
ax1.set_yticklabels(ntickl)

ps=et.varvarPlot(ax2,idata,'Si','mod_silicon','Z',(15,22),'z','m',('mediumseagreen','darkturquoise','navy'))
for el in ps:
    el.set_alpha(.3)
#ax2.legend(handles=ps,bbox_to_anchor=[.65,.4])
ax2.set_xlabel('Observed')
ax2.set_ylabel('Modelled')
ax2.set_title(yy+' dSi ($\mu$M)')
ax2.set_xlim((0,80))
ax2.set_ylim((0,80))
ntick=np.arange(0,81,20)
ntickl=[str(i) for i in ntick]
ax2.set_xticks(ntick)
ax2.set_xticklabels(ntickl)
ax2.set_yticks(ntick)
ax2.set_yticklabels(ntickl)

viz_tools.set_aspect(ax3, coords = 'map')
datanut=idata.loc[(idata.N>=0)|(idata.Si>=0)]
ax3.plot(datanut['Lon'], datanut['Lat'], '.',color='darkorange',label='data',markersize=6)
grid = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
viz_tools.plot_coastline(ax3, grid, coords = 'map')
ax3.set_ylim(48, 50.7)
#ax3.legend()
ax3.set_xlim(-125.8, -122.4);
ax3.axis('off')
test=ax3.get_position()
ax3.set_position(mpl.transforms.Bbox(points=[[.7,.2],[.97,.99]]))

axb.plot(datanut.dtUTC,np.ones(np.shape(datanut.dtUTC)),'.',color='darkorange',markersize=6)
axb.set_yticks(());
yearsFmt = mdates.DateFormatter('%d %b %Y')
axb.xaxis.set_major_formatter(yearsFmt)
axb.xaxis.set_ticks([dt.datetime(int(yy),1,1), dt.datetime(int(yy),4,1),dt.datetime(int(yy),7,1),dt.datetime(int(yy),10,1),dt.datetime(int(yy)+1,1,1)])
#axb.autofmt_xdate(bottom=0.2, rotation=15, ha='right')
#axb.set_title('Sampling Dates')
#plt.setp(plt.xticks()[1], rotation=10, ha='right')
axb.set_title('Sample Dates')
#fig.savefig('/home/eolson/pyCode/notebooks/figs/EvalDFO2016.png',dpi=200,transparent=True)

In [ ]:

Chl¶

In [ ]:

data['l10_obsChl']=np.log10(data['Chlorophyll_Extracted']+0.01)
data['l10_modChl']=np.log10(2*(data['mod_diatoms']+data['mod_ciliates']+data['mod_flagellates'])+0.01)
data['mod_Chl']=2*(data['mod_diatoms']+data['mod_ciliates']+data['mod_flagellates'])

In [ ]:

fig = plt.figure(figsize = (12,4.5))
gs1 = gridspec.GridSpec(1,3,left=.1,right=.96,wspace=.5,bottom=.15,top=.85,width_ratios=[4,4,1])
ax1 = fig.add_subplot(gs1[0])
ax2 = fig.add_subplot(gs1[1])
chtick=np.arange(-2,2.1,1)
chtickl=[str(i) for i in chtick]

data2015=data.loc[(data.dtUTC>=dt.datetime(2015,1,1))&(data.dtUTC<dt.datetime(2016,1,1))]
idata=data2015
yy='2015'
ax1.plot(chtick,chtick,'k-')
ps=et.varvarPlot(ax1,idata,'l10_obsChl','l10_modChl','Z',(5,10,15,20,25),'z','m',('crimson','darkorange','lime','forestgreen','darkturquoise','navy'))
for el in ps:
    el.set_alpha(.5)
#ax1.legend(handles=ps)
ax1.set_xlabel('Observed')
ax1.set_ylabel('Modelled')
ax1.set_title(yy+'\nlog10[Chl ($\mu$g/L)+0.01]')
ax1.set_xticks(chtick)
ax1.set_xticklabels(chtickl)
ax1.set_yticks(chtick)
ax1.set_yticklabels(chtickl)
ax1.set_xlim(-2,1.8)
ax1.set_ylim(-2,1.8)


data2016=data.loc[(data.dtUTC>=dt.datetime(2016,1,1))&(data.dtUTC<dt.datetime(2017,1,1))]
idata=data2016
yy='2016'
ax2.plot(chtick,chtick,'k-')
ax2.plot(np.arange(-.6,1.6,.1),np.arange(-.6,1.6,.1),'k-')
ps=et.varvarPlot(ax2,idata,'l10_obsChl','l10_modChl','Z',(5,10,15,20,25),'z','m',('crimson','darkorange','lime','forestgreen','darkturquoise','navy'))
for el in ps:
    el.set_alpha(.5)
ax2.legend(handles=ps,bbox_to_anchor=[1.08,.8])
ax2.set_xlabel('Observed')
ax2.set_ylabel('Modelled')
ax2.set_title(yy+'\nlog10[Chl ($\mu$g/L)+0.01]')
ax2.set_xticks(chtick)
ax2.set_xticklabels(chtickl)
ax2.set_yticks(chtick)
ax2.set_yticklabels(chtickl)
ax2.set_xlim(-2,1.8)
ax2.set_ylim(-2,1.8)

#fig.savefig('/home/eolson/pyCode/notebooks/figs/EvalDFOChl.png',dpi=200,transparent=True)

In [ ]: