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
%matplotlib inline
In [2]:
PATH= '/data/eolson/MEOPAR/SS36runs/CedarRuns/PAR_1_3/'
start_date = datetime.datetime(2017,1,1)
end_date = datetime.datetime(2017,12,31)
flen=10
namfmt='long'
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}
df1=et.loadPSF()
df1.head()
Out[2]:
| Lat | Lon | NO23 | PO4 | Si | Station | Z | dtUTC | Chl | Depth_m | Phaeo | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 49.608333 | -124.866667 | 18.65 | 1.459 | 68.43 | BS1 | 2.0 | 2015-02-18 23:07:00 | NaN | NaN | NaN |
| 1 | 49.608333 | -124.866667 | 26.62 | 2.175 | 55.73 | BS1 | 20.0 | 2015-02-18 23:07:00 | NaN | NaN | NaN |
| 2 | 49.608333 | -124.866667 | 16.25 | 1.809 | 39.14 | BS1 | 2.0 | 2015-03-21 22:55:00 | NaN | NaN | NaN |
| 3 | 49.608333 | -124.866667 | 16.83 | 2.069 | 39.11 | BS1 | 20.0 | 2015-03-21 22:55:00 | NaN | NaN | NaN |
| 4 | 49.608333 | -124.866667 | 24.27 | 0.551 | 49.90 | BS1 | 20.0 | 2015-04-05 23:10:00 | NaN | NaN | NaN |
In [3]:
data=et.matchData(df1,filemap, fdict, start_date, end_date, namfmt, PATH, flen)
In [4]:
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')
ax.plot(data.loc[(data.Lat>48.8)&(data.Lat<49.5)&(data.Lon>-124),['Lon']],data.loc[(data.Lat>48.8)&(data.Lat<49.5)&(data.Lon>-124),['Lat']],'s',color='c',label='Central SOG')
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);
Nitrate¶
In [6]:
data['N']=data['NO23']
In [7]:
N_s, modmean_s, obsmean_s, bias_s, RMSE_s, WSS_s = et.stats(data.loc[data.Z<15,['N']],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),['N']],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,['N']],data.loc[data.Z>=22,['mod_nitrate']])
N, modmean, obsmean, bias, RMSE, WSS = et.stats(data.loc[:,['N']],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: 573 bias: 1.9083500630181582 RMSE: 4.481007890431355 WSS: 0.8588145737803713 15 m<=z<22 m: N: 568 bias: 7.261366114383572 RMSE: 8.667287582341599 WSS: 0.5127090520837412 z>=22 m: N: 0 bias: nan RMSE: nan WSS: nan all: N: 1141 bias: 4.573129306817943 RMSE: 6.89057701174865 WSS: 0.7874124257784073
/home/eolson/anaconda3/envs/python36/lib/python3.6/site-packages/numpy/core/fromnumeric.py:2957: RuntimeWarning: Mean of empty slice. out=out, **kwargs) /home/eolson/anaconda3/envs/python36/lib/python3.6/site-packages/numpy/core/_methods.py:80: RuntimeWarning: invalid value encountered in double_scalars ret = ret.dtype.type(ret / rcount) /data/eolson/MEOPAR/tools/SalishSeaTools/salishsea_tools/evaltools.py:874: RuntimeWarning: invalid value encountered in double_scalars RMSE=np.sqrt(np.sum((mod-obs)**2)/N) /data/eolson/MEOPAR/tools/SalishSeaTools/salishsea_tools/evaltools.py:875: RuntimeWarning: invalid value encountered in double_scalars WSS=1.0-np.sum((mod-obs)**2)/np.sum((np.abs(mod-obsmean)+np.abs(obs-obsmean))**2)
In [8]:
fig, ax = plt.subplots(figsize = (8,8))
ps=et.varvarPlot(ax,data,'N','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)')
Out[8]:
<matplotlib.text.Text at 0x7f71f833c630>
In [9]:
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(2017,4,1)),:],'N','mod_nitrate',cols=('crimson','darkturquoise','navy'))
ax[0].set_title('Feb-Mar')
ii1=(data.Z < 15)&(data.dtUTC<=dt.datetime(2017,5,1))&(data.dtUTC>dt.datetime(2017,4,1))
ps=et.varvarPlot(ax[1],data.loc[ii1,:],'N','mod_nitrate',cols=('crimson','darkturquoise','navy'))
ax[1].set_title('April')
ii2=(data.Z < 15)&(data.dtUTC<=dt.datetime(2017,9,1))&(data.dtUTC>dt.datetime(2017,5,1))
ps=et.varvarPlot(ax[2],data.loc[ii2,:],'N','mod_nitrate',cols=('crimson','darkturquoise','navy'))
ax[2].set_title('May-Jun')
ii3=(data.Z < 15)&(data.dtUTC<=dt.datetime(2017,12,1))&(data.dtUTC>dt.datetime(2017,9,1))
ps=et.varvarPlot(ax[3],data.loc[ii3,:],'N','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(2017,4,1)),:],'N','mod_nitrate')
print('April:')
et.printstats(data.loc[ii1,:],'N','mod_nitrate')
print('May-Jun:')
et.printstats(data.loc[ii2,:],'N','mod_nitrate')
print('Sep-Oct:')
et.printstats(data.loc[ii3,:],'N','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: 77 bias: 3.0442629656605913 RMSE: 6.481879270684836 WSS: 0.6237513274616784 April: N: 70 bias: 1.0544263214562628 RMSE: 3.9748492087017597 WSS: 0.8273724322099776 May-Jun: N: 360 bias: 1.7237840170481902 RMSE: 3.944310717571498 WSS: 0.6849958265274843 Sep-Oct: N: 66 bias: 2.495518925973524 RMSE: 4.867465847073925 WSS: 0.8324978928074178
Out[9]:
[<matplotlib.lines.Line2D at 0x7f71ed1fe198>]
In [10]:
fig, ax = plt.subplots(figsize = (8,8))
viz_tools.set_aspect(ax, coords = 'map')
ax.plot(data['Lon'], data['Lat'], 'ro',label='data')
dJDF=data.loc[(data.Lon<-123.6)&(data.Lat<48.6)]
ax.plot(dJDF['Lon'],dJDF['Lat'],'b.',label='JDF')
dSJGI=data.loc[(data.Lon>=-123.6)&(data.Lat<48.9)]
ax.plot(dSJGI['Lon'],dSJGI['Lat'],'c.',label='SJGI')
dSOG=data.loc[(data.Lat>=48.9)&(data.Lon>-124.0)]
ax.plot(dSOG['Lon'],dSOG['Lat'],'y.',label='SOG')
dNSOG=data.loc[(data.Lat>=48.9)&(data.Lon<=-124.0)]
ax.plot(dNSOG['Lon'],dNSOG['Lat'],'m.',label='NSOG')
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);
In [11]:
fig, ax = plt.subplots(figsize = (8,8))
ps1=et.varvarPlot(ax,dJDF,'N','mod_nitrate',cols=('b','darkturquoise','navy'),lname='SJDF')
ps2=et.varvarPlot(ax,dSJGI,'N','mod_nitrate',cols=('c','darkturquoise','navy'),lname='SJGI')
ps3=et.varvarPlot(ax,dSOG,'N','mod_nitrate',cols=('y','darkturquoise','navy'),lname='SOG')
ps4=et.varvarPlot(ax,dNSOG,'N','mod_nitrate',cols=('m','darkturquoise','navy'),lname='NSOGF')
ax.legend(handles=[ps1[0][0],ps2[0][0],ps3[0][0],ps4[0][0]])
ax.set_xlabel('Obs')
ax.set_ylabel('Model')
ax.set_title('NO$_3$ ($\mu$M)')
ax.set_xlim(0,35)
ax.set_ylim(0,35)
Out[11]:
(0, 35)
In [12]:
fig, ax = plt.subplots(1,2,figsize = (17,8))
cols=('crimson','red','orangered','darkorange','gold','chartreuse','green','lightseagreen','cyan','darkturquoise','royalblue',
'lightskyblue','blue','darkblue','mediumslateblue','blueviolet','darkmagenta','fuchsia','deeppink','pink')
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,['N']].values, data.loc[iii,['Z']].values,
'.', color = cols[ii],label=str(ii))
ax[1].plot(data.loc[iii,['mod_silicon']].values-data.loc[iii,['Si']].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[0].set_xlim(-20,20)
ax[1].set_xlabel('model - obs Si')
ax[1].set_xlim(-40,20)
Out[12]:
(-40, 20)
Silicate¶
In [13]:
print('Nitrate')
print('z<15 m:')
et.printstats(data.loc[data.Z<15,:],'Si','mod_silicon')
print('15 m<=z<22 m:')
et.printstats(data.loc[(data.Z>=15)&(data.Z<22),:],'Si','mod_silicon')
print('z>=22 m:')
et.printstats(data.loc[data.Z>=22,:],'Si','mod_silicon')
print('all:')
et.printstats(data,'Si','mod_silicon')
print('obs Si < 50:')
et.printstats(data.loc[data.Si<50,:],'Si','mod_silicon')
Nitrate z<15 m: N: 563 bias: 11.123896681617211 RMSE: 17.57222464059505 WSS: 0.6577549687363178 15 m<=z<22 m: N: 558 bias: 0.8036770602303989 RMSE: 13.206223458331541 WSS: 0.5684751489937936 z>=22 m: N: 0 bias: nan RMSE: nan WSS: nan all: N: 1121 bias: 5.986802525744018 RMSE: 15.552922342633765 WSS: 0.6697189610756396 obs Si < 50: N: 892 bias: 10.330617699847483 RMSE: 15.482974743832532 WSS: 0.6368434932345437
/home/eolson/anaconda3/envs/python36/lib/python3.6/site-packages/numpy/core/fromnumeric.py:2957: RuntimeWarning: Mean of empty slice. out=out, **kwargs) /home/eolson/anaconda3/envs/python36/lib/python3.6/site-packages/numpy/core/_methods.py:80: RuntimeWarning: invalid value encountered in double_scalars ret = ret.dtype.type(ret / rcount) /data/eolson/MEOPAR/tools/SalishSeaTools/salishsea_tools/evaltools.py:874: RuntimeWarning: invalid value encountered in double_scalars RMSE=np.sqrt(np.sum((mod-obs)**2)/N) /data/eolson/MEOPAR/tools/SalishSeaTools/salishsea_tools/evaltools.py:875: RuntimeWarning: invalid value encountered in double_scalars WSS=1.0-np.sum((mod-obs)**2)/np.sum((np.abs(mod-obsmean)+np.abs(obs-obsmean))**2)
In [14]:
fig, ax = plt.subplots(figsize = (8,8))
ps=et.varvarPlot(ax,data,'Si','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)')
ax.set_xlim(0,80)
ax.set_ylim(0,80)
Out[14]:
(0, 80)
In [15]:
obsvar='Si'; 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(2017,4,1)),:],obsvar,modvar,cols=('crimson','darkturquoise','navy'))
ax[0].set_title('Feb-Mar')
ii1=(data.Z < 15)&(data.dtUTC<=dt.datetime(2017,5,1))&(data.dtUTC>dt.datetime(2017,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(2017,9,1))&(data.dtUTC>dt.datetime(2017,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(2017,12,1))&(data.dtUTC>dt.datetime(2017,9,1))
ps=et.varvarPlot(ax[3],data.loc[ii3,:],obsvar,modvar,cols=('crimson','darkturquoise','navy'))
ax[3].set_title('Sep-Oct')
ii4=(data.Z < 15)&(data.dtUTC<=dt.datetime(2017,4,1))&(data.dtUTC>dt.datetime(2017,2,1))
ps=et.varvarPlot(ax[0],data.loc[ii4,:],obsvar,modvar,cols=('darkturquoise','navy'))
ii5=(data.Z < 15)&(data.dtUTC<=dt.datetime(2017,5,1))&(data.dtUTC>dt.datetime(2017,4,1))
ps=et.varvarPlot(ax[1],data.loc[ii5,:],obsvar,modvar,cols=('darkturquoise','navy'))
print('Silicate, z<15')
print('Feb-Mar:')
et.printstats(data.loc[(data.Z<15)&(data.dtUTC<=dt.datetime(2017,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: 76 bias: 2.996606863258762 RMSE: 13.798612585063788 WSS: 0.4484065945199933 April: N: 67 bias: 6.053976564417596 RMSE: 16.0502446772289 WSS: 0.6397657842460402 May-Jun: N: 340 bias: 15.509598913147865 RMSE: 19.806464716401642 WSS: 0.5976945137009186 Sep-Oct: N: 80 bias: 4.451645623207099 RMSE: 10.448561257564423 WSS: 0.5127759241128553
Out[15]:
[<matplotlib.lines.Line2D at 0x7f71efb27b38>]
In [16]:
fig, ax = plt.subplots(figsize = (8,8))
ps1=et.varvarPlot(ax,dJDF,obsvar,modvar,cols=('b','darkturquoise','navy'),lname='SJDF')
ps2=et.varvarPlot(ax,dSJGI,obsvar,modvar,cols=('c','darkturquoise','navy'),lname='SJGI')
ps3=et.varvarPlot(ax,dSOG,obsvar,modvar,cols=('y','darkturquoise','navy'),lname='SOG')
ps4=et.varvarPlot(ax,dNSOG,obsvar,modvar,cols=('m','darkturquoise','navy'),lname='NSOGF')
ax.legend(handles=[ps1[0][0],ps2[0][0],ps3[0][0],ps4[0][0]])
ax.set_xlabel('Obs')
ax.set_ylabel('Model')
ax.set_title('Si ($\mu$M)')
ax.set_xlim(0,80)
ax.set_ylim(0,80)
Out[16]:
(0, 80)
Ratios¶
In [17]:
fig,ax=plt.subplots(1,2,figsize=(16,7))
p1=ax[0].plot(dJDF['N'],dJDF['Si'],'b.',label='SJDF')
p2=ax[0].plot(dSJGI['N'],dSJGI['Si'],'c.',label='SJGI')
p3=ax[0].plot(dSOG['N'],dSOG['Si'],'y.',label='SOG')
p4=ax[0].plot(dNSOG['N'],dNSOG['Si'],'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()
p5=ax[1].plot(dJDF['mod_nitrate'],dJDF['mod_silicon'],'b.',label='SJDF')
p6=ax[1].plot(dSJGI['mod_nitrate'],dSJGI['mod_silicon'],'c.',label='SJGI')
p7=ax[1].plot(dSOG['mod_nitrate'],dSOG['mod_silicon'],'y.',label='SOG')
p8=ax[1].plot(dNSOG['mod_nitrate'],dNSOG['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-')
Out[17]:
[<matplotlib.lines.Line2D at 0x7f71ed395400>]
In [18]:
data.loc[data.Si>65,['Month','Lat','Lon','Z','Si']]
/home/eolson/anaconda3/envs/python36/lib/python3.6/site-packages/pandas/core/indexing.py:1472: FutureWarning: Passing list-likes to .loc or [] with any missing label will raise KeyError in the future, you can use .reindex() as an alternative. See the documentation here: https://pandas.pydata.org/pandas-docs/stable/indexing.html#deprecate-loc-reindex-listlike return self._getitem_tuple(key)
Out[18]:
| Month | Lat | Lon | Z | Si | |
|---|---|---|---|---|---|
| 0 | NaN | 49.575000 | -124.155000 | 0.51 | 69.424140 |
| 66 | NaN | 49.206667 | -123.390000 | 0.51 | 71.259510 |
| 78 | NaN | 49.933333 | -124.983333 | 20.00 | 96.245450 |
| 107 | NaN | 49.485000 | -124.605000 | 20.00 | 101.304410 |
| 198 | NaN | 48.755000 | -123.591667 | 0.51 | 73.740000 |
| 256 | NaN | 48.855000 | -123.363333 | 20.00 | 161.472100 |
| 399 | NaN | 49.485000 | -124.605000 | 20.00 | 88.139875 |
| 401 | NaN | 49.483333 | -124.766667 | 20.00 | 75.064360 |
| 540 | NaN | 49.266667 | -123.293333 | 0.51 | 82.525050 |
| 544 | NaN | 49.206667 | -123.390000 | 0.51 | 70.429150 |
| 546 | NaN | 48.945000 | -123.320000 | 0.51 | 81.065200 |
| 554 | NaN | 48.856667 | -123.180000 | 0.51 | 85.444750 |
| 559 | NaN | 49.326667 | -123.305000 | 20.00 | 69.594950 |
| 674 | NaN | 49.266667 | -123.293333 | 0.51 | 75.017250 |
| 855 | NaN | 49.608333 | -124.866667 | 20.00 | 80.400000 |
| 980 | NaN | 49.986667 | -124.821667 | 20.00 | 69.800000 |
| 1050 | NaN | 49.903333 | -124.988333 | 20.00 | 65.500000 |
| 1054 | NaN | 49.608333 | -124.866667 | 20.00 | 77.130000 |
| 1123 | NaN | 49.608333 | -124.866667 | 20.00 | 79.100000 |
| 1169 | NaN | 49.783333 | -124.491667 | 20.00 | 68.200000 |
| 1203 | NaN | 49.691667 | -124.666667 | 20.00 | 73.500000 |
| 1207 | NaN | 49.608333 | -124.866667 | 20.00 | 74.400000 |
| 1209 | NaN | 49.483333 | -124.766667 | 20.00 | 65.300000 |
| 1253 | NaN | 50.020000 | -124.881667 | 20.00 | 67.650000 |
| 1257 | NaN | 49.921667 | -124.920000 | 20.00 | 69.700000 |
| 1269 | NaN | 49.608333 | -124.866667 | 20.00 | 69.300000 |
| 1273 | NaN | 49.783333 | -124.491667 | 20.00 | 70.100000 |
| 1321 | NaN | 48.746667 | -123.595000 | 20.00 | 66.900000 |
| 1340 | NaN | 49.608333 | -124.866667 | 20.00 | 69.900000 |
In [ ]:
Chlorophyll¶
In [21]:
data['l10_obsChl']=np.log10(data['Chl']+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 [28]:
print('log10[Chl+0.01]')
print('z<15 m:')
et.printstats(data.loc[data.Z<15,:],'l10_obsChl','l10_modChl')
print('z>=15 m:')
et.printstats(data.loc[data.Z>=15,:],'l10_obsChl','l10_modChl')
print('all:')
et.printstats(data,'l10_obsChl','l10_modChl')
print('\n')
print('Chl')
print('z<15 m:')
et.printstats(data.loc[data.Z<15,:],'Chl','mod_Chl')
print('z>=15 m:')
et.printstats(data.loc[data.Z>=15,:],'Chl','mod_Chl')
print('all:')
et.printstats(data,'Chl','mod_Chl')
log10[Chl+0.01] z<15 m: N: 0 bias: nan RMSE: nan WSS: nan z>=15 m: N: 122 bias: 0.2362057411871476 RMSE: 0.4332379692439433 WSS: 0.5910697719521386 all: N: 122 bias: 0.2362057411871476 RMSE: 0.4332379692439433 WSS: 0.5910697719521386 Chl z<15 m: N: 0 bias: nan RMSE: nan WSS: nan z>=15 m: N: 122 bias: 1.0309837117753284 RMSE: 3.9271917451621015 WSS: 0.5863013336453302 all: N: 122 bias: 1.0309837117753284 RMSE: 3.9271917451621015 WSS: 0.5863013336453302
/home/eolson/anaconda3/envs/python36/lib/python3.6/site-packages/numpy/core/fromnumeric.py:2957: RuntimeWarning: Mean of empty slice. out=out, **kwargs) /home/eolson/anaconda3/envs/python36/lib/python3.6/site-packages/numpy/core/_methods.py:80: RuntimeWarning: invalid value encountered in double_scalars ret = ret.dtype.type(ret / rcount) /data/eolson/MEOPAR/tools/SalishSeaTools/salishsea_tools/evaltools.py:874: RuntimeWarning: invalid value encountered in double_scalars RMSE=np.sqrt(np.sum((mod-obs)**2)/N) /data/eolson/MEOPAR/tools/SalishSeaTools/salishsea_tools/evaltools.py:875: RuntimeWarning: invalid value encountered in double_scalars WSS=1.0-np.sum((mod-obs)**2)/np.sum((np.abs(mod-obsmean)+np.abs(obs-obsmean))**2)
In [26]:
fig, ax = plt.subplots(1,2,figsize = (14,6))
ax[0].plot(np.arange(-.6,1.6,.1),np.arange(-.6,1.6,.1),'k-')
ps=et.varvarPlot(ax[0],data,'l10_obsChl','l10_modChl','Z',(5,10,15,20,25),'z','m',('crimson','darkorange','lime','mediumseagreen','darkturquoise','navy'))
ax[0].legend(handles=ps)
ax[0].set_xlabel('Obs')
ax[0].set_ylabel('Model')
ax[0].set_title('log10[Chl ($\mu$g/L)+0.01]')
ax[1].plot(np.arange(0,35),np.arange(0,35),'k-')
ps=et.varvarPlot(ax[1],data,'Chl','mod_Chl','Z',(5,10,15,20,25),'z','m',('crimson','darkorange','lime','mediumseagreen','darkturquoise','navy'))
ax[1].legend(handles=ps)
ax[1].set_xlabel('Obs')
ax[1].set_ylabel('Model')
ax[1].set_title('Chl ($\mu$g/L)')
Out[26]:
<matplotlib.text.Text at 0x7f71ed5b1ac8>
In [ ]:
In [ ]:
In [ ]: