In [1]:

import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from salishsea_tools import (tidetools, geo_tools, viz_tools)
import numpy.ma as ma
import pandas as pd
import datetime
import pytz
import os
%matplotlib inline

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

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

In [3]:

mesh = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/mesh_mask201702.nc')

In [4]:

HINDCAST_PATH= '/results/SalishSea/nowcast-green/'

In [5]:

f = pd.read_excel('/ocean/vdo/MEOPAR/2016_Nutrients_20180110_CN_edits.xlsx')
f.keys()

Out[5]:

Index(['Crew', 'Date', 'Time of Sample', 'Lat_reported', 'Long_reported',
       'Latitude', 'Longitude', 'Station', 'Depth', 'NO3+NO', 'PO4',
       'Samples missing!', 'Unnamed: 12', 'Samples in lab.'],
      dtype='object')

In [6]:

f = f.drop(f.keys()[11:], axis=1)

In [7]:

f.shape

Out[7]:

(1553, 11)

In [8]:

f = f.dropna(subset = ['Date', 'Time of Sample', 'Latitude', 'Longitude', 'Depth', 'NO3+NO'])

In [9]:

f[:5]

Out[9]:

	Crew	Date	Time of Sample	Lat_reported	Long_reported	Latitude	Longitude	Station	Depth	NO3+NO	PO4
0	Sentry Shoal	2016-01-22	17:30:00	49° 92.000	125° 00.000	50.533333	-125.000	SS-46131	0	18.491	1.951
1	Sentry Shoal	2016-01-22	17:30:00	49° 92.000	125° 00.000	50.533333	-125.000	SS-46131	0	18.662	1.863
2	Sentry Shoal	2016-01-22	17:30:00	49° 92.000	125° 00.000	50.533333	-125.000	SS-46131	0	19.621	1.561
7	Campbell River	2016-02-16	09:58:00	50° 04.900	125° 15.900	50.081667	-125.265	CR-2	0	24.163	2.237
8	Campbell River	2016-02-16	09:58:00	50° 04.900	125° 15.900	50.081667	-125.265	CR-2	20	24.326	2.787

In [10]:

local = pytz.timezone ("America/Los_Angeles")

In [11]:

import datetime

In [12]:

datetimes = np.array([])
for index in f.index:
    dt = datetime.datetime.combine(pd.to_datetime(pd.Timestamp(f['Date'][index])), 
                          f['Time of Sample'][index])
    datetimes = np.append(datetimes, dt)

In [13]:

f.shape

Out[13]:

(1374, 11)

In [14]:

datetimes.shape

Out[14]:

(1374,)

In [15]:

f = f.assign(datetime = datetimes)

In [16]:

f.shape

Out[16]:

(1374, 12)

In [17]:

f.Crew.unique()

Out[17]:

array(['Sentry Shoal', 'Campbell River', 'Lund', 'Powell River',
       'Steveston', 'Nanaimo/Qualicum', 'Baynes Sound', 'Cowichan Bay',
       'Galiano Island', 'Ladysmith', 'Irvines-Sechelt', 'Malaspina Strait'], dtype=object)

In [63]:

list_of_lons = np.array([])
list_of_lats = np.array([])
list_of_datetimes = np.array([])
list_of_cs_ni = np.array([])
list_of_model_ni = np.array([])
list_of_depths = np.array([])
cb_lons = np.array([])
cb_lats = np.array([])
cb_datetimes = np.array([])
cb_cs_ni = np.array([])
cb_model_ni = np.array([])
cb_depths = np.array([])
for n in f.index:
    Yind, Xind = geo_tools.find_closest_model_point(f.Longitude[n], 
                                                        f.Latitude[n], 
                                                        X, Y, land_mask = bathy.mask)
    if f['Depth'][n] == 0:
        depth = 0
    elif f['Depth'][n] == 20:
        depth = 18
    if mesh.variables['tmask'][0,depth,Yind, Xind] == 1:
        date = local.localize(f['datetime'][n], is_dst=True).astimezone(pytz.utc)
        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)
            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 + 30) / 60
            ni_val = ((1-delta)*(nuts.variables['nitrate'][before.hour, depth, Yind, Xind] ) + 
                       (delta)*(nuts2.variables['nitrate'][after.hour, depth, Yind, Xind] ))
        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_lons = np.append(list_of_lons, f.Longitude[n])
        list_of_lats = np.append(list_of_lats, f.Latitude[n])
        list_of_datetimes = np.append(list_of_datetimes, date)
        if f['NO3+NO'][n] == '<0':
            list_of_cs_ni = np.append(list_of_cs_ni, 0)
        else:
            list_of_cs_ni = np.append(list_of_cs_ni, float(f['NO3+NO'][n]))
        list_of_model_ni = np.append(list_of_model_ni, ni_val)
        list_of_depths = np.append(list_of_depths, depth)
        if f.Crew[n] == 'Cowichan Bay':
            cb_lons = np.append(cb_lons, f.Longitude[n])
            cb_lats = np.append(cb_lats, f.Latitude[n])
            cb_depths = np.append(cb_depths, depth)
            cb_datetimes = np.append(cb_datetimes, date)
            if f['NO3+NO'][n] == '<0':
                cb_cs_ni = np.append(cb_cs_ni, 0)
            else:
                cb_cs_ni = np.append(cb_cs_ni, float(f['NO3+NO'][n]))
            cb_model_ni = np.append(cb_model_ni, ni_val)

In [64]:

list_of_lats.shape

Out[64]:

(1363,)

In [65]:

list_of_model_ni.shape

Out[65]:

(1363,)

In [66]:

cb_cs_ni.shape

Out[66]:

(143,)

In [24]:

import matplotlib as mpl
mpl.rcParams['font.size'] = 12
mpl.rcParams['axes.titlesize'] = 12

In [95]:

fig, ax = plt.subplots(figsize = (10,10))
viz_tools.set_aspect(ax, coords = 'map')
ax.plot(list_of_lons, list_of_lats, 'ro')
viz_tools.plot_coastline(ax, grid, coords = 'map')
ax.set_ylim(48.5, 50.8)
ax.set_xlim(-125.5, -122.5);

In [26]:

list_of_cs_ni.shape

Out[26]:

(1363,)

In [89]:

fig, ax = plt.subplots(figsize = (10,10))
ax.plot(list_of_cs_ni[list_of_depths == 18], list_of_model_ni[list_of_depths == 18],  
        'b.', alpha = 0.5, label = '20m')
ax.plot(cb_cs_ni[cb_depths == 18], cb_model_ni[cb_depths == 18],  
        '.', color = 'purple', alpha = 0.5, label = '20m CB')
ax.plot(list_of_cs_ni[list_of_depths == 0], list_of_model_ni[list_of_depths==0],  
        'g.', alpha = 0.5, label = 'surface')
ax.plot(cb_cs_ni[cb_depths == 0], cb_model_ni[cb_depths == 0],  
        '.', color = 'orange', alpha = 0.5, label = 'surface CB')
ax.plot(np.arange(0,30), color = 'grey')
ax.grid('on')
ax.set_title('Citizen Science Nitrate 2016')
ax.set_xlabel('Citizen Science')
ax.set_ylabel('Nowcast-green');
ax.legend()
print('surface bias =  ' + str(-np.mean(list_of_cs_ni[list_of_depths == 0]) 
                               + np.mean(list_of_model_ni[list_of_depths == 0])))
print('surface RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni[list_of_depths == 0] 
                                              - list_of_cs_ni[list_of_depths == 0])**2) /
                              len(list_of_cs_ni[list_of_depths == 0]))))
xbar = np.mean(list_of_cs_ni[list_of_depths == 0])
print('surface Willmott = ' + str(1-(np.sum((list_of_model_ni[list_of_depths == 0] 
                                             - list_of_cs_ni[list_of_depths == 0])**2)  / 
                             np.sum((np.abs(list_of_model_ni[list_of_depths == 0] - xbar) 
                                     + np.abs(list_of_cs_ni[list_of_depths == 0] - xbar))**2))))

print('20m bias =  ' + str(-np.mean(list_of_cs_ni[list_of_depths == 18]) 
                               + np.mean(list_of_model_ni[list_of_depths == 18])))
print('20m RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni[list_of_depths == 18] 
                                              - list_of_cs_ni[list_of_depths == 18])**2) /
                              len(list_of_cs_ni[list_of_depths == 18]))))
xbar = np.mean(list_of_cs_ni[list_of_depths == 18])
print('20m Willmott = ' + str(1-(np.sum((list_of_model_ni[list_of_depths == 18] 
                                             - list_of_cs_ni[list_of_depths == 18])**2)  / 
                             np.sum((np.abs(list_of_model_ni[list_of_depths == 18] - xbar) 
                                     + np.abs(list_of_cs_ni[list_of_depths == 18] - xbar))**2))))

print('bias =  ' + str(-np.mean(list_of_cs_ni) + np.mean(list_of_model_ni)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni - list_of_cs_ni)**2) /
                              len(list_of_cs_ni))))
xbar = np.mean(list_of_cs_ni)
print('Willmott = ' + str(1-(np.sum((list_of_model_ni - list_of_cs_ni)**2)  / 
                             np.sum((np.abs(list_of_model_ni - xbar) 
                                     + np.abs(list_of_cs_ni - xbar))**2))))

surface bias =  -0.714972513791
surface RMSE = 6.64487594192
surface Willmott = 0.790652334117
20m bias =  0.101739293839
20m RMSE = 8.97831572282
20m Willmott = 0.531037891276
bias =  -0.315305033461
RMSE = 7.87365463663
Willmott = 0.785912376897

In [68]:

fig, ax = plt.subplots(figsize = (20,8))
ax.plot(list_of_lats[list_of_depths == 18], 
        list_of_model_ni[list_of_depths == 18] 
        - list_of_cs_ni[list_of_depths == 18], 'bo', alpha =0.5, label = '20m')
ax.plot(list_of_lats[list_of_depths == 0], 
        list_of_model_ni[list_of_depths == 0] 
        - list_of_cs_ni[list_of_depths == 0], 'go', alpha =0.5, label = 'surface')
ax.plot(cb_lats[cb_depths == 18], 
        cb_model_ni[cb_depths == 18] 
        - cb_cs_ni[cb_depths == 18], 'o', color = 'purple', alpha =0.5, label = '20m CB')
ax.plot(cb_lats[cb_depths == 0], 
        cb_model_ni[cb_depths == 0] 
        - cb_cs_ni[cb_depths == 0], 'o', color = 'orange', alpha =0.5, label = 'surface CB')
ax.legend()
ax.grid('on')
ax.set_xlabel('lat', fontsize = 15)
ax.set_ylabel('Model - Observed',fontsize = 15)
ax.set_title('Ni', fontsize = 20);

In [69]:

fig, ax = plt.subplots(figsize = (20,8))
ax.plot(list_of_lons[list_of_depths == 18], 
        list_of_model_ni[list_of_depths == 18] 
        - list_of_cs_ni[list_of_depths == 18], 'bo', alpha =0.5, label = '20m')
ax.plot(list_of_lons[list_of_depths == 0], 
        list_of_model_ni[list_of_depths == 0] 
        - list_of_cs_ni[list_of_depths == 0], 'go', alpha =0.5, label = 'surface')
ax.plot(cb_lons[cb_depths == 18], 
        cb_model_ni[cb_depths == 18] 
        - cb_cs_ni[cb_depths == 18], 'o', color = 'purple', alpha =0.5, label = '20m CB')
ax.plot(cb_lons[cb_depths == 0], 
        cb_model_ni[cb_depths == 0] 
        - cb_cs_ni[cb_depths == 0], 'o', color = 'orange', alpha =0.5, label = 'surface CB')
ax.legend()
ax.grid('on')
ax.set_xlabel('lon', fontsize = 15)
ax.set_ylabel('Model - Observed',fontsize = 15)
ax.set_title('Ni', fontsize = 20);

In [70]:

fig, ax = plt.subplots(figsize = (20,8))
ax.plot(list_of_datetimes[list_of_depths == 18], 
        list_of_model_ni[list_of_depths == 18] 
        - list_of_cs_ni[list_of_depths == 18], 'bo', alpha =0.5, label = '20m')
ax.plot(list_of_datetimes[list_of_depths == 0], 
        list_of_model_ni[list_of_depths == 0] 
        - list_of_cs_ni[list_of_depths == 0], 'go', alpha =0.5, label = 'surface')
ax.plot(cb_datetimes[cb_depths == 18], 
        cb_model_ni[cb_depths == 18] 
        - cb_cs_ni[cb_depths == 18], 'o', color = 'purple', alpha =0.5, label = '20m CB')
ax.plot(cb_datetimes[cb_depths == 0], 
        cb_model_ni[cb_depths == 0] 
        - cb_cs_ni[cb_depths == 0], 'o', color = 'orange', alpha =0.5, label = 'surface CB')
ax.legend()
ax.grid('on')
ax.set_xlabel('date', fontsize = 15)
ax.set_ylabel('Model - Observed',fontsize = 15)
ax.set_title('Ni', fontsize = 20);