In [32]:
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
In [2]:
import h5py
In [3]:
import matplotlib.gridspec as gridspec
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib as mpl
In [51]:
plt.style.use('/ocean/vdo/MEOPAR/biomodelevalpaper/bioModelEvalPaper.mplstyle')
In [5]:
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>''')
Out[5]:
In [6]:
filepath = '/ocean/eolson/MEOPAR/obs/Hakai/SentryShoal/Full_SUNA_Processed.mat'
In [7]:
arrays = {}
f = h5py.File(filepath)
for k, v in f.items():
    arrays[k] = np.array(v)
In [8]:
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,:]
In [9]:
datenumber = arrays['Full_SUNA_Processed'][-1, :]
In [10]:
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))
In [11]:
#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')
In [12]:
geo_tools.find_closest_model_point(-124 - 59/60, 
                                                    49 + 54.4/60, 
                                                    X, Y, 
                                                    land_mask = bathy.mask)
Out[12]:
(703, 146)
In [13]:
h = nc.Dataset('/results/SalishSea/nowcast-green/01feb18/SalishSea_1d_20180201_20180201_grid_T.nc')
In [14]:
HINDCAST_PATH = '/results/SalishSea/nowcast-green/'
In [16]:
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)
In [18]:
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
In [19]:
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);
In [21]:
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
In [22]:
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
In [105]:
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');
In [106]:
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');
In [107]:
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');
In [108]:
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');