import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
import pandas as pd
import datetime
import xarray as xr
from salishsea_tools import tidetools, geo_tools, viz_tools, nc_tools
import pytz
from matplotlib.colors import LogNorm
import os
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
%matplotlib inline
import h5py
import matplotlib.gridspec as gridspec
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib as mpl
plt.style.use('/ocean/vdo/MEOPAR/biomodelevalpaper/bioModelEvalPaper.mplstyle')
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>''')
filepath = '/ocean/eolson/MEOPAR/obs/Hakai/SentryShoal/Full_SUNA_Processed.mat'
arrays = {}
f = h5py.File(filepath)
for k, v in f.items():
arrays[k] = np.array(v)
dates = arrays['Full_SUNA_Processed'][0,:]
times = arrays['Full_SUNA_Processed'][1,:]
temps = arrays['Full_SUNA_Processed'][2,:]
sals = arrays['Full_SUNA_Processed'][3,:]
nitrate = arrays['Full_SUNA_Processed'][5,:]
datenumber = arrays['Full_SUNA_Processed'][-1, :]
base = datetime.datetime(2000, 1, 1)
py_times = np.array([base for i in range(63721)])
for n in range(63721):
py_times[n] = ((datetime.datetime.fromordinal(int(datenumber[n])))
+ datetime.timedelta(days=datenumber[n]%1)
- datetime.timedelta(days = 366))
#load model grid stuff
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')
geo_tools.find_closest_model_point(-124 - 59/60,
49 + 54.4/60,
X, Y,
land_mask = bathy.mask)
(703, 146)
h = nc.Dataset('/results/SalishSea/nowcast-green/01feb18/SalishSea_1d_20180201_20180201_grid_T.nc')
HINDCAST_PATH = '/results/SalishSea/nowcast-green/'
list_of_model_t = np.array([])
list_of_model_t_west = np.array([])
list_of_model_t_east = np.array([])
list_of_model_t_deep = np.array([])
list_of_model_s = np.array([])
list_of_model_s_west = np.array([])
list_of_model_s_east = np.array([])
list_of_model_s_deep = np.array([])
depth = 0
Yind, Xind = [703, 146]
for n in range(63721):
date = py_times[n]
sub_dir = date.strftime('%d%b%y').lower()
datestr = date.strftime('%Y%m%d')
fname = 'SalishSea_1h_{}_{}_grid_T.nc'.format(datestr, datestr)
nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
if date.minute < 30:
before = datetime.datetime(year = date.year, month = date.month, day = date.day,
hour = (date.hour), minute = 30) - datetime.timedelta(hours=1)
if date.minute >= 30:
before = datetime.datetime(year = date.year, month = date.month, day = date.day,
hour = (date.hour), minute = 30)
after = before + datetime.timedelta(hours=1)
sub_dir2 = after.strftime('%d%b%y').lower()
datestr2 = after.strftime('%Y%m%d')
fname2 = 'SalishSea_1h_{}_{}_grid_T.nc'.format(datestr2, datestr2)
nuts2 = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir2, fname2))
delta = (date.minute) / 60
t_val = (delta*(nuts.variables['votemper'][before.hour, depth, Yind, Xind] ) +
(1- delta)*(nuts2.variables['votemper'][after.hour, depth, Yind, Xind] ))
list_of_model_t = np.append(list_of_model_t, t_val)
t_val_west = (delta*(nuts.variables['votemper'][before.hour, depth, Yind, Xind-3] ) +
(1- delta)*(nuts2.variables['votemper'][after.hour, depth, Yind, Xind-3] ))
t_val_east = (delta*(nuts.variables['votemper'][before.hour, depth, Yind, Xind+3] ) +
(1- delta)*(nuts2.variables['votemper'][after.hour, depth, Yind, Xind+3] ))
t_val_deep = (delta*(nuts.variables['votemper'][before.hour, depth+2, Yind, Xind] ) +
(1- delta)*(nuts2.variables['votemper'][after.hour, depth+2, Yind, Xind] ))
list_of_model_t_west = np.append(list_of_model_t_west, t_val_west)
list_of_model_t_east = np.append(list_of_model_t_east, t_val_east)
list_of_model_t_deep = np.append(list_of_model_t_deep, t_val_deep)
s_val = (delta*(nuts.variables['vosaline'][before.hour, depth, Yind, Xind] ) +
(1- delta)*(nuts2.variables['vosaline'][after.hour, depth, Yind, Xind] ))
list_of_model_s = np.append(list_of_model_s, s_val)
s_val_west = (delta*(nuts.variables['vosaline'][before.hour, depth, Yind, Xind-3] ) +
(1- delta)*(nuts2.variables['vosaline'][after.hour, depth, Yind, Xind-3] ))
s_val_east = (delta*(nuts.variables['vosaline'][before.hour, depth, Yind, Xind+3] ) +
(1- delta)*(nuts2.variables['vosaline'][after.hour, depth, Yind, Xind+3] ))
s_val_deep = (delta*(nuts.variables['vosaline'][before.hour, depth+2, Yind, Xind] ) +
(1- delta)*(nuts2.variables['vosaline'][after.hour, depth+2, Yind, Xind] ))
list_of_model_s_west = np.append(list_of_model_s_west, s_val_west)
list_of_model_s_east = np.append(list_of_model_s_east, s_val_east)
list_of_model_s_deep = np.append(list_of_model_s_deep, s_val_deep)
s2 = np.ma.masked_invalid(sals)
list_of_model_s2 = np.ma.masked_array(list_of_model_s, mask = s2.mask)
t2 = np.ma.masked_invalid(temps)
list_of_model_t2 = np.ma.masked_array(list_of_model_t, mask = t2.mask)
print(np.ma.count(t2), np.ma.count(s2))
63721 63721
fig, ax = plt.subplots(figsize = (8, 8))
ax.pcolormesh(grid.variables['Bathymetry'][:])
ax.plot(Xind, Yind, 'r.')
ax.plot(Xind - 3, Yind, 'r.')
ax.plot(Xind + 3, Yind, 'r.')
ax.set_ylim(650, 800)
ax.set_xlim(100, 200)
viz_tools.set_aspect(ax);
fig, ax = plt.subplots(figsize = (8, 8))
c, xedge, yedge, im = ax.hist2d(s2.compressed(), list_of_model_s2.compressed(),
bins = 100, norm=LogNorm());
ax.plot(np.arange(0,30), np.arange(0,30), color = 'red')
ax.set_title('Sentry Shoal Salinity, depth = 1 m')
ax.set_xlabel('Observed')
ax.set_ylabel('Modelled')
divider = make_axes_locatable(ax)
cax1 = divider.append_axes("right", size="5%", pad=0.15)
cbar = plt.colorbar(im, cax = cax1)
print('bias = ' + str(-np.mean(s2) + np.mean(list_of_model_s2)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_s2 - s2)**2) /
np.ma.count(s2))))
xbar = np.mean(s2)
print('Willmott = ' + str(1-(np.sum((list_of_model_s2 - s2)**2) /
np.sum((np.abs(list_of_model_s2 - xbar)
+ np.abs(s2 - xbar))**2))))
bias = -0.28715771259387424 RMSE = 1.1850808771660581 Willmott = 0.8290720829792324
fig, ax = plt.subplots(figsize = (8, 8))
c, xedge, yedge, im = ax.hist2d(t2.compressed(), list_of_model_t2.compressed(),
bins = 100, norm=LogNorm());
ax.plot(np.arange(0,30), np.arange(0,30), color = 'red')
ax.set_title('Sentry Shoal Temperature, depth = 1 m')
ax.set_xlabel('Observed')
ax.set_ylabel('Modelled')
divider = make_axes_locatable(ax)
cax1 = divider.append_axes("right", size="5%", pad=0.15)
cbar = plt.colorbar(im, cax = cax1)
print('bias = ' + str(-np.mean(t2) + np.mean(list_of_model_t2)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_t2 - t2)**2) /
np.ma.count(t2))))
xbar = np.mean(t2)
print('Willmott = ' + str(1-(np.sum((list_of_model_t2 - t2)**2) /
np.sum((np.abs(list_of_model_t2 - xbar)
+ np.abs(t2 - xbar))**2))))
bias = 0.4049012974465036 RMSE = 1.4453098826255502 Willmott = 0.9424725482806958
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, s2, label = 'observed')
ax.plot(py_times, list_of_model_s_deep, label = 'deep')
ax.plot(py_times, list_of_model_s, label = 'closest model point')
ax.plot(py_times, list_of_model_s_east, label = 'east')
ax.plot(py_times, list_of_model_s_west, label = 'west')
ax.set_title('Spring 2015')
ax.set_xlabel('Date')
ax.set_xlim(py_times[0], datetime.datetime(2015,5,28))
ax.legend()
#ax.grid('on')
ax.set_ylim(20, 30)
ax.set_ylabel('Salinity');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, t2, label = 'observed')
ax.plot(py_times, list_of_model_t_deep, label = 'deep')
ax.plot(py_times, list_of_model_t, label = 'closest model point')
ax.plot(py_times, list_of_model_t_east, label = 'east')
ax.plot(py_times, list_of_model_t_west, label = 'west')
ax.set_title('Spring 2015')
ax.set_xlabel('Date')
ax.set_xlim(py_times[0], datetime.datetime(2015,5,28))
ax.legend()
#ax.grid('on')
ax.set_ylim(9, 20)
ax.set_ylabel('Temperature');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, s2, label = 'observed')
ax.plot(py_times, list_of_model_s_deep, label = 'deep')
ax.plot(py_times, list_of_model_s2, label = 'closest model point')
ax.plot(py_times, list_of_model_s_east, label = 'east')
ax.plot(py_times, list_of_model_s_west, label = 'west')
ax.set_title('Summer 2015')
ax.set_xlabel('Date')
ax.set_xlim(datetime.datetime(2015,6,5), datetime.datetime(2015,9,1))
ax.legend()
#ax.grid('on')
ax.set_ylim(20, 30)
ax.set_ylabel('Salinity');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, t2, label = 'observed')
ax.plot(py_times, list_of_model_t_deep, label = 'deep')
ax.plot(py_times, list_of_model_t2, label = 'closest model point')
ax.plot(py_times, list_of_model_t_east, label = 'east')
ax.plot(py_times, list_of_model_t_west, label = 'west')
ax.set_title('Summer 2015')
ax.set_xlabel('Date')
ax.set_xlim(datetime.datetime(2015,6,5), datetime.datetime(2015,9,1))
ax.legend()
#ax.grid('on')
ax.set_ylim(10, 23.5)
ax.set_ylabel('Temperature');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, s2, label = 'observed')
ax.plot(py_times, list_of_model_s_deep, label = 'deep')
ax.plot(py_times, list_of_model_s, label = 'closest model point')
ax.plot(py_times, list_of_model_s_east, label = 'east')
ax.plot(py_times, list_of_model_s_west, label = 'west')
ax.set_title('Fall 2015')
ax.set_xlabel('Date')
ax.set_xlim(datetime.datetime(2015,9,1), py_times[37652])
ax.legend()
ax.set_ylim(25.5, 30)
#ax.grid('on')
ax.set_ylabel('Salinity');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, t2, label = 'observed')
ax.plot(py_times, list_of_model_t_deep, label = 'deep')
ax.plot(py_times, list_of_model_t, label = 'closest model point')
ax.plot(py_times, list_of_model_t_east, label = 'east')
ax.plot(py_times, list_of_model_t_west, label = 'west')
ax.set_title('Fall 2015')
ax.set_xlabel('Date')
ax.set_xlim(datetime.datetime(2015,9,1), py_times[37652])
ax.legend()
ax.set_ylim(10, 17)
#ax.grid('on')
ax.set_ylabel('Temperature');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, s2, label = 'observed')
ax.plot(py_times, list_of_model_s_deep, label = 'deep')
ax.plot(py_times, list_of_model_s, label = 'closest model point')
ax.plot(py_times, list_of_model_s_east, label = 'east')
ax.plot(py_times, list_of_model_s_west, label = 'west')
ax.set_title('Winter 2016')
ax.set_xlabel('Date')
ax.set_xlim(py_times[37653], datetime.datetime(2016,3,1))
ax.legend()
#ax.grid('on')
ax.set_ylim(25, 29.5)
ax.set_ylabel('Salinity');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, t2, label = 'observed')
ax.plot(py_times, list_of_model_t_deep, label = 'deep')
ax.plot(py_times, list_of_model_t, label = 'closest model point')
ax.plot(py_times, list_of_model_t_east, label = 'east')
ax.plot(py_times, list_of_model_t_west, label = 'west')
ax.set_title('Winter 2016')
ax.set_xlabel('Date')
ax.set_xlim(py_times[37653], datetime.datetime(2016,3,1))
ax.legend()
#ax.grid('on')
ax.set_ylim(7, 8.75)
ax.set_ylabel('Temperature');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, s2, label = 'observed')
ax.plot(py_times, list_of_model_s_deep, label = 'deep')
ax.plot(py_times, list_of_model_s, label = 'closest model point')
ax.plot(py_times, list_of_model_s_east, label = 'east')
ax.plot(py_times, list_of_model_s_west, label = 'west')
ax.set_title('Spring 2016')
ax.set_xlabel('Date')
ax.set_xlim(datetime.datetime(2016,3,1), datetime.datetime(2016,6,1))
ax.legend()
#ax.grid('on')
ax.set_ylim(23, 30)
ax.set_ylabel('Salinity');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, t2, label = 'observed')
ax.plot(py_times, list_of_model_t_deep, label = 'deep')
ax.plot(py_times, list_of_model_t, label = 'closest model point')
ax.plot(py_times, list_of_model_t_east, label = 'east')
ax.plot(py_times, list_of_model_t_west, label = 'west')
ax.set_title('Spring 2016')
ax.set_xlabel('Date')
ax.set_xlim(datetime.datetime(2016,3,1), datetime.datetime(2016,6,1))
ax.legend()
#ax.grid('on')
ax.set_ylim(7.5, 18)
ax.set_ylabel('Temperature');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, s2, label = 'observed')
ax.plot(py_times, list_of_model_s_deep, label = 'deep')
ax.plot(py_times, list_of_model_s, label = 'closest model point')
ax.plot(py_times, list_of_model_s_east, label = 'east')
ax.plot(py_times, list_of_model_s_west, label = 'west')
ax.set_title('Summer 2016')
ax.set_xlabel('Date')
ax.set_xlim(datetime.datetime(2016,6,1), datetime.datetime(2016,9,1))
ax.legend()
#ax.grid('on')
ax.set_ylim(20, 30)
ax.set_ylabel('Salinity');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, t2, label = 'observed')
ax.plot(py_times, list_of_model_t_deep, label = 'deep')
ax.plot(py_times, list_of_model_t, label = 'closest model point')
ax.plot(py_times, list_of_model_t_east, label = 'east')
ax.plot(py_times, list_of_model_t_west, label = 'west')
ax.set_title('Summer 2016')
ax.set_xlabel('Date')
ax.set_xlim(datetime.datetime(2016,6,1), datetime.datetime(2016,9,1))
ax.legend()
#ax.grid('on')
ax.set_ylim(10, 22)
ax.set_ylabel('Temperature');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, s2, label = 'observed')
ax.plot(py_times, list_of_model_s_deep, label = 'deep')
ax.plot(py_times, list_of_model_s, label = 'closest model point')
ax.plot(py_times, list_of_model_s_east, label = 'east')
ax.plot(py_times, list_of_model_s_west, label = 'west')
ax.set_title('Fall 2016')
ax.set_xlabel('Date')
ax.set_xlim(datetime.datetime(2016,9,1), datetime.datetime(2016,11,29))
ax.legend()
#ax.grid('on')
ax.set_ylim(25, 30)
ax.set_ylabel('Salinity');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, t2, label = 'observed')
ax.plot(py_times, list_of_model_t_deep, label = 'deep')
ax.plot(py_times, list_of_model_t, label = 'closest model point')
ax.plot(py_times, list_of_model_t_east, label = 'east')
ax.plot(py_times, list_of_model_t_west, label = 'west')
ax.set_title('Fall 2016')
ax.set_xlabel('Date')
ax.set_xlim(datetime.datetime(2016,9,1), datetime.datetime(2016,11,29))
ax.legend()
#ax.grid('on')
ax.set_ylim(8, 18)
ax.set_ylabel('Temperature');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, s2, label = 'observed')
ax.plot(py_times, list_of_model_s_deep, label = 'deep')
ax.plot(py_times, list_of_model_s, label = 'closest model point')
ax.plot(py_times, list_of_model_s_east, label = 'east')
ax.plot(py_times, list_of_model_s_west, label = 'west')
ax.set_title('Spring 2017')
ax.set_xlabel('Date')
ax.set_xlim(py_times[55684], datetime.datetime(2017,6,1))
ax.legend()
#ax.grid('on')
ax.set_ylim(19, 30)
ax.set_ylabel('Salinity');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, t2, label = 'observed')
ax.plot(py_times, list_of_model_t_deep, label = 'deep')
ax.plot(py_times, list_of_model_t, label = 'closest model point')
ax.plot(py_times, list_of_model_t_east, label = 'east')
ax.plot(py_times, list_of_model_t_west, label = 'west')
ax.set_title('Spring 2017')
ax.set_xlabel('Date')
ax.set_xlim(py_times[55684], datetime.datetime(2017,6,1))
ax.legend()
#ax.grid('on')
ax.set_ylim(6, 19)
ax.set_ylabel('Temperature');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, s2, label = 'observed')
ax.plot(py_times, list_of_model_s_deep, label = 'deep')
ax.plot(py_times, list_of_model_s, label = 'closest model point')
ax.plot(py_times, list_of_model_s_east, label = 'east')
ax.plot(py_times, list_of_model_s_west, label = 'west')
ax.set_title('Summer 2017')
ax.set_xlabel('Date')
ax.set_xlim(datetime.datetime(2017,6,1), py_times[-1])
ax.legend()
#ax.grid('on')
ax.set_ylim(19, 29)
ax.set_ylabel('Salinity');
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, t2, label = 'observed')
ax.plot(py_times, list_of_model_t_deep, label = 'deep')
ax.plot(py_times, list_of_model_t, label = 'closest model point')
ax.plot(py_times, list_of_model_t_east, label = 'east')
ax.plot(py_times, list_of_model_t_west, label = 'west')
ax.set_title('Summer 2017')
ax.set_xlabel('Date')
ax.set_xlim(datetime.datetime(2017,6,1), py_times[-1])
ax.legend()
#ax.grid('on')
ax.set_ylim(10, 22.5)
ax.set_ylabel('Temperature');