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
PATH= '/data/eolson/MEOPAR/SS36runs/CedarRuns/spring2015_T3/'
start_date = datetime.datetime(2015,2,6)
#end_date = datetime.datetime(2015,10,4)
end_date = datetime.datetime(2015,11,12)
flen=10
varmap={'N':'nitrate','Si':'silicon','Ammonium':'ammonium'}
filemap={'nitrate':'ptrc_T','silicon':'ptrc_T','ammonium':'ptrc_T'}
gridmap={'nitrate':'tmask','silicon':'tmask','ammonium':'tmask'}
fdict={'ptrc_T':1,'grid_T':1}
df1=et.loadDFO()
df1.head()
| 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 | None | None | None | NaN | mg/m^3 | 15.31 | 32.14 | umol/L | 29.227507 | 9.7647 | 'deg_C_(ITS90)' | 1.884 | 2015-02-11 11:04:07 |
| 1 | 2015.0 | 2.0 | 11.0 | 11.068611 | 48.300833 | -124.000333 | 6.6 | None | None | None | 2.57 | mg/m^3 | 17.13 | 33.90 | umol/L | 29.484341 | 9.6880 | 'deg_C_(ITS90)' | 6.54434 | 2015-02-11 11:04:07 |
| 2 | 2015.0 | 2.0 | 11.0 | 11.068611 | 48.300833 | -124.000333 | 6.7 | None | None | None | NaN | mg/m^3 | NaN | NaN | umol/L | 29.484839 | 9.6828 | 'deg_C_(ITS90)' | 6.6435 | 2015-02-11 11:04:07 |
| 3 | 2015.0 | 2.0 | 11.0 | 11.068611 | 48.300833 | -124.000333 | 11.0 | None | None | None | NaN | mg/m^3 | NaN | NaN | umol/L | 30.144549 | 9.3646 | 'deg_C_(ITS90)' | 10.9071 | 2015-02-11 11:04:07 |
| 4 | 2015.0 | 2.0 | 11.0 | 11.068611 | 48.300833 | -124.000333 | 11.0 | None | None | None | NaN | mg/m^3 | 20.62 | 37.65 | umol/L | 30.157913 | 9.3586 | 'deg_C_(ITS90)' | 10.9071 | 2015-02-11 11:04:07 |
data,varmap=et.matchData(df1, varmap, filemap, fdict, start_date, end_date, 'long', PATH, flen)
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-7-3d9c623e9b0c> in <module>() ----> 1 data,varmap=et.matchData(df1, varmap, filemap, fdict, start_date, end_date, 'long', PATH, flen) /data/eolson/MEOPAR/tools/SalishSeaTools/salishsea_tools/evaltools.py in matchData(data, filemap, fdict, mod_start, mod_end, mod_nam_fmt, mod_basedir, mod_flen, method, deltat, deltad, meshPath, maskName) 117 118 # adjustments to data dataframe --> 119 data=data.loc[(data.dtUTC>=mod_start)&(data.dtUTC<mod_end)].copy(deep=True) 120 data=data.dropna(how='any',subset=['dtUTC','Lat','Lon','Z']) #.dropna(how='all',subset=[*varmap.keys()]) 121 data['j']=np.zeros((len(data))).astype(int) ~/anaconda3/envs/python36/lib/python3.6/site-packages/pandas/core/ops.py in wrapper(self, other, axis) 859 860 with np.errstate(all='ignore'): --> 861 res = na_op(values, other) 862 if is_scalar(res): 863 raise TypeError('Could not compare %s type with Series' % ~/anaconda3/envs/python36/lib/python3.6/site-packages/pandas/core/ops.py in na_op(x, y) 791 else: 792 mask = isnull(x) | isnull(y) --> 793 y = y.view('i8') 794 x = x.view('i8') 795 AttributeError: 'dict' object has no attribute 'view'
fig, ax = plt.subplots(figsize = (10,10))
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')
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);
list_of_cs_ni=data['N']
list_of_model_ni=data['mod_nitrate']
list_of_depths=data['Z']
list_of_lons=data['Lon']
list_of_lats=data['Lat']
list_of_datetimes=data['dtUTC']
list_of_cs_si=data['Si']
list_of_model_si=data['mod_silicon']
iii=np.logical_and(~np.isnan(list_of_cs_ni),~np.isnan(list_of_model_ni))
test=list_of_cs_ni[iii]
test2=list_of_model_ni[iii]
fig, ax = plt.subplots(figsize = (8,8))
ax.plot(list_of_cs_ni[list_of_depths < 15], list_of_model_ni[list_of_depths < 15],
'.', color = 'Crimson', alpha = 0.5, label = 'surface')
ax.plot(list_of_cs_ni[(22 > list_of_depths) & (list_of_depths >= 15)],
list_of_model_ni[(22 > list_of_depths) & (list_of_depths >= 15)],
'.', color = 'DarkOrange', alpha = 0.5, label = 'intermediate')
ax.plot(list_of_cs_ni[list_of_depths >= 22], list_of_model_ni[list_of_depths >= 22],
'.', color = 'darkblue', alpha = 0.5, label = 'deep')
ax.plot(list_of_cs_ni[(list_of_depths < 15) & (list_of_lons < -123.5) & (list_of_lats < 48.6)],
list_of_model_ni[(list_of_depths < 15) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
'.', color = 'ForestGreen', alpha = 0.5, label = ' Juan de Fuca surface')
ax.plot(list_of_cs_ni[(22 > list_of_depths) & (list_of_depths >= 15)
& (list_of_lons < -123.5) & (list_of_lats < 48.6)],
list_of_model_ni[(22 > list_of_depths) & (list_of_depths >= 15)
& (list_of_lons < -123.5) & (list_of_lats < 48.6)],
'.', color = 'CornflowerBlue', alpha = 0.5, label = ' Juan de Fuca intermediate')
ax.plot(list_of_cs_ni[(list_of_depths >= 22) & (list_of_lons < -123.5) & (list_of_lats < 48.6)],
list_of_model_ni[(list_of_depths >= 22) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
'.', color = 'DarkOrchid', alpha = 0.5, label = ' Juan de Fuca deep')
ax.plot(np.arange(0,35), color = 'grey')
ax.grid('on')
ax.set_title('DFO Nitrate')
ax.set_xlabel('DFO')
ax.set_ylabel('Nowcast-green');
ax.legend(loc=4)
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']])
print('N_s:', N_s,'modmean_s:',modmean_s, 'obsmean_s:',obsmean_s, 'bias_s"',bias_s, 'RMSE_s:',RMSE_s, 'WSS_s:', WSS_s)
print('surface bias = ' + str(-np.mean(list_of_cs_ni[list_of_depths < 15])
+ np.mean(list_of_model_ni[list_of_depths < 15])))
print('surface RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni[list_of_depths < 15]
- list_of_cs_ni[list_of_depths < 15])**2) /
len(list_of_cs_ni[list_of_depths < 15]))))
xbar = np.mean(list_of_cs_ni[list_of_depths < 15])
print('surface Willmott = ' + str(1-(np.sum((list_of_model_ni[list_of_depths < 15]
- list_of_cs_ni[list_of_depths < 15])**2) /
np.sum((np.abs(list_of_model_ni[list_of_depths < 15] - xbar)
+ np.abs(list_of_cs_ni[list_of_depths < 15] - xbar))**2))))
print('intermediate bias = ' + str(-np.mean(list_of_cs_ni[(22 > list_of_depths)
& (list_of_depths >= 15)])
+ np.mean(list_of_model_ni[(22 > list_of_depths)
& (list_of_depths >= 15)])))
print('intermediate RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni[(22 > list_of_depths)
& (list_of_depths >= 15)]
- list_of_cs_ni[(22 > list_of_depths)
& (list_of_depths >= 15)])**2) /
len(list_of_cs_ni[(22 > list_of_depths) & (list_of_depths >= 15)]))))
xbar = np.mean(list_of_cs_ni[(22 > list_of_depths) & (list_of_depths >= 15)])
print('intermediate Willmott = ' + str(1-(np.sum((list_of_model_ni[(22 > list_of_depths)
& (list_of_depths >= 15)]
- list_of_cs_ni[(22 > list_of_depths)
& (list_of_depths >= 15)])**2) /
np.sum((np.abs(list_of_model_ni[(22 > list_of_depths)
& (list_of_depths >= 15)] - xbar)
+ np.abs(list_of_cs_ni[(22 > list_of_depths)
& (list_of_depths >= 15)] - xbar))**2))))
print('deep bias = ' + str(-np.mean(list_of_cs_ni[list_of_depths >= 22])
+ np.mean(list_of_model_ni[list_of_depths >= 22])))
print('deep RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni[list_of_depths >= 22]
- list_of_cs_ni[list_of_depths >= 22])**2) /
len(list_of_cs_ni[list_of_depths >= 22]))))
xbar = np.mean(list_of_cs_ni[list_of_depths >= 22])
print('deep Willmott = ' + str(1-(np.sum((list_of_model_ni[list_of_depths >= 22]
- list_of_cs_ni[list_of_depths >= 22])**2) /
np.sum((np.abs(list_of_model_ni[list_of_depths >= 22] - xbar)
+ np.abs(list_of_cs_ni[list_of_depths >= 22] - xbar))**2))))
print('bias = ' + str(-np.mean(list_of_cs_ni) + np.mean(list_of_model_ni)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni - list_of_cs_ni)**2) /
len(list_of_cs_ni))))
xbar = np.mean(list_of_cs_ni)
print('Willmott = ' + str(1-(np.sum((list_of_model_ni - list_of_cs_ni)**2) /
np.sum((np.abs(list_of_model_ni - xbar)
+ np.abs(list_of_cs_ni - xbar))**2))))
fig, ax = plt.subplots(figsize = (8,8))
ps,ls=et.varvarPlot(ax,data,'N','mod_nitrate')
fig, ax = plt.subplots(figsize = (8,8))
ps,ls=et2.varvarPlot(ax,data,'N','mod_nitrate','Z',(15,22),'z','m',('mediumseagreen','skyblue','darkslateblue'))
ls
ps
ii=23
print(ii,listtemp[ii*10:(ii+1)*10])
data.loc[data.N==23.84]
with nc.Dataset('/ocean/eolson/MEOPAR/NEMO-forcing/grid/mesh_mask201702_noLPE.nc') as f:
dep=np.copy(f.variables['gdept_1d'])
tmask=np.copy(f.variables['tmask'])
lon=np.copy(f.variables['nav_lon'])
lat=np.copy(f.variables['nav_lat'])
tmask[0,:,344,292]
for ii in range(0,40):
print(ii,dep[0,ii],tmask[0,ii,344,292])
fig,ax=plt.subplots(1,4,figsize=(24,6))
ii=(list_of_depths < 15)&(list_of_datetimes<=dt.datetime(2015,4,1))
ax[0].plot(list_of_cs_ni[ii], list_of_model_ni[ii],
'.', color = 'Crimson', alpha = 0.5, label = 'surface')
ax[0].set_title('Feb-Mar')
print('Feb-Mar')
print('bias = ' + str(-np.mean(list_of_cs_ni[ii]) + np.mean(list_of_model_ni[ii])))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni[ii] - list_of_cs_ni[ii])**2) /
len(list_of_cs_ni[ii]))))
xbar = np.mean(list_of_cs_ni[ii])
print('Willmott = ' + str(1-(np.sum((list_of_model_ni[ii] - list_of_cs_ni[ii])**2) /
np.sum((np.abs(list_of_model_ni[ii] - xbar)
+ np.abs(list_of_cs_ni[ii] - xbar))**2))))
ii=(list_of_depths < 15)&(list_of_datetimes<=dt.datetime(2015,5,1))&(list_of_datetimes>dt.datetime(2015,4,1))
ax[1].plot(list_of_cs_ni[ii], list_of_model_ni[ii],
'.', color = 'Crimson', alpha = 0.5, label = 'surface')
ax[1].set_title('April')
print('April')
print('bias = ' + str(-np.mean(list_of_cs_ni[ii]) + np.mean(list_of_model_ni[ii])))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni[ii] - list_of_cs_ni[ii])**2) /
len(list_of_cs_ni[ii]))))
xbar = np.mean(list_of_cs_ni[ii])
print('Willmott = ' + str(1-(np.sum((list_of_model_ni[ii] - list_of_cs_ni[ii])**2) /
np.sum((np.abs(list_of_model_ni[ii] - xbar)
+ np.abs(list_of_cs_ni[ii] - xbar))**2))))
ii=(list_of_depths < 15)&(list_of_datetimes<=dt.datetime(2015,9,1))&(list_of_datetimes>dt.datetime(2015,5,1))
cs=ax[2].scatter(list_of_cs_ni[ii], list_of_model_ni[ii], c=np.array([np.float((ii-dt.datetime(2015,3,1)).days) for ii in list_of_datetimes[ii]]), alpha = 0.5, label = 'surface')
ax[2].set_title('May-Jun-Jul')
plt.colorbar(cs)
print('May-Jun-Jul')
print('bias = ' + str(-np.mean(list_of_cs_ni[ii]) + np.mean(list_of_model_ni[ii])))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni[ii] - list_of_cs_ni[ii])**2) /
len(list_of_cs_ni[ii]))))
xbar = np.mean(list_of_cs_ni[ii])
print('Willmott = ' + str(1-(np.sum((list_of_model_ni[ii] - list_of_cs_ni[ii])**2) /
np.sum((np.abs(list_of_model_ni[ii] - xbar)
+ np.abs(list_of_cs_ni[ii] - xbar))**2))))
ii=(list_of_depths < 15)&(list_of_datetimes<=dt.datetime(2015,12,1))&(list_of_datetimes>dt.datetime(2015,9,1))
#ax[3].plot(list_of_cs_ni[ii], list_of_model_ni[ii],
# '.', color = 'Crimson', alpha = 0.5, label = 'surface')
#ax[3].set_title('Sep-Oct')
#print('Sep-Oct')
#print('bias = ' + str(-np.mean(list_of_cs_ni[ii]) + np.mean(list_of_model_ni[ii])))
#print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni[ii] - list_of_cs_ni[ii])**2) /
# len(list_of_cs_ni[ii]))))
#xbar = np.mean(list_of_cs_ni[ii])
#print('Willmott = ' + str(1-(np.sum((list_of_model_ni[ii] - list_of_cs_ni[ii])**2) /
# np.sum((np.abs(list_of_model_ni[ii] - xbar)
# + np.abs(list_of_cs_ni[ii] - xbar))**2))))
for ii in range(0,4):
ax[ii].plot(np.arange(0,35), color = 'grey')
fig,ax=plt.subplots(1,1,figsize=(24,1))
plt.plot(list_of_datetimes,np.ones(np.shape(list_of_datetimes)),'k.')
fig, ax = plt.subplots(figsize = (8,8))
ax.plot(list_of_cs_ni[list_of_depths < 15], list_of_model_ni[list_of_depths < 15],
'.', color = 'Crimson', alpha = 0.8, label = 'surface')
ax.plot(list_of_cs_ni[(22 > list_of_depths) & (list_of_depths >= 15)],
list_of_model_ni[(22 > list_of_depths) & (list_of_depths >= 15)],
'.', color = 'DarkOrange', alpha = 0.8, label = 'intermediate')
ax.plot(list_of_cs_ni[(list_of_depths < 15) & (list_of_lons < -123.5) & (list_of_lats < 48.6)],
list_of_model_ni[(list_of_depths < 15) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
'.', color = 'ForestGreen', alpha = 0.8, label = ' Juan de Fuca surface')
ax.plot(list_of_cs_ni[(22 > list_of_depths) & (list_of_depths >= 15)
& (list_of_lons < -123.5) & (list_of_lats < 48.6)],
list_of_model_ni[(22 > list_of_depths) & (list_of_depths >= 15)
& (list_of_lons < -123.5) & (list_of_lats < 48.6)],
'.', color = 'CornflowerBlue', alpha = 0.8, label = ' Juan de Fuca intermediate')
ax.plot(np.arange(0,35), color = 'grey')
ax.grid('on')
ax.set_title('DFO Nitrate')
ax.set_xlabel('DFO')
ax.set_ylabel('Nowcast-green');
ax.legend(loc=4)
np.shape(list_of_depths)
fig, ax = plt.subplots(figsize = (8,8))
ax.plot(list_of_cs_si[list_of_depths < 15], list_of_model_si[list_of_depths < 15],
'.', color = 'Crimson', alpha = 0.5, label = 'surface')
ax.plot(list_of_cs_si[(22 > list_of_depths) & (list_of_depths >= 15)],
list_of_model_si[(22 > list_of_depths) & (list_of_depths >= 15)],
'.', color = 'DarkOrange', alpha = 0.5, label = 'intermediate')
ax.plot(list_of_cs_si[list_of_depths >= 22], list_of_model_si[list_of_depths >= 22],
'.', color = 'darkblue', alpha = 0.5, label = 'deep')
ax.plot(list_of_cs_si[(list_of_depths < 15) & (list_of_lons < -123.5) & (list_of_lats < 48.6)],
list_of_model_si[(list_of_depths < 15) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
'.', color = 'ForestGreen', alpha = 0.5, label = ' Juan de Fuca surface')
ax.plot(list_of_cs_si[(22 > list_of_depths) & (list_of_depths >= 15)
& (list_of_lons < -123.5) & (list_of_lats < 48.6)],
list_of_model_si[(22 > list_of_depths) & (list_of_depths >= 15)
& (list_of_lons < -123.5) & (list_of_lats < 48.6)],
'.', color = 'CornflowerBlue', alpha = 0.5, label = ' Juan de Fuca intermediate')
ax.plot(list_of_cs_si[(list_of_depths >= 22) & (list_of_lons < -123.5) & (list_of_lats < 48.6)],
list_of_model_si[(list_of_depths >= 22) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
'.', color = 'DarkOrchid', alpha = 0.5, label = ' Juan de Fuca deep')
ax.plot(np.arange(0,56), color = 'grey')
ax.grid('on')
ax.set_title('DFO Silicon')
ax.set_xlabel('DFO')
ax.set_ylabel('Nowcast-green');
ax.legend(loc=4)
print('surface bias = ' + str(-np.mean(list_of_cs_si[list_of_depths < 15])
+ np.mean(list_of_model_si[list_of_depths < 15])))
print('surface RMSE = ' + str(np.sqrt(np.sum((list_of_model_si[list_of_depths < 15]
- list_of_cs_si[list_of_depths < 15])**2) /
len(list_of_cs_si[list_of_depths < 15]))))
xbar = np.mean(list_of_cs_si[list_of_depths < 15])
print('surface Willmott = ' + str(1-(np.sum((list_of_model_si[list_of_depths < 15]
- list_of_cs_si[list_of_depths < 15])**2) /
np.sum((np.abs(list_of_model_si[list_of_depths < 15] - xbar)
+ np.abs(list_of_cs_si[list_of_depths < 15] - xbar))**2))))
print('intermediate bias = ' + str(-np.mean(list_of_cs_si[(22 > list_of_depths)
& (list_of_depths >= 15)])
+ np.mean(list_of_model_si[(22 > list_of_depths)
& (list_of_depths >= 15)])))
print('intermediate RMSE = ' + str(np.sqrt(np.sum((list_of_model_si[(22 > list_of_depths)
& (list_of_depths >= 15)]
- list_of_cs_si[(22 > list_of_depths)
& (list_of_depths >= 15)])**2) /
len(list_of_cs_si[(22 > list_of_depths) & (list_of_depths >= 15)]))))
xbar = np.mean(list_of_cs_si[(22 > list_of_depths) & (list_of_depths >= 15)])
print('intermediate Willmott = ' + str(1-(np.sum((list_of_model_si[(22 > list_of_depths)
& (list_of_depths >= 15)]
- list_of_cs_si[(22 > list_of_depths)
& (list_of_depths >= 15)])**2) /
np.sum((np.abs(list_of_model_si[(22 > list_of_depths)
& (list_of_depths >= 15)] - xbar)
+ np.abs(list_of_cs_si[(22 > list_of_depths)
& (list_of_depths >= 15)] - xbar))**2))))
print('deep bias = ' + str(-np.mean(list_of_cs_si[list_of_depths >= 22])
+ np.mean(list_of_model_si[list_of_depths >= 22])))
print('deep RMSE = ' + str(np.sqrt(np.sum((list_of_model_si[list_of_depths >= 22]
- list_of_cs_si[list_of_depths >= 22])**2) /
len(list_of_cs_si[list_of_depths >= 22]))))
xbar = np.mean(list_of_cs_si[list_of_depths >= 22])
print('deep Willmott = ' + str(1-(np.sum((list_of_model_si[list_of_depths >= 22]
- list_of_cs_si[list_of_depths >= 22])**2) /
np.sum((np.abs(list_of_model_si[list_of_depths >= 22] - xbar)
+ np.abs(list_of_cs_si[list_of_depths >= 22] - xbar))**2))))
print('bias = ' + str(-np.mean(list_of_cs_si) + np.mean(list_of_model_si)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_si - list_of_cs_si)**2) /
len(list_of_cs_si))))
xbar = np.mean(list_of_cs_si)
print('Willmott = ' + str(1-(np.sum((list_of_model_si - list_of_cs_si)**2) /
np.sum((np.abs(list_of_model_si - xbar)
+ np.abs(list_of_cs_si - xbar))**2))))
fig, ax = plt.subplots(1,2, figsize = (16,8))
ax[0].plot(list_of_cs_ni[list_of_depths >= 22], list_of_cs_si[list_of_depths >=22],
'.', color = 'darkblue', alpha = 0.5)
ax[0].plot(list_of_cs_ni[list_of_depths < 15], list_of_cs_si[list_of_depths < 15],
'.', color = 'Crimson', alpha = 0.5)
ax[0].plot(list_of_cs_ni[(22 > list_of_depths) & (list_of_depths >= 15)],
list_of_cs_si[(22 > list_of_depths) & (list_of_depths >= 15)],
'.', color = 'DarkOrange', alpha = 0.5)
ax[0].plot(list_of_cs_ni[(list_of_depths >= 22) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
list_of_cs_si[(list_of_depths >= 22) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
'.', color = 'DarkOrchid', alpha = 0.5)
ax[0].plot(list_of_cs_ni[(list_of_depths < 15) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
list_of_cs_si[(list_of_depths < 15) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
'.', color = 'ForestGreen', alpha = 0.5)
ax[0].plot(list_of_cs_ni[(22 > list_of_depths) & (list_of_depths >= 15)
& (list_of_lons < -123.5) & (list_of_lats < 48.6)],
list_of_cs_si[(22 > list_of_depths) & (list_of_depths >= 15)
& (list_of_lons < -123.5) & (list_of_lats < 48.6)],
'.', color = 'CornflowerBlue', alpha = 0.5)
#ax[0].plot(np.unique(list_of_cs_ni),
# np.poly1d(np.polyfit(list_of_cs_ni, list_of_cs_si, 1))(np.unique(list_of_cs_ni)))
x = np.arange(0,40)
ax[0].plot(x, 2*x, 'g-', alpha = 0.3)
ax[0].plot(x, 2*x+30, 'y-', alpha = 0.3)
ax[1].plot(list_of_model_ni[list_of_depths < 15], list_of_model_si[list_of_depths < 15],
'.', color = 'Crimson', alpha = 0.5, label = 'surface')
ax[1].plot(list_of_model_ni[(22 > list_of_depths) & (list_of_depths >= 15)],
list_of_model_si[(22 > list_of_depths) & (list_of_depths >= 15)],
'.', color = 'DarkOrange', alpha = 0.5, label = 'intermediate')
ax[1].plot(list_of_model_ni[list_of_depths >= 22], list_of_model_si[list_of_depths >=22],
'.', color = 'darkblue', alpha = 0.5, label = 'deep')
ax[1].plot(list_of_model_ni[(list_of_depths < 15) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
list_of_model_si[(list_of_depths < 15) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
'.', color = 'ForestGreen', alpha = 0.5, label = 'Juan de Fuca surface')
ax[1].plot(list_of_model_ni[(22 > list_of_depths) & (list_of_depths >= 15)
& (list_of_lons < -123.5) & (list_of_lats < 48.6)],
list_of_model_si[(22 > list_of_depths) & (list_of_depths >= 15)
& (list_of_lons < -123.5) & (list_of_lats < 48.6)],
'.', color = 'CornflowerBlue', alpha = 0.5, label = 'Juan de Fuca intermediate')
ax[1].plot(list_of_model_ni[(list_of_depths >= 22) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
list_of_model_si[(list_of_depths >= 22) & (list_of_lons < -123.5)
& (list_of_lats < 48.6)],
'.', color = 'DarkOrchid', alpha = 0.5, label = 'Juan de Fuca deep')
ax[1].plot(np.unique(list_of_model_ni),
np.poly1d(np.polyfit(list_of_model_ni,
list_of_model_si, 1))(np.unique(list_of_model_ni)))
x = np.arange(0,40)
ax[1].plot(x, 2*x, 'g-', alpha = 0.3, label = 'slope = 2')
ax[1].plot(x, 2*x+30, 'y-', alpha = 0.3, label = '2*N + 30')
ax[0].set_title('DFO', fontsize = 16)
ax[1].set_title('Model', fontsize = 16)
for ax in ax:
ax.grid('on')
ax.set_ylabel('Si', fontsize = 14)
ax.set_xlabel('N', fontsize = 14)
ax.set_ylim(0,95)
ax.set_xlim(0,37)
plt.legend();
fig, ax = plt.subplots(1,2, figsize = (16,8))
ax[0].plot(list_of_cs_ni, list_of_cs_si,
'.', color = 'darkblue', alpha = 0.5)
x = np.arange(0,40)
ax[0].plot(x, 2*x, 'g-', alpha = 0.3)
ax[0].plot(x, 2*x+30, 'y-', alpha = 0.3)
ax[1].plot(list_of_model_ni, list_of_model_si,
'.', color = 'Crimson', alpha = 0.5)
x = np.arange(0,40)
ax[1].plot(x, 2*x, 'g-', alpha = 0.3, label = 'slope = 2')
ax[1].plot(x, 2*x+30, 'y-', alpha = 0.3, label = '2*N + 30')
ax[0].set_title('DFO', fontsize = 16)
ax[1].set_title('Model', fontsize = 16)
for ax in ax:
ax.grid('on')
ax.set_ylabel('Si', fontsize = 14)
ax.set_xlabel('N', fontsize = 14)
ax.set_ylim(0,95)
ax.set_xlim(0,37)
plt.legend();
fig,ax=plt.subplots(1,4,figsize=(24,6))
ii=(list_of_depths < 15)&(list_of_datetimes<=dt.datetime(2015,4,1))
ax[0].plot(list_of_cs_si[ii], list_of_model_si[ii],
'.', color = 'Crimson', alpha = 0.5, label = 'surface')
ax[0].set_title('Feb-Mar')
print('Feb-Mar')
print('bias = ' + str(-np.mean(list_of_cs_si[ii]) + np.mean(list_of_model_si[ii])))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_si[ii] - list_of_cs_si[ii])**2) /
len(list_of_cs_si[ii]))))
xbar = np.mean(list_of_cs_si[ii])
print('Willmott = ' + str(1-(np.sum((list_of_model_si[ii] - list_of_cs_si[ii])**2) /
np.sum((np.abs(list_of_model_si[ii] - xbar)
+ np.abs(list_of_cs_si[ii] - xbar))**2))))
ii=(list_of_depths < 15)&(list_of_datetimes<=dt.datetime(2015,5,1))&(list_of_datetimes>dt.datetime(2015,4,1))
ax[1].plot(list_of_cs_si[ii], list_of_model_si[ii],
'.', color = 'Crimson', alpha = 0.5, label = 'surface')
ax[1].set_title('April')
print('April')
print('bias = ' + str(-np.mean(list_of_cs_si[ii]) + np.mean(list_of_model_si[ii])))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_si[ii] - list_of_cs_si[ii])**2) /
len(list_of_cs_si[ii]))))
xbar = np.mean(list_of_cs_si[ii])
print('Willmott = ' + str(1-(np.sum((list_of_model_si[ii] - list_of_cs_si[ii])**2) /
np.sum((np.abs(list_of_model_si[ii] - xbar)
+ np.abs(list_of_cs_si[ii] - xbar))**2))))
ii=(list_of_depths < 15)&(list_of_datetimes<=dt.datetime(2015,9,1))&(list_of_datetimes>dt.datetime(2015,5,1))
ax[2].plot(list_of_cs_si[ii], list_of_model_si[ii],
'.', color = 'Crimson', alpha = 0.5, label = 'surface')
ax[2].set_title('May-Jun-Jul')
print('May-Jun-Jul')
print('bias = ' + str(-np.mean(list_of_cs_si[ii]) + np.mean(list_of_model_si[ii])))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_si[ii] - list_of_cs_si[ii])**2) /
len(list_of_cs_si[ii]))))
xbar = np.mean(list_of_cs_si[ii])
print('Willmott = ' + str(1-(np.sum((list_of_model_si[ii] - list_of_cs_si[ii])**2) /
np.sum((np.abs(list_of_model_si[ii] - xbar)
+ np.abs(list_of_cs_si[ii] - xbar))**2))))
ii=(list_of_depths < 15)&(list_of_datetimes<=dt.datetime(2015,12,1))&(list_of_datetimes>dt.datetime(2015,9,1))
#ax[3].plot(list_of_cs_si[ii], list_of_model_si[ii],
# '.', color = 'Crimson', alpha = 0.5, label = 'surface')
#ax[3].set_title('Sep-Oct')
#print('Sep-Oct')
#print('bias = ' + str(-np.mean(list_of_cs_si[ii]) + np.mean(list_of_model_si[ii])))
#print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_si[ii] - list_of_cs_si[ii])**2) /
# len(list_of_cs_si[ii]))))
#xbar = np.mean(list_of_cs_si[ii])
#print('Willmott = ' + str(1-(np.sum((list_of_model_si[ii] - list_of_cs_si[ii])**2) /
# np.sum((np.abs(list_of_model_si[ii] - xbar)
# + np.abs(list_of_cs_si[ii] - xbar))**2))))
for ii in range(0,4):
ax[ii].plot(np.arange(0,35), color = 'grey')
fig, ax = plt.subplots(1,2,figsize = (17,8))
iii=(list_of_datetimes>dt.datetime(2015,2,1))
cols=('crimson','orangered','darkorange','gold','chartreuse','green','lightseagreen','cyan',
'lightskyblue','blue','mediumslateblue','blueviolet','darkmagenta','fuchsia')
for ii in range(2,12):
iii=(list_of_datetimes>=dt.datetime(2015,ii,1))&(list_of_datetimes<dt.datetime(2015,ii+1,1))
ax[0].plot(list_of_model_ni[iii]-list_of_cs_ni[iii], list_of_depths[iii],
'.', color = cols[ii],label=str(ii))
ax[1].plot(list_of_model_si[iii]-list_of_cs_si[iii], list_of_depths[iii],
'.', 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[0].set_xlim(-20,20)
ax[1].set_xlim(-40,40)