import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from salishsea_tools import (tidetools, geo_tools, viz_tools, gsw_calls)
import numpy.ma as ma
import pandas as pd
import datetime
import pytz
import os
%matplotlib inline
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
grid = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
bathy, X, Y = tidetools.get_bathy_data(grid)
mesh = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/mesh_mask201702.nc')
HINDCAST_PATH= '/results/SalishSea/nowcast-green/'
f = pd.read_csv('/ocean/eolson/MEOPAR/obs/NANOOS_SJPEF/ticket2180-csv/Discrete.csv')
f.keys()
Index(['year', 'month', 'day', 'time', 'depth', 'station', 'chl', 'phaeo', 'o2', 'po4', 'sio4', 'no3', 'no2', 'nh4'], dtype='object')
f.shape
(1031, 14)
f.station.unique()
array(['I', 'O', 'N', 'S', 'T'], dtype=object)
f.depth.unique()
array([ 30., 20., 15., 10., 5., 0., 80., 50., 125., 90., 70., 95., 120., 105., 85., 110., 100., 60., 40., 29., 115., 96., 86., 88., 75.])
g = pd.read_excel('/ocean/eolson/MEOPAR/obs/NANOOS_SJPEF/SanJuansPEF_dataset/SanJuansPEF_dataset_Inventory.xls',
sheetname = 1)
stations = {}
for station in f.station.unique():
Yind, Xind = geo_tools.find_closest_model_point(g[g.code == station].lon.values[0],
g[g.code == station].lat.values[0],
X, Y, land_mask = bathy.mask)
stations[station] = (Yind, Xind)
f.dropna(subset=['year', 'month', 'day', 'depth', 'station', 'chl']).shape
(784, 14)
f.dropna(subset=['year', 'month', 'day', 'depth', 'station', 'sio4']).shape
(150, 14)
f.dropna(subset=['year', 'month', 'day', 'depth', 'station', 'no3']).shape
(150, 14)
import datetime
deptht = (nc.Dataset(
'/results/SalishSea/nowcast-green/01jan18/SalishSea_1d_20180101_20180101_dia2_T.nc')
.variables['deptht'][:])
f2 = f.dropna(subset=['month', 'day', 'depth', 'station', 'no3', 'sio4'])
f2.station.unique()
array(['I', 'O', 'S', 'N'], dtype=object)
f2.year.unique()
array([2008, 2006, 2009, 2005, 2007])
f2.shape
(150, 14)
datetimes = np.array([datetime.date(2017,f2.month[n], f2.day[n]) for n in f2.index])
datetimes.min()
datetime.date(2017, 9, 27)
datetimes.max()
datetime.date(2017, 10, 23)
list_of_Yinds = np.array([])
list_of_Xinds = np.array([])
list_of_datetimes = np.array([])
list_of_cs_ni = np.array([])
list_of_cs_si = np.array([])
list_of_model_ni = np.array([])
list_of_model_si = np.array([])
list_of_depths = np.array([])
for n in f2.index:
Yind, Xind = stations[f2.station[n]]
depth = np.argmin(np.abs(deptht - f2.depth[n]))
if mesh.variables['tmask'][0,depth,Yind, Xind] == 1:
for year in [2014,2015,2016,2017]:
date = datetime.date(year, f2.month[n], f2.day[n])
sub_dir = date.strftime('%d%b%y').lower()
datestr = date.strftime('%Y%m%d')
fname = 'SalishSea_1d_{}_{}_ptrc_T.nc'.format(datestr, datestr)
nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
ni_val = nuts.variables['nitrate'][0, depth, Yind, Xind]
si_val = nuts.variables['silicon'][0, depth, Yind, Xind]
list_of_Yinds = np.append(list_of_Yinds, Yind)
list_of_Xinds = np.append(list_of_Xinds, Xind)
list_of_datetimes = np.append(list_of_datetimes, date)
list_of_cs_ni = np.append(list_of_cs_ni, float(f2['no3'][n]))
list_of_cs_si = np.append(list_of_cs_si, float(f2['sio4'][n]))
list_of_model_ni = np.append(list_of_model_ni, ni_val)
list_of_model_si = np.append(list_of_model_si, si_val)
list_of_depths = np.append(list_of_depths, depth)
list_of_years = np.array([])
for n in f2.index:
Yind, Xind = stations[f2.station[n]]
depth = np.argmin(np.abs(deptht - f2.depth[n]))
if mesh.variables['tmask'][0,depth,Yind, Xind] == 1:
for year in [2014,2015,2016,2017]:
list_of_years = np.append(list_of_years, int(f2.year[n]))
list_of_cs_ni.shape
(564,)
list_of_model_ni.shape
(564,)
564/4
141.0
fig, ax = plt.subplots(figsize = ((10,8)))
viz_tools.plot_coastline(ax, grid)
viz_tools.set_aspect(ax)
ax.set_ylim(250,350)
ax.set_xlim(200,300)
ax.plot(list_of_Xinds, list_of_Yinds, 'ro');
np.unique(list_of_depths)
array([ 0., 4., 9., 18., 21., 23., 24., 25., 26., 27.])
np.unique(list_of_depths)
array([ 0., 4., 9., 18., 21., 23., 24., 25., 26., 27.])
np.unique(list_of_depths2)
array([ 0., 4., 9., 14., 18., 21., 23., 24., 25.])
import matplotlib.cm as cm
import cmocean
fig, ax = plt.subplots(figsize = (8,8))
colours = cmocean.cm.deep(np.linspace(0,1,11))
for depth, colour in zip(np.unique(list_of_depths), colours[1:]):
ax.plot(list_of_cs_si[list_of_depths == depth], list_of_model_si[list_of_depths==depth],
'o', color = colour, alpha = 0.5, label = 'depth = '+str(depth))
ax.plot(np.arange(40,65), np.arange(40,65), color = 'grey')
ax.set_title('Si')
ax.set_xlabel('Observed')
ax.set_ylabel('Modelled')
ax.legend();
np.unique(list_of_years)
array([2005., 2006., 2007., 2008., 2009.])
p.year.unique()
array([2014, 2015, 2017, 2016])
fig, ax = plt.subplots(figsize = (8,8))
colours = cm.tab10(np.linspace(0,1,10))
for year, colour in zip(np.unique(list_of_years), colours):
ax.plot(list_of_cs_si[list_of_years == year], list_of_model_si[list_of_years==year],
'o', color = colour, alpha = 0.5, label = 'year = '+str(year))
ax.plot(np.arange(40,65), np.arange(40,65), color = 'grey')
ax.set_title('Si')
ax.set_xlabel('Observed')
ax.set_ylabel('Modelled')
ax.legend();
fig, ax = plt.subplots(figsize = (8,8))
colours = cmocean.cm.deep(np.linspace(0,1,11))
for depth, colour in zip(np.unique(list_of_depths), colours[1:]):
ax.plot(list_of_cs_ni[list_of_depths == depth], list_of_model_ni[list_of_depths==depth],
'o', color = colour, alpha = 0.5, label = 'depth = '+str(depth))
ax.plot(np.arange(18,33), np.arange(18,33), color = 'grey')
ax.set_title('N')
ax.set_xlabel('Observed')
ax.set_ylabel('Modelled')
ax.legend();
fig, ax = plt.subplots(figsize = (8,8))
colours = cm.tab10(np.linspace(0,1,10))
for year, colour in zip(np.unique(list_of_years), colours):
ax.plot(list_of_cs_ni[list_of_years == year], list_of_model_ni[list_of_years==year],
'o', color = colour, alpha = 0.5, label = 'year = '+str(year))
ax.plot(np.arange(18,33), np.arange(18,33), color = 'grey')
ax.set_title('N')
ax.set_xlabel('Observed')
ax.set_ylabel('Modelled')
ax.legend();
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(list_of_depths, list_of_cs_ni - list_of_model_ni,'.', alpha = 0.5, label = 'N')
ax.plot(list_of_depths, list_of_cs_si - list_of_model_si,'.', alpha = 0.5, label = 'Si')
ax.grid('on')
ax.set_xlabel('Model depth level')
ax.set_ylabel('Observed - Modelled')
ax.legend();
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(list_of_cs_ni - list_of_model_ni,list_of_Yinds,'.', alpha = 0.5, label = 'N')
ax.plot(list_of_cs_si - list_of_model_si,list_of_Yinds,'.', alpha = 0.5, label = 'Si')
ax.grid('on')
ax.set_ylabel('Yind')
ax.set_xlabel('Observed - Modelled')
ax.legend();
fig, ax = plt.subplots(4,1, figsize = (12,12))
ax[0].plot(list_of_datetimes[list_of_datetimes < datetime.date(2015,1,1)],
list_of_cs_ni[list_of_datetimes < datetime.date(2015,1,1)]
- list_of_model_ni[list_of_datetimes < datetime.date(2015,1,1)],
'.', alpha = 0.5, label = '2014 N')
ax[0].plot(list_of_datetimes[list_of_datetimes < datetime.date(2015,1,1)],
list_of_cs_si[list_of_datetimes < datetime.date(2015,1,1)]
- list_of_model_si[list_of_datetimes < datetime.date(2015,1,1)],
'.', alpha = 0.5, label = '2014 Si')
ax[0].grid('on')
ax[0].set_ylabel('Observed - Modelled')
ax[0].legend();
ax[1].plot(list_of_datetimes[(list_of_datetimes < datetime.date(2016,1,1))
& (list_of_datetimes > datetime.date(2015,1,1))],
list_of_cs_ni[(list_of_datetimes < datetime.date(2016,1,1))
& (list_of_datetimes > datetime.date(2015,1,1))]
- list_of_model_ni[(list_of_datetimes < datetime.date(2016,1,1))
& (list_of_datetimes > datetime.date(2015,1,1))],
'.', alpha = 0.5, label = '2015 N')
ax[1].plot(list_of_datetimes[(list_of_datetimes < datetime.date(2016,1,1))
& (list_of_datetimes > datetime.date(2015,1,1))],
list_of_cs_si[(list_of_datetimes < datetime.date(2016,1,1))
& (list_of_datetimes > datetime.date(2015,1,1))]
- list_of_model_si[(list_of_datetimes < datetime.date(2016,1,1))
& (list_of_datetimes > datetime.date(2015,1,1))],
'.', alpha = 0.5, label = '2015 Si')
ax[1].grid('on')
ax[1].set_ylabel('Observed - Modelled')
ax[1].legend();
ax[2].plot(list_of_datetimes[(list_of_datetimes < datetime.date(2017,1,1))
& (list_of_datetimes > datetime.date(2016,1,1))],
list_of_cs_ni[(list_of_datetimes < datetime.date(2017,1,1))
& (list_of_datetimes > datetime.date(2016,1,1))]
- list_of_model_ni[(list_of_datetimes < datetime.date(2017,1,1))
& (list_of_datetimes > datetime.date(2016,1,1))],
'.', alpha = 0.5, label = '2016 N')
ax[2].plot(list_of_datetimes[(list_of_datetimes < datetime.date(2017,1,1))
& (list_of_datetimes > datetime.date(2016,1,1))],
list_of_cs_si[(list_of_datetimes < datetime.date(2017,1,1))
& (list_of_datetimes > datetime.date(2016,1,1))]
- list_of_model_si[(list_of_datetimes < datetime.date(2017,1,1))
& (list_of_datetimes > datetime.date(2016,1,1))],
'.', alpha = 0.5, label = '2016 Si')
ax[2].grid('on')
ax[2].set_ylabel('Observed - Modelled')
ax[2].legend();
ax[3].plot(list_of_datetimes[(list_of_datetimes > datetime.date(2017,1,1))],
list_of_cs_ni[(list_of_datetimes > datetime.date(2017,1,1))]
- list_of_model_ni[(list_of_datetimes > datetime.date(2017,1,1))],
'.', alpha = 0.5, label = '2017 N')
ax[3].plot(list_of_datetimes[(list_of_datetimes > datetime.date(2017,1,1))],
list_of_cs_si[ (list_of_datetimes > datetime.date(2017,1,1))]
- list_of_model_si[(list_of_datetimes > datetime.date(2017,1,1))],
'.', alpha = 0.5, label = '2017 Si')
ax[3].grid('on')
ax[3].set_xlabel('Date')
ax[3].set_ylabel('Observed - Modelled')
ax[3].legend();
for a in ax:
a.set_ylim(-5, 35)
f.dropna(subset=['month', 'day', 'depth', 'station', 'chl']).shape
(784, 14)
f3 = f.dropna(subset=['month', 'day', 'depth', 'station', 'chl'])
np.array([datetime.date(2017, f3.month[n], f3.day[n]) for n in f3.index]).min()
datetime.date(2017, 9, 27)
np.array([datetime.date(2017, f3.month[n], f3.day[n]) for n in f3.index]).max()
datetime.date(2017, 11, 19)
f3.year.unique()
array([2008, 2006, 2009, 2007])
list_of_Yinds2 = np.array([])
list_of_Xinds2 = np.array([])
list_of_datetimes2 = np.array([])
list_of_chl = np.array([])
list_of_model_chl = np.array([])
list_of_depths2 = np.array([])
for n in f3.index:
Yind, Xind = stations[f3.station[n]]
depth = np.argmin(np.abs(deptht - f3.depth[n]))
if mesh.variables['tmask'][0,depth,Yind, Xind] == 1:
for year in [2014,2015,2016,2017]:
date = datetime.date(year, f3.month[n], f3.day[n])
sub_dir = date.strftime('%d%b%y').lower()
datestr = date.strftime('%Y%m%d')
fname = 'SalishSea_1d_{}_{}_ptrc_T.nc'.format(datestr, datestr)
nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
chl_val = 1.6*(nuts.variables['diatoms'][0, depth, Yind, Xind]
+ nuts.variables['ciliates'][0, depth, Yind, Xind]
+ nuts.variables['flagellates'][0, depth, Yind, Xind])
list_of_Yinds2 = np.append(list_of_Yinds2, Yind)
list_of_Xinds2 = np.append(list_of_Xinds2, Xind)
list_of_datetimes2 = np.append(list_of_datetimes2, date)
list_of_chl = np.append(list_of_chl, float(f3['chl'][n]))
list_of_model_chl = np.append(list_of_model_chl, chl_val)
list_of_depths2 = np.append(list_of_depths2, depth)
list_of_years2 = np.array([])
for n in f3.index:
Yind, Xind = stations[f3.station[n]]
depth = np.argmin(np.abs(deptht - f3.depth[n]))
if mesh.variables['tmask'][0,depth,Yind, Xind] == 1:
for year in [2014,2015,2016,2017]:
list_of_years2 = np.append(list_of_years2, int(f3.year[n]))
np.unique(list_of_depths2)
array([ 0., 4., 9., 14., 18., 21., 23., 24., 25.])
fig, ax = plt.subplots(figsize = (8,8))
colours = cmocean.cm.dense(np.linspace(0,1,10))
for depth, colour in zip(np.unique(list_of_depths2), colours[1:]):
ax.plot(list_of_chl[list_of_depths2 == depth], list_of_model_chl[list_of_depths2==depth],
'o', color = colour, alpha = 0.5, label = 'depth = '+str(depth))
ax.plot(np.arange(0,6), np.arange(0,6), color = 'grey')
ax.set_title('Chl')
ax.set_xlabel('Observed')
ax.set_ylabel('Modelled')
ax.legend();
np.unique(list_of_years2)
array([2006., 2007., 2008., 2009.])
fig, ax = plt.subplots(figsize = (8,8))
colours = cm.Set1(np.linspace(0,1,9))
for year, colour in zip(np.unique(list_of_years2), colours):
ax.plot(list_of_chl[list_of_years2 == year], list_of_model_chl[list_of_years2==year],
'o', color = colour, alpha = 0.5, label = 'year = '+str(year))
ax.plot(np.arange(0,6), np.arange(0,6), color = 'grey')
ax.set_title('Chl')
ax.set_xlabel('Observed')
ax.set_ylabel('Modelled')
ax.legend();
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(list_of_depths2, list_of_chl - list_of_model_chl,'.', alpha = 0.3)
ax.grid('on')
ax.set_xlabel('Model depth level')
ax.set_ylabel('Observed - Modelled')
ax.set_title('Chl');
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(list_of_chl - list_of_model_chl,list_of_Yinds2,'.', alpha = 0.3, label = 'N')
ax.grid('on')
ax.set_ylabel('Yind')
ax.set_xlabel('Observed - Modelled')
ax.set_title('Chl');
fig, ax = plt.subplots(4,1, figsize = (12,12))
ax[0].plot(list_of_datetimes2[list_of_datetimes2 < datetime.date(2015,1,1)],
list_of_chl[list_of_datetimes2 < datetime.date(2015,1,1)]
- list_of_model_chl[list_of_datetimes2 < datetime.date(2015,1,1)],
'.', alpha = 0.5, label = '2014 Chl')
ax[0].grid('on')
ax[0].set_ylabel('Observed - Modelled')
ax[0].legend();
ax[1].plot(list_of_datetimes2[(list_of_datetimes2 < datetime.date(2016,1,1))
& (list_of_datetimes2 > datetime.date(2015,1,1))],
list_of_chl[(list_of_datetimes2 < datetime.date(2016,1,1))
& (list_of_datetimes2 > datetime.date(2015,1,1))]
- list_of_model_chl[(list_of_datetimes2 < datetime.date(2016,1,1))
& (list_of_datetimes2 > datetime.date(2015,1,1))],
'.', alpha = 0.5, label = '2015 Chl')
ax[1].grid('on')
ax[1].set_ylabel('Observed - Modelled')
ax[1].legend();
ax[2].plot(list_of_datetimes2[(list_of_datetimes2 < datetime.date(2017,1,1))
& (list_of_datetimes2 > datetime.date(2016,1,1))],
list_of_chl[(list_of_datetimes2 < datetime.date(2017,1,1))
& (list_of_datetimes2 > datetime.date(2016,1,1))]
- list_of_model_chl[(list_of_datetimes2 < datetime.date(2017,1,1))
& (list_of_datetimes2 > datetime.date(2016,1,1))],
'.', alpha = 0.5, label = '2016 Chl')
ax[2].grid('on')
ax[2].set_ylabel('Observed - Modelled')
ax[2].legend();
ax[3].plot(list_of_datetimes2[(list_of_datetimes2 > datetime.date(2017,1,1))],
list_of_chl[(list_of_datetimes2 > datetime.date(2017,1,1))]
- list_of_model_chl[(list_of_datetimes2 > datetime.date(2017,1,1))],
'.', alpha = 0.5, label = '2017 Chl')
ax[3].grid('on')
ax[3].set_xlabel('Date')
ax[3].set_ylabel('Observed - Modelled')
ax[3].legend();
for a in ax:
a.set_ylim(-5, 10)
p = pd.read_csv('/ocean/eolson/MEOPAR/obs/NANOOS_SJPEF/ticket2180-csv/Oceanography.csv')
p.head()
/home/vdo/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py:2717: DtypeWarning: Columns (6) have mixed types. Specify dtype option on import or set low_memory=False. interactivity=interactivity, compiler=compiler, result=result)
year | month | day | station | cast | cast_time | source_file | depth | temperature | salinity | ... | o2v | conductivity | pressure | fluorescence | corr_fluorescence | density | descent | nbins | scancnt | flag | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2012 | 9.0 | 28.0 | N | 6.0 | NaN | NaN | 108.0 | 10.02 | 30.69 | ... | NaN | 3.39 | 108.94 | 2.35 | NaN | 23.59 | 0.52 | NaN | NaN | NaN |
1 | 2012 | 9.0 | 28.0 | N | 6.0 | NaN | NaN | 108.5 | 10.02 | 30.69 | ... | NaN | 3.39 | 109.45 | 2.49 | NaN | 23.59 | 0.52 | NaN | NaN | NaN |
2 | 2012 | 9.0 | 28.0 | N | 6.0 | NaN | NaN | 109.0 | 10.02 | 30.69 | ... | NaN | 3.39 | 109.95 | 2.79 | NaN | 23.59 | 0.51 | NaN | NaN | NaN |
3 | 2012 | 9.0 | 28.0 | N | 6.0 | NaN | NaN | 109.5 | 10.03 | 30.69 | ... | NaN | 3.39 | 110.46 | 2.49 | NaN | 23.58 | 0.51 | NaN | NaN | NaN |
4 | 2012 | 9.0 | 28.0 | N | 6.0 | NaN | NaN | 110.0 | 10.03 | 30.69 | ... | NaN | 3.39 | 110.96 | 2.51 | NaN | 23.58 | 0.52 | NaN | NaN | NaN |
5 rows × 24 columns
p = p[p.year >= 2014]
p.keys()
Index(['year', 'month', 'day', 'station', 'cast', 'cast_time', 'source_file', 'depth', 'temperature', 'salinity', 'o2', 'corr_o2', 'o2satper', 'o2satmgl', 'o2v', 'conductivity', 'pressure', 'fluorescence', 'corr_fluorescence', 'density', 'descent', 'nbins', 'scancnt', 'flag'], dtype='object')
p.shape
(17980, 24)
p = p.dropna(subset = ['year', 'day', 'month', 'station', 'depth',
'temperature', 'salinity', 'pressure'])
p.shape
(14708, 24)
p.station.unique()
array(['C', 'B', 'A', 'N', 'S'], dtype=object)
for station in p.station.unique()[:3]:
Yind, Xind = geo_tools.find_closest_model_point(g[g.code == station].lon.values[0],
g[g.code == station].lat.values[0],
X, Y, land_mask = bathy.mask)
stations[station] = (Yind, Xind)
list_of_Yinds3 = np.array([])
list_of_Xinds3 = np.array([])
list_of_datetimes3 = np.array([])
list_of_cs_sal = np.array([])
list_of_cs_tem = np.array([])
list_of_model_sal = np.array([])
list_of_model_tem = np.array([])
list_of_depths3 = np.array([])
list_of_pressure = np.array([])
for n in p.index:
Yind, Xind = stations[p.station[n]]
depth = np.argmin(np.abs(deptht - p.depth[n]))
if mesh.variables['tmask'][0,depth,Yind, Xind] == 1:
date = datetime.date(int(p.year[n]), int(p.month[n]), int(p.day[n]))
sub_dir = date.strftime('%d%b%y').lower()
datestr = date.strftime('%Y%m%d')
fname = 'SalishSea_1d_{}_{}_grid_T.nc'.format(datestr, datestr)
nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
sal_val = nuts.variables['vosaline'][0, depth, Yind, Xind]
tem_val = nuts.variables['votemper'][0, depth, Yind, Xind]
list_of_Yinds3 = np.append(list_of_Yinds3, Yind)
list_of_Xinds3 = np.append(list_of_Xinds3, Xind)
list_of_datetimes3 = np.append(list_of_datetimes3, date)
list_of_cs_sal = np.append(list_of_cs_sal, float(p['salinity'][n]))
list_of_cs_tem = np.append(list_of_cs_tem, float(p['temperature'][n]))
list_of_model_sal = np.append(list_of_model_sal, sal_val)
list_of_model_tem = np.append(list_of_model_tem, tem_val)
list_of_depths3 = np.append(list_of_depths3, depth)
list_of_pressure = np.append(list_of_pressure, p['pressure'][n])
s = gsw_calls.generic_gsw_caller('gsw_SR_from_SP.m', [list_of_cs_sal])
t = gsw_calls.generic_gsw_caller('gsw_CT_from_t.m',
[s, list_of_cs_tem, list_of_pressure])
fig, ax = plt.subplots(figsize = ((10,8)))
viz_tools.plot_coastline(ax, grid)
viz_tools.set_aspect(ax)
ax.set_ylim(250,350)
ax.set_xlim(200,300)
ax.plot(list_of_Xinds3, list_of_Yinds3, 'ro');
from matplotlib.colors import LogNorm
s.shape
(14579,)
fig, ax = plt.subplots(figsize = (10,10))
c, xedge, yedge, im = ax.hist2d(s, list_of_model_sal, bins = 100, norm=LogNorm())
im
ax.plot(np.arange(28,35), np.arange(28,35), color = 'grey')
fig.colorbar(im, ax=ax)
ax.set_xlabel('NANOOS')
ax.set_ylabel('Nowcast-green')
ax.set_title('NANOOS vs Nowcast-green salinity')
print('bias = ' + str(-np.mean(s) + np.mean(list_of_model_sal)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_sal - s)**2)
/ len(list_of_model_sal))))
xbar = np.mean(s)
print('Willmott = ' + str(1-(np.sum((list_of_model_sal - s)**2) /
np.sum((np.abs(list_of_model_sal - xbar)
+ np.abs(s - xbar))**2))))
bias = 0.13353483468716476 RMSE = 0.41932263796326535 Willmott = 0.8433595642923585
fig, ax = plt.subplots(figsize = (10,10))
c, xedge, yedge, im = ax.hist2d(t, list_of_model_tem, bins = 100, norm=LogNorm())
im
ax.plot(np.arange(35), np.arange(35), color = 'grey')
fig.colorbar(im, ax=ax)
ax.set_xlabel('NANOOS')
ax.set_ylabel('Nowcast-green')
ax.set_title('NANOOS vs Nowcast-green temperature')
print('bias = ' + str(-np.mean(t) + np.mean(list_of_model_tem)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_tem - t)**2)
/ len(list_of_model_tem))))
xbar = np.mean(t)
print('Willmott = ' + str(1-(np.sum((list_of_model_tem - t)**2) /
np.sum((np.abs(list_of_model_tem - xbar)
+ np.abs(t - xbar))**2))))
bias = -0.4300265998100876 RMSE = 0.5960042238649761 Willmott = 0.7320908904807624
np.unique(list_of_depths3)
array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27.])
model_depths3 = np.array([deptht[int(n)] for n in list_of_depths3])
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(model_depths3, s - list_of_model_sal,'.', alpha = 0.3)
ax.grid('on')
ax.set_xlabel('Depth (m)')
ax.set_ylabel('Observed - Modelled')
ax.set_title('Salinity');
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(model_depths3, t - list_of_model_tem,'.', color = 'ForestGreen', alpha = 0.3)
ax.grid('on')
ax.set_xlabel('Depth (m)')
ax.set_ylabel('Observed - Modelled')
ax.set_title('Temperature');
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(list_of_Yinds3, s - list_of_model_sal,'.', alpha = 0.3)
ax.grid('on')
ax.set_xlabel('Yind')
ax.set_ylabel('Observed - Modelled')
ax.set_title('Salinity');
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(list_of_Yinds3, t - list_of_model_tem,'.', color = 'ForestGreen', alpha = 0.3)
ax.grid('on')
ax.set_xlabel('Yind')
ax.set_ylabel('Observed - Modelled')
ax.set_title('Temperature');
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(list_of_datetimes3, s - list_of_model_sal,'.', alpha = 0.3)
ax.grid('on')
ax.set_xlabel('Date')
ax.set_ylabel('Observed - Modelled')
ax.set_title('Salinity');
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(list_of_datetimes3, t - list_of_model_tem,'.', color = 'ForestGreen', alpha = 0.3)
ax.grid('on')
ax.set_xlabel('Date')
ax.set_ylabel('Observed - Modelled')
ax.set_title('Temperature');
p = pd.read_csv('/ocean/eolson/MEOPAR/obs/NANOOS_SJPEF/ticket2180-csv/Oceanography.csv')
p.head()
/home/vdo/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py:2717: DtypeWarning: Columns (6) have mixed types. Specify dtype option on import or set low_memory=False. interactivity=interactivity, compiler=compiler, result=result)
year | month | day | station | cast | cast_time | source_file | depth | temperature | salinity | ... | o2v | conductivity | pressure | fluorescence | corr_fluorescence | density | descent | nbins | scancnt | flag | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2012 | 9.0 | 28.0 | N | 6.0 | NaN | NaN | 108.0 | 10.02 | 30.69 | ... | NaN | 3.39 | 108.94 | 2.35 | NaN | 23.59 | 0.52 | NaN | NaN | NaN |
1 | 2012 | 9.0 | 28.0 | N | 6.0 | NaN | NaN | 108.5 | 10.02 | 30.69 | ... | NaN | 3.39 | 109.45 | 2.49 | NaN | 23.59 | 0.52 | NaN | NaN | NaN |
2 | 2012 | 9.0 | 28.0 | N | 6.0 | NaN | NaN | 109.0 | 10.02 | 30.69 | ... | NaN | 3.39 | 109.95 | 2.79 | NaN | 23.59 | 0.51 | NaN | NaN | NaN |
3 | 2012 | 9.0 | 28.0 | N | 6.0 | NaN | NaN | 109.5 | 10.03 | 30.69 | ... | NaN | 3.39 | 110.46 | 2.49 | NaN | 23.58 | 0.51 | NaN | NaN | NaN |
4 | 2012 | 9.0 | 28.0 | N | 6.0 | NaN | NaN | 110.0 | 10.03 | 30.69 | ... | NaN | 3.39 | 110.96 | 2.51 | NaN | 23.58 | 0.52 | NaN | NaN | NaN |
5 rows × 24 columns
p2 = p[p.year>=2006]
p2 = p2[p2.year <=2009]
p2 = p2.dropna(subset = ['year', 'day', 'month', 'station', 'depth',
'temperature', 'salinity', 'pressure'])
p2.shape
(20461, 24)
p2.station.unique()
array(['S', 'N', 'B', 'A', 'C'], dtype=object)
for station in p2.station.unique():
Yind, Xind = geo_tools.find_closest_model_point(g[g.code == station].lon.values[0],
g[g.code == station].lat.values[0],
X, Y, land_mask = bathy.mask)
stations[station] = (Yind, Xind)
stations.keys()
dict_keys(['I', 'O', 'N', 'S', 'T', 'A', 'B', 'C'])
list_of_Yinds4 = np.array([])
list_of_Xinds4 = np.array([])
list_of_datetimes4 = np.array([])
list_of_cs_sal2 = np.array([])
list_of_cs_tem2 = np.array([])
list_of_model_sal2 = np.array([])
list_of_model_tem2 = np.array([])
list_of_depths4 = np.array([])
list_of_pressure2 = np.array([])
list_of_years4 = np.array([])
for n in p2.index:
Yind, Xind = stations[p2.station[n]]
depth = np.argmin(np.abs(deptht - p2.depth[n]))
if mesh.variables['tmask'][0,depth,Yind, Xind] == 1:
for year in [2014, 2015, 2016, 2017]:
date = datetime.date(year, int(p2.month[n]), int(p2.day[n]))
sub_dir = date.strftime('%d%b%y').lower()
datestr = date.strftime('%Y%m%d')
fname = 'SalishSea_1d_{}_{}_grid_T.nc'.format(datestr, datestr)
nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
sal_val = nuts.variables['vosaline'][0, depth, Yind, Xind]
tem_val = nuts.variables['votemper'][0, depth, Yind, Xind]
list_of_Yinds4 = np.append(list_of_Yinds4, Yind)
list_of_Xinds4 = np.append(list_of_Xinds4, Xind)
list_of_datetimes4 = np.append(list_of_datetimes4, date)
list_of_cs_sal2 = np.append(list_of_cs_sal2, float(p2['salinity'][n]))
list_of_cs_tem2 = np.append(list_of_cs_tem2, float(p2['temperature'][n]))
list_of_model_sal2 = np.append(list_of_model_sal2, sal_val)
list_of_model_tem2 = np.append(list_of_model_tem2, tem_val)
list_of_depths4 = np.append(list_of_depths4, depth)
list_of_pressure2 = np.append(list_of_pressure2, p2['pressure'][n])
list_of_years4 = np.append(list_of_years4, int(p2.year[n]))
s2 = gsw_calls.generic_gsw_caller('gsw_SR_from_SP.m', [list_of_cs_sal2])
t2 = gsw_calls.generic_gsw_caller('gsw_CT_from_t.m',
[s2, list_of_cs_tem2, list_of_pressure2])
t2.shape
(80852,)
fig, ax = plt.subplots(figsize = (10,10))
c, xedge, yedge, im = ax.hist2d(s2, list_of_model_sal2, bins = 100,
norm=LogNorm())
im
ax.plot(np.arange(28,35), np.arange(28,35), color = 'grey')
fig.colorbar(im, ax=ax)
ax.set_xlabel('NANOOS')
ax.set_ylabel('Nowcast-green')
ax.set_title('NANOOS vs Nowcast-green salinity')
print('bias = ' + str(-np.mean(s2) + np.mean(list_of_model_sal2)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_sal2 - s2)**2)
/ len(list_of_model_sal2))))
xbar = np.mean(s2)
print('Willmott = ' + str(1-(np.sum((list_of_model_sal2 - s2)**2) /
np.sum((np.abs(list_of_model_sal2 - xbar)
+ np.abs(s2 - xbar))**2))))
bias = -0.25460206199699087 RMSE = 0.5188934979996966 Willmott = 0.7816300601990559
fig, ax = plt.subplots(figsize = (10,10))
c, xedge, yedge, im = ax.hist2d(t2, list_of_model_tem2, bins = 100,
norm=LogNorm())
im
ax.plot(np.arange(35), np.arange(35), color = 'grey')
fig.colorbar(im, ax=ax)
ax.set_xlabel('NANOOS')
ax.set_ylabel('Nowcast-green')
ax.set_title('NANOOS vs Nowcast-green temperature')
print('bias = ' + str(-np.mean(t2) + np.mean(list_of_model_tem2)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_tem2 - t2)**2)
/ len(list_of_model_tem2))))
xbar = np.mean(t2)
print('Willmott = ' + str(1-(np.sum((list_of_model_tem2 - t2)**2) /
np.sum((np.abs(list_of_model_tem2 - xbar)
+ np.abs(t2 - xbar))**2))))
bias = 0.20463601998433134 RMSE = 0.41351953088856835 Willmott = 0.8455840533273167
model_depths4 = np.array([deptht[int(n)] for n in list_of_depths4])
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(model_depths4, s2 - list_of_model_sal2,'.', alpha = 0.3)
ax.grid('on')
ax.set_xlabel('Depth (m)')
ax.set_ylabel('Observed - Modelled')
ax.set_title('Salinity');
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(model_depths4, t2 - list_of_model_tem2,
'.', color = 'ForestGreen', alpha = 0.3)
ax.grid('on')
ax.set_xlabel('Depth (m)')
ax.set_ylabel('Observed - Modelled')
ax.set_title('Temperature');
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(list_of_Yinds4, s2 - list_of_model_sal2,'.', alpha = 0.3)
ax.grid('on')
ax.set_xlabel('Yind')
ax.set_ylabel('Observed - Modelled')
ax.set_title('Salinity');
fig, ax = plt.subplots(figsize = (12,4))
ax.plot(list_of_Yinds4, t2 - list_of_model_tem2,
'.', color = 'ForestGreen', alpha = 0.3)
ax.grid('on')
ax.set_xlabel('Yind')
ax.set_ylabel('Observed - Modelled')
ax.set_title('Temperature');
fig, ax = plt.subplots(4,1, figsize = (12,12))
ax[0].plot(list_of_datetimes4[list_of_datetimes4 < datetime.date(2015,1,1)],
list_of_cs_sal2[list_of_datetimes4 < datetime.date(2015,1,1)]
- list_of_model_sal2[list_of_datetimes4 < datetime.date(2015,1,1)],
'.', alpha = 0.2, label = '2014 salinity')
ax[0].grid('on')
ax[0].set_ylabel('Observed - Modelled')
ax[0].legend();
ax[1].plot(list_of_datetimes4[(list_of_datetimes4 < datetime.date(2016,1,1))
& (list_of_datetimes4 > datetime.date(2015,1,1))],
list_of_cs_sal2[(list_of_datetimes4 < datetime.date(2016,1,1))
& (list_of_datetimes4 > datetime.date(2015,1,1))]
- list_of_model_sal2[(list_of_datetimes4 < datetime.date(2016,1,1))
& (list_of_datetimes4 > datetime.date(2015,1,1))],
'.', alpha = 0.2, label = '2015 salinity')
ax[1].grid('on')
ax[1].set_ylabel('Observed - Modelled')
ax[1].legend();
ax[2].plot(list_of_datetimes4[(list_of_datetimes4 < datetime.date(2017,1,1))
& (list_of_datetimes4 > datetime.date(2016,1,1))],
list_of_cs_sal2[(list_of_datetimes4 < datetime.date(2017,1,1))
& (list_of_datetimes4 > datetime.date(2016,1,1))]
- list_of_model_sal2[(list_of_datetimes4 < datetime.date(2017,1,1))
& (list_of_datetimes4 > datetime.date(2016,1,1))],
'.', alpha = 0.2, label = '2016 salinity')
ax[2].grid('on')
ax[2].set_ylabel('Observed - Modelled')
ax[2].legend();
ax[3].plot(list_of_datetimes4[(list_of_datetimes4 > datetime.date(2017,1,1))],
list_of_cs_sal2[(list_of_datetimes4 > datetime.date(2017,1,1))]
- list_of_model_sal2[(list_of_datetimes4 > datetime.date(2017,1,1))],
'.', alpha = 0.2, label = '2017 salinity')
ax[3].grid('on')
ax[3].set_xlabel('Date')
ax[3].set_ylabel('Observed - Modelled')
ax[3].legend();
for a in ax:
a.set_ylim(-3,4)
fig, ax = plt.subplots(4,1, figsize = (12,12))
ax[0].plot(list_of_datetimes4[list_of_datetimes4 < datetime.date(2015,1,1)],
list_of_cs_tem2[list_of_datetimes4 < datetime.date(2015,1,1)]
- list_of_model_tem2[list_of_datetimes4 < datetime.date(2015,1,1)],
'.', color = 'ForestGreen', alpha = 0.2, label = '2014 temperature')
ax[0].grid('on')
ax[0].set_ylabel('Observed - Modelled')
ax[0].legend();
ax[1].plot(list_of_datetimes4[(list_of_datetimes4 < datetime.date(2016,1,1))
& (list_of_datetimes4 > datetime.date(2015,1,1))],
list_of_cs_tem2[(list_of_datetimes4 < datetime.date(2016,1,1))
& (list_of_datetimes4 > datetime.date(2015,1,1))]
- list_of_model_tem2[(list_of_datetimes4 < datetime.date(2016,1,1))
& (list_of_datetimes4 > datetime.date(2015,1,1))],
'.', color = 'ForestGreen', alpha = 0.2, label = '2015 temperature')
ax[1].grid('on')
ax[1].set_ylabel('Observed - Modelled')
ax[1].legend();
ax[2].plot(list_of_datetimes4[(list_of_datetimes4 < datetime.date(2017,1,1))
& (list_of_datetimes4 > datetime.date(2016,1,1))],
list_of_cs_tem2[(list_of_datetimes4 < datetime.date(2017,1,1))
& (list_of_datetimes4 > datetime.date(2016,1,1))]
- list_of_model_tem2[(list_of_datetimes4 < datetime.date(2017,1,1))
& (list_of_datetimes4 > datetime.date(2016,1,1))],
'.', color = 'ForestGreen', alpha = 0.2, label = '2016 temperature')
ax[2].grid('on')
ax[2].set_ylabel('Observed - Modelled')
ax[2].legend();
ax[3].plot(list_of_datetimes4[(list_of_datetimes4 > datetime.date(2017,1,1))],
list_of_cs_tem2[(list_of_datetimes4 > datetime.date(2017,1,1))]
- list_of_model_tem2[(list_of_datetimes4 > datetime.date(2017,1,1))],
'.', color = 'ForestGreen', alpha = 0.2, label = '2017 temperature')
ax[3].grid('on')
ax[3].set_xlabel('Date')
ax[3].set_ylabel('Observed - Modelled')
ax[3].legend();
for a in ax:
a.set_ylim(-3,2)