import numpy as np
import matplotlib.pyplot as plt
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, geo_tools
import glob
import gsw
import pytz
%matplotlib inline
f0 = pd.read_excel('/ocean/eolson/MEOPAR/obs/Hakai/Dosser20180911/2018-09-11_144804_HakaiData_nutrients.xlsx',
sheetname = 'Hakai Data')
f0.drop(['ACTION','Lat', 'Long', 'Collection Method', 'Installed', 'Lab Technician', 'NH4+', 'NO2+NO3 (ug/L)',
'no2_no3_units', 'TP', 'TDP', 'TN', 'TDN', 'SRP', 'Project Specific ID', 'Hakai ID', 'Source',
'po4pfilt', 'no3nfilt', 'po4punfl', 'no3nunfl', 'nh4nunfl', 'NH4+ Flag',
'TP FLag', 'TDP FLag', 'TN Flag', 'TDN FLag','Volume (ml)',
'SRP Flag', 'PO4 Flag', 'po4pfilt_flag', 'no3nfilt_flag','Preserved', 'Analyzed',
'po4punfl_flag', 'no3nunfl_flag', 'nh4nunfl_flag', 'Analyzing Lab', 'Sample Status',
'Quality Level', 'Comments', 'Quality Log'], axis = 1, inplace = True)
dts0=[pytz.timezone('Canada/Pacific').localize(i).astimezone(pytz.utc).replace(tzinfo=None)
for i in f0['Collected']]
f0['dtUTC']=dts0
f0.keys()
/home/eolson/anaconda3/envs/python36/lib/python3.6/site-packages/pandas/io/excel.py:329: FutureWarning: The `sheetname` keyword is deprecated, use `sheet_name` instead **kwds)
Index(['no', 'event_pk', 'Replicate Number', 'Date', 'Work Area', 'Survey', 'Sampling Bout', 'Site ID', 'Gather Lat', 'Gather Long', 'Line Out Depth', 'Pressure Transducer Depth (m)', 'Filter Type', 'Collected', 'NO2+NO3 (uM)', 'PO4', 'SiO2', 'NO2+NO3 Flag', 'SiO2 Flag', 'dtUTC'], dtype='object')
f0.head()
no | event_pk | Replicate Number | Date | Work Area | Survey | Sampling Bout | Site ID | Gather Lat | Gather Long | Line Out Depth | Pressure Transducer Depth (m) | Filter Type | Collected | NO2+NO3 (uM) | PO4 | SiO2 | NO2+NO3 Flag | SiO2 Flag | dtUTC | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 3483 | 4455 | 1 | 2014-11-07 | QUADRA | QOMA | 1 | QU5 | NaN | NaN | 0.0 | NaN | 0.45nm | 2014-11-07 11:30:39 | 17.528459 | 1.480 | 37.597988 | NaN | NaN | 2014-11-07 19:30:39 |
1 | 3484 | 4467 | 1 | 2014-11-07 | QUADRA | QOMA | 1 | QU17 | NaN | NaN | 0.0 | NaN | 0.45nm | 2014-11-07 10:30:30 | 19.534208 | 1.637 | 39.536296 | NaN | NaN | 2014-11-07 18:30:30 |
2 | 3485 | 4460 | 1 | 2014-11-07 | QUADRA | QOMA | 1 | QU24 | NaN | NaN | 0.0 | NaN | 0.45nm | 2014-11-07 09:30:32 | 20.243301 | 1.671 | 38.212329 | NaN | NaN | 2014-11-07 17:30:32 |
3 | 3486 | 4455 | 1 | 2014-11-07 | QUADRA | QOMA | 1 | QU5 | NaN | NaN | 5.0 | NaN | 0.45nm | 2014-11-07 11:30:39 | 19.414687 | 1.621 | 38.879558 | NaN | NaN | 2014-11-07 19:30:39 |
4 | 3487 | 4467 | 1 | 2014-11-07 | QUADRA | QOMA | 1 | QU17 | NaN | NaN | 5.0 | NaN | 0.45nm | 2014-11-07 10:30:30 | 19.728054 | 1.646 | 39.173135 | NaN | NaN | 2014-11-07 18:30:30 |
fc = pd.read_csv('/ocean/eolson/MEOPAR/obs/Hakai/Dosser20180911/ctd-bulk-1536702711696.csv',
usecols=['Cast PK','Cruise','Station', 'Drop number','Start time', 'Bottom time',
'Latitude', 'Longitude', 'Depth (m)', 'Temperature (deg C)', 'Temperature flag', 'Pressure (dbar)',
'Pressure flag', 'PAR', 'PAR flag', 'Fluorometry Chlorophyll (ug/L)', 'Fluorometry Chlorophyll flag',
'Turbidity (FTU)', 'Turbidity flag',
'Salinity (PSU)', 'Salinity flag'],
dtype={'Drop number':np.float64,'PAR flag':str,'Fluorometry Chlorophyll flag':str},na_values='null')
## fix apparent typos:
# reversed lats and lons
iii=fc['Latitude']>90
lons=-1*fc.loc[iii,'Latitude'].values
lats=-1*fc.loc[iii,'Longitude'].values
fc.loc[iii,'Longitude']=lons
fc.loc[iii,'Latitude']=lats
# remove data with missing lats and lons
nans=fc.loc[(fc['Latitude'].isnull())|(fc['Longitude'].isnull())]
fc=fc.drop(nans.index)
# apparently bad lats/lons
QU16bad=fc.loc[(fc['Station']=='QU16')&(fc['Latitude']>50.3)]
fc=fc.drop(QU16bad.index)
QU36bad=fc.loc[(fc['Station']=='QU36')&(fc['Latitude']>50.2)]
fc=fc.drop(QU36bad.index)
QU37bad=fc.loc[(fc['Station']=='QU37')&(fc['Longitude']<-125.1)]
fc=fc.drop(QU37bad.index)
QU38bad=fc.loc[(fc['Station']=='QU38')&(fc['Longitude']>-125.2)]
fc=fc.drop(QU38bad.index)
QU5bad=fc.loc[(fc['Station']=='QU5')&(fc['Longitude']>-125.18)]
fc=fc.drop(QU5bad.index)
#np.unique(fc.loc[(fc['Station']==sta)&(fc['Longitude']>-125.18)]['Cruise'].values)
fc.head()
Cast PK | Cruise | Station | Drop number | Start time | Bottom time | Latitude | Longitude | Depth (m) | Temperature (deg C) | ... | Pressure (dbar) | Pressure flag | PAR | PAR flag | Fluorometry Chlorophyll (ug/L) | Fluorometry Chlorophyll flag | Turbidity (FTU) | Turbidity flag | Salinity (PSU) | Salinity flag | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | QOM A | QU40 | 1.0 | 2015-11-03 12:37:20.000 | 2015-11-03 12:41:03.000 | 50.092435 | -125.164387 | 0.991 | 10.6879 | ... | 1 | NaN | 232.550 | NaN | 0.28830 | NaN | 0.742 | NaN | 28.1286 | NaN |
1 | 1 | QOM A | QU40 | 1.0 | 2015-11-03 12:37:20.000 | 2015-11-03 12:41:03.000 | 50.092435 | -125.164387 | 1.983 | 10.8076 | ... | 2 | NaN | 146.550 | NaN | 0.36906 | NaN | 0.827 | NaN | 28.2369 | NaN |
2 | 1 | QOM A | QU40 | 1.0 | 2015-11-03 12:37:20.000 | 2015-11-03 12:41:03.000 | 50.092435 | -125.164387 | 2.974 | 10.8952 | ... | 3 | NaN | 108.020 | NaN | 0.45831 | NaN | 0.742 | NaN | 28.3184 | NaN |
3 | 1 | QOM A | QU40 | 1.0 | 2015-11-03 12:37:20.000 | 2015-11-03 12:41:03.000 | 50.092435 | -125.164387 | 3.966 | 10.9617 | ... | 4 | NaN | 92.601 | NaN | 0.49377 | NaN | 0.776 | NaN | 28.3926 | NaN |
4 | 1 | QOM A | QU40 | 1.0 | 2015-11-03 12:37:20.000 | 2015-11-03 12:41:03.000 | 50.092435 | -125.164387 | 4.957 | 11.0704 | ... | 5 | NaN | 81.387 | NaN | 0.57567 | NaN | 0.718 | NaN | 28.5540 | NaN |
5 rows × 21 columns
fc.keys()
Index(['Cast PK', 'Cruise', 'Station', 'Drop number', 'Start time', 'Bottom time', 'Latitude', 'Longitude', 'Depth (m)', 'Temperature (deg C)', 'Temperature flag', 'Pressure (dbar)', 'Pressure flag', 'PAR', 'PAR flag', 'Fluorometry Chlorophyll (ug/L)', 'Fluorometry Chlorophyll flag', 'Turbidity (FTU)', 'Turbidity flag', 'Salinity (PSU)', 'Salinity flag'], dtype='object')
fc['dt']=[dt.datetime.strptime(i.split('.')[0],'%Y-%m-%d %H:%M:%S') for i in fc['Start time']]
fc['dt'][0].hour
12
hh=[ii.hour for ii in fc['dt']]
plt.hist(hh)
# looks like times are local vancouver time since people usually sample during daylight....
(array([ 2769., 11727., 26700., 68449., 86168., 83460., 67215., 37648., 11945., 3597.]), array([ 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16.]), <a list of 10 Patch objects>)
fc['dt']=[dt.datetime.strptime(i.split('.')[0],'%Y-%m-%d %H:%M:%S') for i in fc['Start time']]
dts=[pytz.timezone('Canada/Pacific').localize(dt.datetime.strptime(i.split('.')[0],'%Y-%m-%d %H:%M:%S')).astimezone(pytz.utc).replace(tzinfo=None)
for i in fc['Start time']]
fc['dtUTC']=dts
dloc=[dt.datetime(i.year,i.month,i.day) for i in fc['dt']]
fc['dloc']=dloc
48/8
6.0
fig,ax=plt.subplots(4,2,figsize=(18,36))
axf=ax.flatten()
grid = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
for axi in axf:
viz_tools.set_aspect(axi, coords = 'map')
viz_tools.plot_coastline(axi, grid, coords = 'map')
axi.set_xlim(-125.5,-124.8)
axi.set_ylim(50,50.5)
cols=['red','darkorange','goldenrod','chartreuse','forestgreen','turquoise','darkviolet']
count=0
axi=0
namelist=list()
ps=list()
for name, group in fc.groupby('Station'):
meanlon=np.nanmean(group['Longitude'].values)
meanlat=np.nanmean(group['Latitude'].values)
if count==6:
axf[axi].legend(handles=ps)
axi=axi+1
count=0
namelist=list()
ps=list()
# print the name of the regiment
ll=np.unique([(i,j) for i,j in zip(group['Latitude'].values,group['Longitude'].values)],axis=0)
for el in ll:
p,=axf[axi].plot(el[1],el[0],'.',color=cols[count%6],label=name)
if not name in namelist:
namelist.append(name)
ps.append(p)
axf[axi].plot(meanlon,meanlat,'*',color=cols[count%6],alpha=.5,markersize=20)
count=count+1
axf[axi].legend(handles=ps)
<matplotlib.legend.Legend at 0x7f1d03e56940>
fcS=fc.loc[:,['Latitude','Longitude']].groupby([fc['Station'],fc['dloc']]).mean().reset_index()
fcS.head()
Station | dloc | Latitude | Longitude | |
---|---|---|---|---|
0 | QU 23 | 2014-08-28 | 50.030000 | -125.120000 |
1 | QU 23 | 2014-10-09 | 50.030000 | -125.120000 |
2 | QU 23 | 2014-10-29 | 50.030000 | -125.120000 |
3 | QU08 | 2016-12-01 | 50.146282 | -125.175323 |
4 | QU10 | 2014-08-15 | 50.232300 | -125.120700 |
f0['Station']=f0['Site ID']
#f0['dt']=[dt.datetime.strptime(i,'%Y-%m-%d %H:%M:%S') for i in f0['Collected']]
dloc0=[dt.datetime(i.year,i.month,i.day) for i in f0['Collected']]
f0['dloc']=dloc0
fdata=f0.merge(fcS,how='left')
len(fdata.loc[fdata['Latitude'].isnull()]),len(fdata.loc[fdata['Latitude']>=0])
(247, 3997)
fdata['Lat']=fdata['Latitude']
fdata['Lon']=fdata['Longitude']
fdata['Z']=fdata['Line Out Depth']
PATH= '/results/SalishSea/hindcast.201812/'
start_date = dt.datetime(2016,1,1)
end_date = dt.datetime(2016,9,20)
flen=1
namfmt='nowcast'
filemap={'nitrate':'ptrc_T','silicon':'ptrc_T','ammonium':'ptrc_T','diatoms':'ptrc_T','ciliates':'ptrc_T','flagellates':'ptrc_T'}
fdict={'ptrc_T':1,'grid_T':1}
fdata2=fdata.loc[(fdata['dtUTC']>start_date)&(fdata['dtUTC']<end_date)&(fdata['Z']>=0)&(fdata['Z']<440)&(fdata['Lon']<360)&(fdata['Lat']<=90)].copy(deep=True).reset_index()
data=et.matchData(fdata2,filemap, fdict, start_date, end_date, namfmt, PATH, flen)
data.keys()
Index(['index', 'no', 'event_pk', 'Replicate Number', 'Date', 'Work Area', 'Survey', 'Sampling Bout', 'Site ID', 'Gather Lat', 'Gather Long', 'Line Out Depth', 'Pressure Transducer Depth (m)', 'Filter Type', 'Collected', 'NO2+NO3 (uM)', 'PO4', 'SiO2', 'NO2+NO3 Flag', 'SiO2 Flag', 'dtUTC', 'Station', 'dloc', 'Latitude', 'Longitude', 'Lat', 'Lon', 'Z', 'j', 'i', 'mod_nitrate', 'mod_silicon', 'mod_ammonium', 'mod_diatoms', 'mod_ciliates', 'mod_flagellates'], dtype='object')
N_s, modmean_s, obsmean_s, bias_s, RMSE_s, WSS_s = et.stats(data.loc[data.Z<15,['NO2+NO3 (uM)']],data.loc[data.Z<15,['mod_nitrate']])
N_i, modmean_i, obsmean_i, bias_i, RMSE_i, WSS_i = et.stats(data.loc[(data.Z>=15)&(data.Z<22),['NO2+NO3 (uM)']],data.loc[(data.Z>=15)&(data.Z<22),['mod_nitrate']])
N_d, modmean_d, obsmean_d, bias_d, RMSE_d, WSS_d = et.stats(data.loc[data.Z>=22,['NO2+NO3 (uM)']],data.loc[data.Z>=22,['mod_nitrate']])
N, modmean, obsmean, bias, RMSE, WSS = et.stats(data.loc[:,['NO2+NO3 (uM)']],data.loc[:,['mod_nitrate']])
print('Nitrate')
print('z<15 m:')
print(' N: {}\n bias: {}\n RMSE: {}\n WSS: {}'.format(N_s,bias_s,RMSE_s,WSS_s))
print('15 m<=z<22 m:')
print(' N: {}\n bias: {}\n RMSE: {}\n WSS: {}'.format(N_i,bias_i,RMSE_i,WSS_i))
print('z>=22 m:')
print(' N: {}\n bias: {}\n RMSE: {}\n WSS: {}'.format(N_d,bias_d,RMSE_d,WSS_d))
print('all:')
print(' N: {}\n bias: {}\n RMSE: {}\n WSS: {}'.format(N,bias,RMSE,WSS))
Nitrate z<15 m: N: 333 bias: -3.952714515595355 RMSE: 6.70460106001521 WSS: 0.8784532665618973 15 m<=z<22 m: N: 36 bias: -0.6949500264641557 RMSE: 4.887436118313972 WSS: 0.7310261390693128 z>=22 m: N: 335 bias: -2.935951217408032 RMSE: 3.5906368098522243 WSS: 0.8155143932461633 all: N: 704 bias: -3.3022951597693933 RMSE: 5.349687543883768 WSS: 0.9043167952085497
fig, ax = plt.subplots(figsize = (8,8))
ps=et.varvarPlot(ax,data,'NO2+NO3 (uM)','mod_nitrate','Z',(15,22),'z','m',('mediumseagreen','darkturquoise','navy'))
ax.legend(handles=ps)
ax.set_xlabel('Obs')
ax.set_ylabel('Model')
ax.set_title('NO$_3$ ($\mu$M)')
<matplotlib.text.Text at 0x7f1cfe4bc160>
nvar='NO2+NO3 (uM)'
fig, ax = plt.subplots(1,4,figsize = (24,6))
for axi in ax:
axi.plot(np.arange(0,30),np.arange(0,30),'k-')
ps=et.varvarPlot(ax[0],data.loc[(data.Z<15)&(data.dtUTC<=dt.datetime(2016,4,1)),:],nvar,'mod_nitrate',cols=('crimson','darkturquoise','navy'))
ax[0].set_title('Feb-Mar')
ii1=(data.Z < 15)&(data.dtUTC<=dt.datetime(2016,5,1))&(data.dtUTC>dt.datetime(2016,4,1))
ps=et.varvarPlot(ax[1],data.loc[ii1,:],nvar,'mod_nitrate',cols=('crimson','darkturquoise','navy'))
ax[1].set_title('April')
ii2=(data.Z < 15)&(data.dtUTC<=dt.datetime(2016,9,1))&(data.dtUTC>dt.datetime(2016,5,1))
ps=et.varvarPlot(ax[2],data.loc[ii2,:],nvar,'mod_nitrate',cols=('crimson','darkturquoise','navy'))
ax[2].set_title('May-Jun')
ii3=(data.Z < 15)&(data.dtUTC<=dt.datetime(2016,12,1))&(data.dtUTC>dt.datetime(2016,9,1))
ps=et.varvarPlot(ax[3],data.loc[ii3,:],nvar,'mod_nitrate',cols=('crimson','darkturquoise','navy'))
ax[3].set_title('Sep-Oct')
print('Nitrate, z<15')
print('Feb-Mar:')
et.printstats(data.loc[(data.Z<15)&(data.dtUTC<=dt.datetime(2016,4,1)),:],nvar,'mod_nitrate')
print('April:')
et.printstats(data.loc[ii1,:],nvar,'mod_nitrate')
print('May-Jun:')
et.printstats(data.loc[ii2,:],nvar,'mod_nitrate')
print('Sep-Oct:')
et.printstats(data.loc[ii3,:],nvar,'mod_nitrate')
fig,ax=plt.subplots(1,1,figsize=(24,1))
plt.plot(data.dtUTC,np.ones(np.shape(data.dtUTC)),'k.')
Nitrate, z<15 Feb-Mar: N: 191 bias: -7.09526086726769 RMSE: 7.818988259088877 WSS: 0.3892853308052122 April: N: 35 bias: 0.2855394217410421 RMSE: 5.386738762951721 WSS: 0.7990566637058589 May-Jun: N: 91 bias: 0.06683711461279174 RMSE: 4.614664443438714 WSS: 0.901603908205006 Sep-Oct: N: 16 bias: 1.4290521727610122 RMSE: 4.598719714508531 WSS: 0.850352982559323
[<matplotlib.lines.Line2D at 0x7f1cfc339780>]
sivar='SiO2'
fig, ax = plt.subplots(1,2,figsize = (17,8))
cols=('crimson','orangered','darkorange','gold','chartreuse','green','lightseagreen','cyan',
'lightskyblue','blue','mediumslateblue','blueviolet','darkmagenta','fuchsia')
ii0=start_date
for ii in range(0,int((end_date-start_date).days/30)):
iii=(data.dtUTC>=(start_date+dt.timedelta(days=ii*30)))&(data.dtUTC<(start_date+dt.timedelta(days=(ii+1)*30)))
ax[0].plot(data.loc[iii,['mod_nitrate']].values-data.loc[iii,[nvar]].values, data.loc[iii,['Z']].values,
'.', color = cols[ii],label=str(ii))
ax[1].plot(data.loc[iii,['mod_silicon']].values-data.loc[iii,[sivar]].values, data.loc[iii,['Z']].values,
'.', color = cols[ii],label=str(ii))
for axi in (ax[0],ax[1]):
axi.legend(loc=4)
axi.set_ylim(400,0)
axi.set_ylabel('depth (m)')
ax[0].set_xlabel('model - obs N')
ax[1].set_xlabel('model - obs Si')
ax[1].set_xlim(-40,20)
(-40, 20)
print('Nitrate')
print('z<15 m:')
et.printstats(data.loc[data.Z<15,:],sivar,'mod_silicon')
print('15 m<=z<22 m:')
et.printstats(data.loc[(data.Z>=15)&(data.Z<22),:],sivar,'mod_silicon')
print('z>=22 m:')
et.printstats(data.loc[data.Z>=22,:],sivar,'mod_silicon')
print('all:')
et.printstats(data,sivar,'mod_silicon')
Nitrate z<15 m: N: 331 bias: 2.492966389072329 RMSE: 11.365341568768184 WSS: 0.8350741856727155 15 m<=z<22 m: N: 36 bias: 5.570504666604677 RMSE: 10.32896988683681 WSS: 0.6658440932374314 z>=22 m: N: 335 bias: -2.1107273596716 RMSE: 6.227630055184728 WSS: 0.7509855281035017 all: N: 702 bias: 0.4538694833201262 RMSE: 9.213267918398328 WSS: 0.8490746296048648
fig, ax = plt.subplots(figsize = (8,8))
ps=et.varvarPlot(ax,data,sivar,'mod_silicon','Z',(15,22),'z','m',('mediumseagreen','darkturquoise','navy'))
ax.legend(handles=ps)
ax.set_xlabel('Obs')
ax.set_ylabel('Model')
ax.set_title('dSi ($\mu$M)')
<matplotlib.text.Text at 0x7f1d0edbfc18>
obsvar=sivar; modvar='mod_silicon'
fig, ax = plt.subplots(1,4,figsize = (24,6))
for axi in ax:
axi.plot(np.arange(0,70),np.arange(0,70),'k-')
ps=et.varvarPlot(ax[0],data.loc[(data.Z<15)&(data.dtUTC<=dt.datetime(2016,4,1)),:],obsvar,modvar,cols=('crimson','darkturquoise','navy'))
ax[0].set_title('Feb-Mar')
ii1=(data.Z < 15)&(data.dtUTC<=dt.datetime(2016,5,1))&(data.dtUTC>dt.datetime(2016,4,1))
ps=et.varvarPlot(ax[1],data.loc[ii1,:],obsvar,modvar,cols=('crimson','darkturquoise','navy'))
ax[1].set_title('April')
ii2=(data.Z < 15)&(data.dtUTC<=dt.datetime(2016,9,1))&(data.dtUTC>dt.datetime(2016,5,1))
ps=et.varvarPlot(ax[2],data.loc[ii2,:],obsvar,modvar,cols=('crimson','darkturquoise','navy'))
ax[2].set_title('May-Jun')
ii3=(data.Z < 15)&(data.dtUTC<=dt.datetime(2016,12,1))&(data.dtUTC>dt.datetime(2016,9,1))
ps=et.varvarPlot(ax[3],data.loc[ii3,:],obsvar,modvar,cols=('crimson','darkturquoise','navy'))
ax[3].set_title('Sep-Oct')
print('Silicate, z<15')
print('Feb-Mar:')
et.printstats(data.loc[(data.Z<15)&(data.dtUTC<=dt.datetime(2016,4,1)),:],obsvar,modvar)
print('April:')
et.printstats(data.loc[ii1,:],obsvar,modvar)
print('May-Jun:')
et.printstats(data.loc[ii2,:],obsvar,modvar)
print('Sep-Oct:')
et.printstats(data.loc[ii3,:],obsvar,modvar)
fig,ax=plt.subplots(1,1,figsize=(24,1))
plt.plot(data.dtUTC,np.ones(np.shape(data.dtUTC)),'k.')
Silicate, z<15 Feb-Mar: N: 189 bias: -4.524216587050645 RMSE: 6.6633051687484635 WSS: 0.35892182139320417 April: N: 35 bias: 10.942657450570477 RMSE: 14.649967485631025 WSS: 0.6442833174958666 May-Jun: N: 91 bias: 11.066768090494993 RMSE: 14.99628617031127 WSS: 0.7174794378844713 Sep-Oct: N: 16 bias: 18.136243920656305 RMSE: 19.980380272194136 WSS: 0.4619698161539917
[<matplotlib.lines.Line2D at 0x7f1cfc128630>]
fig,ax=plt.subplots(1,2,figsize=(16,7))
p4=ax[0].plot(data[nvar],data[sivar],'m.',label='NSOG')
ax[0].set_title('Observed')
ax[0].set_xlabel('N')
ax[0].set_ylabel('Si')
ax[0].set_xlim(0,40)
ax[0].set_ylim(0,85)
ax[0].legend()
p8=ax[1].plot(data['mod_nitrate'],data['mod_silicon'],'m.',label='NSOG')
ax[1].set_title('Model')
ax[1].set_xlabel('N')
ax[1].set_ylabel('Si')
ax[1].set_xlim(0,40)
ax[1].set_ylim(0,85)
ax[1].legend()
ax[0].plot(np.arange(0,35),1.3*np.arange(0,35),'k-')
ax[1].plot(np.arange(0,35),1.3*np.arange(0,35),'k-')
[<matplotlib.lines.Line2D at 0x7f1d0e8dd780>]