In [1]:
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
from salishsea_tools import geo_tools, nc_tools, tidetools, viz_tools
import xarray as xr
import datetime
import pandas as pd
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:90% !important; }</style>"))

%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]:
ferry = pd.read_csv('https://salishsea.eos.ubc.ca/erddap/tabledap/ubcONCTWDP1mV18-01.csv?time%2Clongitude%2Clatitude%2Cchlorophyll&time%3E=2017-01-01T00%3A00%3A00Z&time%3C=2017-12-31T23%3A59%3A00Z&longitude%3E=-123.78116&longitude%3C=-123.694065')
ferry = ferry.drop(ferry.index[0])
In [4]:
ferry.shape
Out[4]:
(22878, 4)
In [5]:
ferry[:3]
Out[5]:
time longitude latitude chlorophyll
1 2017-02-15T00:36:00Z -123.69407221833332 49.19252778166666 NaN
2 2017-02-15T00:37:00Z -123.70168333 49.195111114999996 NaN
3 2017-02-15T00:38:00Z -123.70926666500002 49.19772222500001 NaN
In [6]:
ferry[-3:]
Out[6]:
time longitude latitude chlorophyll
22876 2017-12-31T21:16:00Z -123.71265555166667 49.18855556 NaN
22877 2017-12-31T21:17:00Z -123.70532221833334 49.18605555833333 NaN
22878 2017-12-31T21:18:00Z -123.69793333000001 49.18355556 NaN
In [7]:
import pickle
In [8]:
HINDCAST_PATH= '/results/SalishSea/nowcast-green/'
In [9]:
import os
In [10]:
ferry = ferry.dropna()
In [11]:
ferry.shape
Out[11]:
(19836, 4)
In [12]:
list_of_ferry_chl = np.array([])
list_of_date = np.array([])
for n in ferry.index:
    date = datetime.datetime.strptime(ferry.time[n][:-1], '%Y-%m-%dT%H:%M:%S')
    list_of_ferry_chl = np.append(list_of_ferry_chl, 
                                  float(ferry.chlorophyll[n]))
    list_of_date = np.append(list_of_date, date)
In [13]:
mask = np.zeros((898,398))
for n in ferry.index[:1000]:
    Yind, Xind = geo_tools.find_closest_model_point(float(ferry.longitude[n]), 
                                                    float(ferry.latitude[n]), 
                                                    X, Y, 
                                                    land_mask = bathy.mask)
    mask[Yind, Xind] = 1
In [14]:
t0=datetime.datetime(2017,1,1)
model_average = np.array([])
for date in [t0+datetime.timedelta(days=ii*1) for ii in range(365)]:
    sub_dir = date.strftime('%d%b%y').lower()
    datestr = date.strftime('%Y%m%d')
    fname = 'SalishSea_1d_{}_{}_ptrc_T.nc'.format(datestr, datestr)
    nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
    chl = (1.6 * np.ma.masked_equal((nuts.variables['diatoms'][0, 1, ...] 
               + nuts.variables['ciliates'][0,1,...] 
               + nuts.variables['flagellates'][0,1,...]) * mask, 0)).mean()
    model_average = np.append(model_average, chl)
In [15]:
t0=datetime.datetime(2017,1,1)
log_model_average = np.array([])
for date in [t0+datetime.timedelta(days=ii*1) for ii in range(365)]:
    sub_dir = date.strftime('%d%b%y').lower()
    datestr = date.strftime('%Y%m%d')
    fname = 'SalishSea_1d_{}_{}_ptrc_T.nc'.format(datestr, datestr)
    nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
    chl = np.log10((1.6 * np.ma.masked_equal((nuts.variables['diatoms'][0, 1, ...] 
               + nuts.variables['ciliates'][0,1,...] 
               + nuts.variables['flagellates'][0,1,...]) * mask, 0))).mean()
    log_model_average = np.append(log_model_average, chl)
In [16]:
t0=datetime.datetime(2015,1,1)
model_average_2015 = np.array([])
for date in [t0+datetime.timedelta(days=ii*1) for ii in range(365)]:
    sub_dir = date.strftime('%d%b%y').lower()
    datestr = date.strftime('%Y%m%d')
    fname = 'SalishSea_1d_{}_{}_ptrc_T.nc'.format(datestr, datestr)
    nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
    chl = (1.6 * np.ma.masked_equal((nuts.variables['diatoms'][0, 1, ...] 
               + nuts.variables['ciliates'][0,1,...] 
               + nuts.variables['flagellates'][0,1,...]) * mask, 0)).mean()
    model_average_2015 = np.append(model_average_2015, chl)
