In [1]:

import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
import scipy.io as sio
import pandas as pd
import datetime
import xarray as xr
from salishsea_tools import tidetools, geo_tools, viz_tools
from matplotlib.colors import LinearSegmentedColormap
%matplotlib inline

from IPython.core.display import display, HTML
display(HTML("<style>.container { width:90% !important; }</style>"))

In [2]:

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[2]:

In [3]:

buoy = sio.loadmat('/ocean/rich/home/metro/ecbuoy/matlab/ECBUOY.mat')

In [7]:

buoy = buoy['ecbuoy']

In [9]:

buoy_dtype = buoy.dtype
buoydata = {n: buoy[n][0, 0] for n in buoy_dtype.names}

In [13]:

buoydata.keys()

Out[13]:

dict_keys(['comment', 'lat', 'lon', 'tz', 'mtime', 'temp', 'specc', 'turb', 'salR', 'vel', 'airtemp', 'cwind', 'wspd', 'wdir', 'wu', 'wv', 'pH', 'O2', 'O2sat', 'dailyturb'])

In [45]:

temps = buoydata['temp'][:,0]
sals = buoydata['salR'][:,0]
mtimes = buoydata['mtime'][:,0]

In [46]:

mtimes.shape

Out[46]:

(76535,)

In [47]:

base = datetime.datetime(2000, 1, 1)
py_buoy_times = np.array([base for i in range(76535)])
print(py_buoy_times.shape)
for n in range(76535):
    py_buoy_times[n] = ((datetime.datetime.fromordinal(int(mtimes[n])))
                       + datetime.timedelta(days=mtimes[n]%1)
                       - datetime.timedelta(days = 366))

(76535,)

In [30]:

nowcast = xr.open_dataset('https://salishsea.eos.ubc.ca/erddap/griddap/ubcSSg3DTracerFields1hV17-02')

In [32]:

grid = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
bathy, X, Y = tidetools.get_bathy_data(grid)

In [35]:

buoydata['lon'][0][0]

Out[35]:

-123.03856589999999

In [37]:

geo_tools.find_closest_model_point(buoydata['lon'][0][0], buoydata['lat'][0][0], 
                                   X, Y, land_mask = bathy.mask)

Out[37]:

(415, 337)

In [48]:

temps = temps[py_buoy_times > datetime.datetime(2014, 9, 12)]
sals = sals[py_buoy_times > datetime.datetime(2014, 9, 12)]
py_buoy_times = py_buoy_times[py_buoy_times > datetime.datetime(2014, 9, 12)]

In [49]:

print(temps.shape)

(26462,)

In [56]:

sal_surface_model = np.array([])
t_surface_model = np.array([])
sal_1_5_model = np.array([])
t_1_5_model = np.array([])
for n in range(26462):
    date = py_buoy_times[n]
    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)
        delta = (date - before).seconds / 3600
        s_val0 = (delta * ((nowcast.salinity
               .sel(time = before, depth = 0.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values) + 
               (1- delta)*((nowcast.salinity
               .sel(time = before + datetime.timedelta(hours=1), 
                    depth = 0.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values))
        t_val0 = (delta * ((nowcast.temperature
               .sel(time = before, depth = 0.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values) + 
               (1- delta)*((nowcast.temperature
               .sel(time = before + datetime.timedelta(hours=1), 
                    depth = 0.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values))
        s_val1 = (delta * ((nowcast.salinity
               .sel(time = before, depth = 1.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values) + 
               (1- delta)*((nowcast.salinity
               .sel(time = before + datetime.timedelta(hours=1), 
                    depth = 0.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values))
        t_val1 = (delta * ((nowcast.temperature
               .sel(time = before, depth = 1.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values) + 
               (1- delta)*((nowcast.temperature
               .sel(time = before + datetime.timedelta(hours=1), 
                    depth = 0.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values))
    if date.minute >= 30:
        before = datetime.datetime(year = date.year, month = date.month, day = date.day, 
                               hour = (date.hour), minute = 30)
        delta = (date - before).seconds / 3600
        s_val0 = (delta * ((nowcast.salinity
               .sel(time = before, depth = 0.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values) + 
               (1- delta)*((nowcast.salinity
               .sel(time = before + datetime.timedelta(hours=1), 
                    depth = 0.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values))
        t_val0 = (delta * ((nowcast.temperature
               .sel(time = before, depth = 0.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values) + 
               (1- delta)*((nowcast.temperature
               .sel(time = before + datetime.timedelta(hours=1), 
                    depth = 0.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values))
        s_val1 = (delta * ((nowcast.salinity
               .sel(time = before, depth = 1.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values) + 
               (1- delta)*((nowcast.salinity
               .sel(time = before + datetime.timedelta(hours=1), 
                    depth = 0.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values))
        t_val1 = (delta * ((nowcast.temperature
               .sel(time = before, depth = 1.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values) + 
               (1- delta)*((nowcast.temperature
               .sel(time = before + datetime.timedelta(hours=1), 
                    depth = 0.5, method = 'nearest')
               .isel(gridY = 415, gridX = 337)).values))
    sal_surface_model = np.append(sal_surface_model, s_val0)
    t_surface_model = np.append(t_surface_model, t_val0)
    sal_1_5_model = np.append(sal_1_5_model, s_val1)
    t_1_5_model = np.append(t_1_5_model, t_val1)

In [70]:

fig, ax = plt.subplots(figsize = (24,6))
ax.plot(py_buoy_times, temps, 'r-', alpha = 0.5, label = 'buoy')
ax.plot(py_buoy_times, t_1_5_model, 'g-', alpha = 0.5, linewidth = 7, label = 'model - 1.5m')
ax.plot(py_buoy_times, t_surface_model, 'b-', alpha = 0.5, label = 'model - surface')
ax.legend()
ax.set_title('Temperature');

In [74]:

temps = np.ma.masked_invalid(temps)
t_1_5_model = np.ma.masked_array(t_1_5_model, mask = temps.mask)
t_surface_model = np.ma.masked_array(t_surface_model, mask = temps.mask)
print(np.ma.count(temps))

In [76]:

from matplotlib.colors import LogNorm

In [78]:

temps = temps.compressed()
t_1_5_model = t_1_5_model.compressed()
t_surface_model = t_surface_model.compressed()

In [80]:

fig, ax = plt.subplots(figsize = (10,10))
c, xedge, yedge, im = ax.hist2d(temps, t_surface_model,  
                                bins = 100, norm=LogNorm())
fig.colorbar(im, ax=ax)
ax.plot(np.arange(0,22), 'r-')
ax.grid('on')
ax.set_title('Water Quality Buoy vs Nowcast-green: Temperature, model at surface')
ax.set_xlabel('Buoy')
ax.set_ylabel('Nowcast-green');
print('bias =  ' + str(-np.mean(temps) + np.mean(t_surface_model)))
print('RMSE = ' + str(np.sqrt(np.sum((t_surface_model - temps)**2) /
                              len(t_surface_model))))
xbar = np.mean(temps)
print('Willmott = ' + str(1-(np.sum((t_surface_model - temps)**2)  / 
                             np.sum((np.abs(t_surface_model - xbar) 
                                     + np.abs(temps - xbar))**2))))

bias =  -0.834625751208
RMSE = 1.48700936448
Willmott = 0.983537985575

In [82]:

fig, ax = plt.subplots(figsize = (10,10))
c, xedge, yedge, im = ax.hist2d(temps, t_1_5_model,  
                                bins = 100, norm=LogNorm())
fig.colorbar(im, ax=ax)
ax.plot(np.arange(0,22), 'r-')
ax.grid('on')
ax.set_title('Water Quality Buoy vs Nowcast-green: Temperature, model at 1.5m')
ax.set_xlabel('Buoy')
ax.set_ylabel('Nowcast-green');
print('bias =  ' + str(-np.mean(temps) + np.mean(t_1_5_model)))
print('RMSE = ' + str(np.sqrt(np.sum((t_1_5_model - temps)**2) /
                              len(t_1_5_model))))
xbar = np.mean(temps)
print('Willmott = ' + str(1-(np.sum((t_1_5_model - temps)**2)  / 
                             np.sum((np.abs(t_1_5_model - xbar) 
                                     + np.abs(temps - xbar))**2))))

bias =  -0.837083078946
RMSE = 1.48756389894
Willmott = 0.983513130584

In [83]:

fig, ax = plt.subplots(figsize = (24,6))
ax.plot(py_buoy_times, sals, 'r-', alpha = 0.5, label = 'buoy')
ax.plot(py_buoy_times, sal_1_5_model, 'g-', alpha = 0.5, linewidth = 7, label = 'model - 1.5m')
ax.plot(py_buoy_times, sal_surface_model, 'b-', alpha = 0.5, label = 'model - surface')
ax.legend()
ax.set_title('Salinity');

In [89]:

sals = np.ma.masked_invalid(sals)
sal_1_5_model = np.ma.masked_array(sal_1_5_model, mask = sals.mask)
sal_surface_model = np.ma.masked_array(sal_surface_model, mask = sals.mask)
sals = sals.compressed()
sal_1_5_model = sal_1_5_model.compressed()
sal_surface_model = sal_surface_model.compressed()

In [90]:

len(sals)

Out[90]:

In [91]:

fig, ax = plt.subplots(figsize = (10,10))
c, xedge, yedge, im = ax.hist2d(sals, sal_surface_model,  
                                bins = 100, norm=LogNorm())
fig.colorbar(im, ax=ax)
ax.plot(np.arange(0,22), 'r-')
ax.grid('on')
ax.set_title('Water Quality Buoy vs Nowcast-green: Salinity, model at surface')
ax.set_xlabel('Buoy')
ax.set_ylabel('Nowcast-green');
print('bias =  ' + str(-np.mean(sals) + np.mean(sal_surface_model)))
print('RMSE = ' + str(np.sqrt(np.sum((sal_surface_model - sals)**2) /
                              len(sal_surface_model))))
xbar = np.mean(sals)
print('Willmott = ' + str(1-(np.sum((sal_surface_model - sals)**2)  / 
                             np.sum((np.abs(sal_surface_model - xbar) 
                                     + np.abs(sals - xbar))**2))))

bias =  -0.22612423275
RMSE = 0.648830427331
Willmott = 0.277567562816

In [94]:

fig, ax = plt.subplots(figsize = (10,10))
c, xedge, yedge, im = ax.hist2d(sals, sal_1_5_model,  
                                bins = 100, norm=LogNorm())
fig.colorbar(im, ax=ax)
ax.plot(np.arange(0,22), 'r-')
ax.grid('on')
ax.set_title('Water Quality Buoy vs Nowcast-green: Salinity, model at 1.5m')
ax.set_xlabel('Buoy')
ax.set_ylabel('Nowcast-green');
print('bias =  ' + str(-np.mean(sals) + np.mean(sal_1_5_model)))
print('RMSE = ' + str(np.sqrt(np.sum((sal_1_5_model - sals)**2) /
                              len(sal_1_5_model))))
xbar = np.mean(sals)
print('Willmott = ' + str(1-(np.sum((sal_1_5_model - sals)**2)  / 
                             np.sum((np.abs(sal_1_5_model - xbar) 
                                     + np.abs(sals - xbar))**2))))

bias =  -0.22296377628
RMSE = 0.64780516164
Willmott = 0.285052396749

In [ ]: