In [48]:
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 [3]:
import h5py
In [4]:
import matplotlib.gridspec as gridspec
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib as mpl
In [5]:
plt.style.use('/ocean/vdo/MEOPAR/biomodelevalpaper/bioModelEvalPaper.mplstyle')
In [6]:
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[6]:
In [7]:
filepath = '/ocean/eolson/MEOPAR/obs/Hakai/SentryShoal/Full_SUNA_Processed.mat'
In [8]:
arrays = {}
f = h5py.File(filepath)
for k, v in f.items():
arrays[k] = np.array(v)
In [9]:
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 [10]:
datenumber = arrays['Full_SUNA_Processed'][-1, :]
In [11]:
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 [12]:
#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 [13]:
geo_tools.find_closest_model_point(-124 - 59/60,
49 + 54.4/60,
X, Y,
land_mask = bathy.mask)
Out[13]:
(703, 146)
In [14]:
h = nc.Dataset('/results/SalishSea/nowcast-green/01feb18/SalishSea_1d_20180201_20180201_grid_T.nc')
In [15]:
HINDCAST_PATH = '/results/SalishSea/nowcast-green/'
In [16]:
list_of_model_ni = np.array([])
list_of_model_ni_west = np.array([])
list_of_model_ni_east = np.array([])
list_of_model_ni_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_{}_{}_ptrc_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_{}_{}_ptrc_T.nc'.format(datestr2, datestr2)
nuts2 = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir2, fname2))
delta = (date.minute) / 60
ni_val = (delta*(nuts.variables['nitrate'][before.hour, depth, Yind, Xind] ) +
(1- delta)*(nuts2.variables['nitrate'][after.hour, depth, Yind, Xind] ))
list_of_model_ni = np.append(list_of_model_ni, ni_val)
ni_val_west = (delta*(nuts.variables['nitrate'][before.hour, depth, Yind, Xind-3] ) +
(1- delta)*(nuts2.variables['nitrate'][after.hour, depth, Yind, Xind-3] ))
ni_val_east = (delta*(nuts.variables['nitrate'][before.hour, depth, Yind, Xind+3] ) +
(1- delta)*(nuts2.variables['nitrate'][after.hour, depth, Yind, Xind+3] ))
ni_val_deep = (delta*(nuts.variables['nitrate'][before.hour, depth+2, Yind, Xind] ) +
(1- delta)*(nuts2.variables['nitrate'][after.hour, depth+2, Yind, Xind] ))
list_of_model_ni_west = np.append(list_of_model_ni_west, ni_val_west)
list_of_model_ni_east = np.append(list_of_model_ni_east, ni_val_east)
list_of_model_ni_deep = np.append(list_of_model_ni_deep, ni_val_deep)
In [17]:
list_of_model_ni.shape
Out[17]:
(63721,)
In [18]:
nitrate2 = np.ma.masked_invalid(nitrate)
list_of_model_ni2 = np.ma.masked_array(list_of_model_ni, mask = nitrate2.mask)
In [19]:
list_of_model_ni_deep2 = np.ma.masked_array(list_of_model_ni_deep, mask = nitrate2.mask)
list_of_model_ni_west2 = np.ma.masked_array(list_of_model_ni_west, mask = nitrate2.mask)
list_of_model_ni_east2 = np.ma.masked_array(list_of_model_ni_east, mask = nitrate2.mask)
In [20]:
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]:
np.ma.count(nitrate2)
Out[21]:
60951
In [22]:
fig, ax = plt.subplots(figsize = (8, 8))
c, xedge, yedge, im = ax.hist2d(nitrate2.compressed(), list_of_model_ni2.compressed(),
bins = 100, norm=LogNorm());
ax.plot(np.arange(0,30), np.arange(0,30), color = 'red')
ax.set_title('Sentry Shoal Nitrate, 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(nitrate2) + np.mean(list_of_model_ni2)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni2 - nitrate2)**2) /
np.ma.count(nitrate2))))
xbar = np.mean(nitrate2)
print('Willmott = ' + str(1-(np.sum((list_of_model_ni2 - nitrate2)**2) /
np.sum((np.abs(list_of_model_ni2 - xbar)
+ np.abs(nitrate2 - xbar))**2))))
bias = -1.8818928507550012 RMSE = 5.801877196234966 Willmott = 0.8469078330784996
In [49]:
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, nitrate2, label = 'observed')
ax.plot(py_times, list_of_model_ni_deep, label = 'deep')
ax.plot(py_times, list_of_model_ni, label = 'closest model point')
ax.plot(py_times, list_of_model_ni_east, label = 'east')
ax.plot(py_times, list_of_model_ni_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.set_ylim(-1, 24)
ax.set_ylabel('N');
In [50]:
fig, ax = plt.subplots(figsize = (24, 6))
ax.plot(py_times, nitrate2, label = 'observed')
ax.plot(py_times, list_of_model_ni_deep, label = 'deep')
ax.plot(py_times, list_of_model_ni, label = 'closest model point')
ax.plot(py_times, list_of_model_ni_east, label = 'east')
ax.plot(py_times, list_of_model_ni_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.set_ylim(-1, 23)
ax.set_ylabel('N');