In [17]:
t0=datetime.datetime(2015,1,1)
log_model_average_2015 = np.array([])
for date in [t0+datetime.timedelta(days=ii*1) for ii in range(365)]:
    sub_dir = date.strftime('%d%b%y').lower()
    datestr = date.strftime('%Y%m%d')
    fname = 'SalishSea_1d_{}_{}_ptrc_T.nc'.format(datestr, datestr)
    nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
    chl = np.log10((1.6 * np.ma.masked_equal((nuts.variables['diatoms'][0, 1, ...] 
               + nuts.variables['ciliates'][0,1,...] 
               + nuts.variables['flagellates'][0,1,...]) * mask, 0))).mean()
    log_model_average_2015 = np.append(log_model_average_2015, chl)
In [18]:
t0=datetime.datetime(2016,1,1)
model_average_2016 = np.array([])
for date in [t0+datetime.timedelta(days=ii*1) for ii in range(366)]:
    sub_dir = date.strftime('%d%b%y').lower()
    datestr = date.strftime('%Y%m%d')
    fname = 'SalishSea_1d_{}_{}_ptrc_T.nc'.format(datestr, datestr)
    nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
    chl = (1.6 * np.ma.masked_equal((nuts.variables['diatoms'][0, 1, ...] 
               + nuts.variables['ciliates'][0,1,...] 
               + nuts.variables['flagellates'][0,1,...]) * mask, 0)).mean()
    model_average_2016 = np.append(model_average_2016, chl)
In [19]:
t0=datetime.datetime(2016,1,1)
log_model_average_2016 = np.array([])
for date in [t0+datetime.timedelta(days=ii*1) for ii in range(366)]:
    sub_dir = date.strftime('%d%b%y').lower()
    datestr = date.strftime('%Y%m%d')
    fname = 'SalishSea_1d_{}_{}_ptrc_T.nc'.format(datestr, datestr)
    nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
    chl = np.log10((1.6 * np.ma.masked_equal((nuts.variables['diatoms'][0, 1, ...] 
               + nuts.variables['ciliates'][0,1,...] 
               + nuts.variables['flagellates'][0,1,...]) * mask, 0))).mean()
    log_model_average_2016 = np.append(log_model_average_2016, chl)
In [20]:
list_of_days = np.array([list_of_date[n].date() for n in range(list_of_date.shape[0])])
list_of_days.shape
Out[20]:
(19836,)
In [21]:
log_chl = np.log10(list_of_ferry_chl)
daily_ferry_logchl = np.array([])
for day in np.unique(list_of_days):
    daily_ferry_logchl = np.append(daily_ferry_logchl, 
                                   log_chl[list_of_days == day].mean())
In [22]:
daily_ferry_chl = np.array([])
for day in np.unique(list_of_days):
    daily_ferry_chl = np.append(daily_ferry_chl, 
                                list_of_ferry_chl[list_of_days == day].mean())
In [23]:
ferry_2015 = pickle.load(open('daily_ferry_chl_2015.pkl', 'rb'))
days_2015 = pickle.load(open('days_2015.pkl', 'rb'))
ferry_2016 = pickle.load(open('daily_ferry_chl_2016.pkl', 'rb'))
days_2016 = pickle.load(open('days_2016.pkl', 'rb'))
In [24]:
log_ferry_2015 = pickle.load(open('daily_ferry_logchl_2015.pkl', 'rb'))
log_ferry_2016 = pickle.load(open('daily_ferry_logchl_2016.pkl', 'rb'))
In [25]:
import matplotlib as mpl
mpl.rcParams['font.size'] = 18
mpl.rcParams['axes.titlesize'] = 18
In [26]:
fig, ax = plt.subplots(figsize = (12,12))
viz_tools.plot_coastline(ax, grid, coords = 'map')
ax.plot((-123.78116, -123.78116), (48, 50), '--', color = 'grey')
ax.plot((-123.694065, -123.694065), (48, 50), '--', color = 'grey')
ax.set_xlim(-124.25, -122.75)
ax.set_ylim(48, 49.5);

Arithmetic Mean - Full

In [27]:
fig, ax = plt.subplots(5,1, figsize = (24,24))
ax[0].plot(np.unique(list_of_days), daily_ferry_chl, 's' , color = 'SeaGreen', label = '2017 ferry')
ax[0].plot([datetime.date(2017,1,1)+datetime.timedelta(days=ii*1) for ii in range(365)], 
                model_average, color = 'SeaGreen', label = '2017 model')
ax[0].set_ylim(0, 30)
ax[0].set_xlim(datetime.date(2017,1,1), datetime.date(2017,12,31))
ax[0].set_title('2017')
ax[0].legend();

ax[1].plot(days_2016, ferry_2016, 's' , color = 'SeaGreen', label = '2016 ferry')
ax[1].plot([datetime.date(2016,1,1)+datetime.timedelta(days=ii*1) for ii in range(366)], 
                model_average_2016, color = 'SeaGreen', label = '2016 model')
ax[1].set_ylim(0, 30)
ax[1].set_xlim(datetime.date(2016,1,1), datetime.date(2016,12,31))
ax[1].set_title('2016')
ax[1].legend();

ax[2].plot(days_2015, ferry_2015, 's', color = 'SeaGreen',
          label = '2015 ferry')
ax[2].plot([datetime.date(2015,1,1)+datetime.timedelta(days=ii*1) for ii in range(365)], 
                model_average_2015,color = 'SeaGreen', label = '2015 model')
ax[2].set_ylim(0, 30)
ax[2].set_xlim(datetime.date(2015,1,1), datetime.date(2015,12,31))
ax[2].set_title('2015')
ax[2].legend();

ax[3].plot([(np.unique(list_of_days) - datetime.date(2017,1,1))[n].days 
            for n in range(np.unique(list_of_days).shape[0])], 
           daily_ferry_chl, 'o' , color = 'DeepSkyBlue', label = '2017 ferry')
ax[3].plot([(days_2016 - datetime.date(2016,1,1))[n].days 
            for n in range(days_2016.shape[0])], 
           ferry_2016, 'o' , color = 'Coral', label = '2016 ferry')
ax[3].plot([(days_2015 - datetime.date(2015,1,1))[n].days 
            for n in range(days_2015.shape[0])], 
           ferry_2015, 'o' , color = 'DarkViolet', label = '2015 ferry')
ax[3].set_ylim(0, 30)
ax[3].set_xlim(0,366)
ax[3].set_title('Ferry')
ax[3].legend();

ax[4].plot(range(365), model_average, color = 'DeepSkyBlue', label = '2017 model')
ax[4].plot(range(366), model_average_2016, color = 'Coral', label = '2016 model')
ax[4].plot(range(365), model_average_2015, color = 'DarkViolet', label = '2015 model')

ax[4].set_ylim(0, 30)
ax[4].set_xlim(0,366)
ax[4].set_title('Model Chl')
ax[4].legend();

Arithmetic Mean - Fall

In [28]:
fig, ax = plt.subplots(5,1, figsize = (24,24))
ax[0].plot(np.unique(list_of_days), daily_ferry_chl, 's' , color = 'SeaGreen', label = '2017 ferry')
ax[0].plot([datetime.date(2017,1,1)+datetime.timedelta(days=ii*1) for ii in range(365)], 
                model_average, color = 'SeaGreen', label = '2017 model')
ax[0].set_ylim(0, 30)
ax[0].set_xlim(datetime.date(2017,8,1), datetime.date(2017,12,31))
ax[0].set_title('2017')
ax[0].legend();

ax[1].plot(days_2016, ferry_2016, 's' , color = 'SeaGreen', label = '2016 ferry')
ax[1].plot([datetime.date(2016,1,1)+datetime.timedelta(days=ii*1) for ii in range(366)], 
                model_average_2016, color = 'SeaGreen', label = '2016 model')
ax[1].set_ylim(0, 30)
ax[1].set_xlim(datetime.date(2016,8,1), datetime.date(2016,12,31))
ax[1].set_title('2016')
ax[1].legend();

ax[2].plot(days_2015, ferry_2015, 's', color = 'SeaGreen',
          label = '2015 ferry')
ax[2].plot([datetime.date(2015,1,1)+datetime.timedelta(days=ii*1) for ii in range(365)], 
                model_average_2015,color = 'SeaGreen', label = '2015 model')
ax[2].set_ylim(0, 30)
ax[2].set_xlim(datetime.date(2015,8,1), datetime.date(2015,12,31))
ax[2].set_title('2015')
ax[2].legend();

ax[3].plot([(np.unique(list_of_days) - datetime.date(2017,1,1))[n].days 
            for n in range(np.unique(list_of_days).shape[0])], 
           daily_ferry_chl, 'o' , color = 'DeepSkyBlue', label = '2017 ferry')
ax[3].plot([(days_2016 - datetime.date(2016,1,1))[n].days 
            for n in range(days_2016.shape[0])], 
           ferry_2016, 'o' , color = 'Coral', label = '2016 ferry')
ax[3].plot([(days_2015 - datetime.date(2015,1,1))[n].days 
            for n in range(days_2015.shape[0])], 
           ferry_2015, 'o' , color = 'DarkViolet', label = '2015 ferry')
ax[3].set_ylim(0, 30)
ax[3].set_xlim(212,366)
ax[3].set_title('Ferry')
ax[3].legend();

ax[4].plot(range(365), model_average, color = 'DeepSkyBlue', label = '2017 model')
ax[4].plot(range(366), model_average_2016, color = 'Coral', label = '2016 model')
ax[4].plot(range(365), model_average_2015, color = 'DarkViolet', label = '2015 model')

ax[4].set_ylim(0, 30)
ax[4].set_xlim(212,366)
ax[4].set_title('Model Chl')
ax[4].legend();

Geometric Mean - Full

In [29]:
fig, ax = plt.subplots(5,1, figsize = (24,24))
ax[0].plot(np.unique(list_of_days), np.power(10,daily_ferry_logchl), 
           's' , color = 'SeaGreen', label = '2017 ferry log cgl')
ax[0].plot([datetime.date(2017,1,1)+datetime.timedelta(days=ii*1) for ii in range(365)], 
                np.power(10,log_model_average), color = 'SeaGreen', label = '2017 model')
ax[0].set_ylim(0, 30)
ax[0].set_xlim(datetime.date(2017,1,1), datetime.date(2017,12,31))
ax[0].set_title('2017')
ax[0].legend();

ax[1].plot(days_2016, np.power(10,log_ferry_2016), 's' , color = 'SeaGreen', label = '2016 ferry log chl')
ax[1].plot([datetime.date(2016,1,1)+datetime.timedelta(days=ii*1) for ii in range(366)], 
                np.power(10,log_model_average_2016), color = 'SeaGreen', label = '2016 model')
ax[1].set_ylim(0, 30)
ax[1].set_xlim(datetime.date(2016,1,1), datetime.date(2016,12,31))
ax[1].set_title('2016')
ax[1].legend();

ax[2].plot(days_2015, np.power(10,log_ferry_2015), 's', color = 'SeaGreen',
          label = '2015 ferry')
ax[2].plot([datetime.date(2015,1,1)+datetime.timedelta(days=ii*1) for ii in range(365)], 
                np.power(10,log_model_average_2015),color = 'SeaGreen', label = '2015 model')
ax[2].set_ylim(0, 30)
ax[2].set_xlim(datetime.date(2015,1,1), datetime.date(2015,12,31))
ax[2].set_title('2015')
ax[2].legend();

ax[3].plot([(np.unique(list_of_days) - datetime.date(2017,1,1))[n].days 
            for n in range(np.unique(list_of_days).shape[0])], 
           np.power(10,daily_ferry_logchl), 'o' , color = 'DeepSkyBlue', label = '2017 ferry log chl')
ax[3].plot([(days_2016 - datetime.date(2016,1,1))[n].days 
            for n in range(days_2016.shape[0])], 
           np.power(10,log_ferry_2016), 'o' , color = 'Coral', label = '2016 ferry log chl')
ax[3].plot([(days_2015 - datetime.date(2015,1,1))[n].days 
            for n in range(days_2015.shape[0])], 
           np.power(10,log_ferry_2015), 'o' , color = 'DarkViolet', label = '2015 ferry log chl')
ax[3].set_ylim(0, 30)
ax[3].set_xlim(0,366)
ax[3].set_title('Ferry')
ax[3].legend();

ax[4].plot(range(365), np.power(10,log_model_average), color = 'DeepSkyBlue', label = '2017 model')
ax[4].plot(range(366), np.power(10,log_model_average_2016), color = 'Coral', label = '2016 model')
ax[4].plot(range(365), np.power(10,log_model_average_2015), color = 'DarkViolet', label = '2015 model')

ax[4].set_ylim(0, 30)
ax[4].set_xlim(0,366)
ax[4].set_title('Model Chl')
ax[4].legend();