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
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>''')
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')
HINDCAST_PATH= '/results/SalishSea/nowcast-green/'
f = pd.read_excel('/ocean/eolson/MEOPAR/obs/PSFCitSci/Chla_2015PSFSalish_Sea_22.01.2018vers_8_CN_edits.xls',header=6)
f = f.drop(f.index[[0]])
f.shape
(137, 9)
f = f.dropna(subset = ['Date sampled (mm/dd/yyyy)', 'Time of Day (Local)',
'Latitude', 'Longitude', 'Depth', 'Chl a'])
f.shape
(51, 9)
local = pytz.timezone ("America/Los_Angeles")
f.loc[f['Date sampled (mm/dd/yyyy)'] == '07/15//2015', 'Date sampled (mm/dd/yyyy)'] = '07/15/2015'
f.loc[f['Date sampled (mm/dd/yyyy)'] == '05/23//2015', 'Date sampled (mm/dd/yyyy)'] = '05/23/2015'
f.loc[f['Date sampled (mm/dd/yyyy)'] == '05/17//2015', 'Date sampled (mm/dd/yyyy)'] = '05/17/2015'
f['Latitude'].unique()[0][4:]
'44.200'
f.head()
Date sampled (mm/dd/yyyy) | Station Name | Time of Day (Local) | Latitude | Longitude | Depth | Chl a | Phaeophytin | Comments | Date sampled (mm/dd/yyy) | |
---|---|---|---|---|---|---|---|---|---|---|
1 | 09/16/2015 | CBE-2 | 11:08:00 | 48° 44.200 | 123° 34.300 | 5 | 3.08 | 0.32 | Averaged duplicates | NaN |
2 | 09/14/2015 | LS-2 | 12:19:00 | 48° 58.400 | 123° 46.000 | 5 | 3.93919 | 0.957782 | Averaged duplicates | NaN |
3 | 08/09/2015 | CR-8 | 15:31:00 | 49° 57.700 | 125° 08.800 | 5 | 0.420117 | 0.238546 | Averaged duplicates | NaN |
5 | 06/16/2015 | CR-8 | 10:15:00 | 49° 57.700 | 125° 08.800 | 5 | 0.947415 | 0.527436 | Averaged duplicates | NaN |
6 | 06/23/2015 | CR-8 | 10:15:00 | 49° 57.700 | 125° 08.800 | 5 | 0.493184 | 0.438475 | Averaged duplicates | NaN |
list_of_lons = np.array([])
list_of_lats = np.array([])
list_of_datetimes = np.array([])
list_of_cs_chl = np.array([])
list_of_model_chl = np.array([])
list_of_depths = np.array([])
for n in f.index:
decLat = float(f.Latitude[n][:2]) + float(f['Latitude'][n][4:])/60
decLon = (float(f.Longitude[n][:3]) + float(f['Longitude'][n][5:])/60) * -1
Yind, Xind = geo_tools.find_closest_model_point(decLon, decLat,
X, Y, land_mask = bathy.mask)
if f['Depth'][n] == 5:
depth = 4
elif f['Depth'][n] == 20:
depth = 18
elif f['Depth'][n] == 10:
depth = 9
if mesh.variables['tmask'][0,depth,Yind, Xind] == 1:
try:
local_datetime = datetime.datetime.combine(datetime.datetime.strptime(
f['Date sampled (mm/dd/yyyy)'][n], '%m/%d/%Y'),f['Time of Day (Local)'][n])
except (TypeError):
local_datetime = (datetime.datetime.combine(f['Date sampled (mm/dd/yyyy)'][n],
f['Time of Day (Local)'][n]))
date = local.localize(local_datetime, is_dst=True).astimezone(pytz.utc)
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)
sub_dir = before.strftime('%d%b%y').lower()
datestr = before.strftime('%Y%m%d')
fname = 'SalishSea_1h_{}_{}_ptrc_T.nc'.format(datestr, datestr)
nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
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
chl_val = 1.6*((1-delta)*(nuts.variables['diatoms'][before.hour, depth, Yind, Xind]
+ nuts.variables['ciliates'][before.hour,depth,Yind, Xind]
+ nuts.variables['flagellates'][before.hour,depth,Yind,Xind]) +
(delta)*(nuts2.variables['diatoms'][after.hour, depth, Yind, Xind]
+ nuts2.variables['ciliates'][after.hour,depth,Yind, Xind]
+ nuts2.variables['flagellates'][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)
sub_dir = before.strftime('%d%b%y').lower()
datestr = before.strftime('%Y%m%d')
fname = 'SalishSea_1h_{}_{}_ptrc_T.nc'.format(datestr, datestr)
nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
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
chl_val = 1.6*(delta*(nuts.variables['diatoms'][before.hour, depth, Yind, Xind]
+ nuts.variables['ciliates'][before.hour,depth,Yind, Xind]
+ nuts.variables['flagellates'][before.hour,depth,Yind,Xind]) +
(1- delta)*(nuts2.variables['diatoms'][after.hour, depth, Yind, Xind]
+ nuts2.variables['ciliates'][after.hour,depth,Yind, Xind]
+ nuts2.variables['flagellates'][after.hour,depth,Yind,Xind]))
list_of_lons = np.append(list_of_lons, decLon)
list_of_lats = np.append(list_of_lats, decLat)
list_of_datetimes = np.append(list_of_datetimes, date)
list_of_cs_chl = np.append(list_of_cs_chl, f['Chl a'][n])
list_of_model_chl = np.append(list_of_model_chl, chl_val)
list_of_depths = np.append(list_of_depths, depth)
import matplotlib as mpl
mpl.rcParams['font.size'] = 12
mpl.rcParams['axes.titlesize'] = 12
fig, ax = plt.subplots(figsize = (8,8))
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.2, 50.5)
ax.set_xlim(-125.5, -122.5);
list_of_cs_chl.shape
(51,)
fig, ax = plt.subplots(figsize = (8,8))
ax.plot(list_of_cs_chl, list_of_model_chl, 'b.', alpha = 0.5)
ax.plot(np.arange(0,10), color = 'grey')
ax.grid('on')
ax.set_title('Citizen Science Chl 2015, depth = 5m')
ax.set_xlabel('Citizen Science')
ax.set_ylabel('Nowcast-green');
print('bias = ' + str(-np.mean(list_of_cs_chl) + np.mean(list_of_model_chl)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_chl - list_of_cs_chl)**2) /
len(list_of_cs_chl))))
xbar = np.mean(list_of_cs_chl)
print('Willmott = ' + str(1-(np.sum((list_of_model_chl - list_of_cs_chl)**2) /
np.sum((np.abs(list_of_model_chl - xbar)
+ np.abs(list_of_cs_chl - xbar))**2))))
bias = 4.09138969437 RMSE = 5.16042682101 Willmott = 0.279270060415
fig, ax = plt.subplots(figsize = (20,8))
ax.plot(list_of_lats, list_of_model_chl - list_of_cs_chl, 'ro', alpha =0.5)
ax.grid('on')
ax.set_xlabel('lat', fontsize = 15)
ax.set_ylabel('Model - Observed',fontsize = 15)
ax.set_title('Chl', fontsize = 20)
<matplotlib.text.Text at 0x7f6a80836ac8>
fig, ax = plt.subplots(figsize = (20,8))
ax.plot(list_of_lons, list_of_model_chl - list_of_cs_chl, 'ro', alpha =0.5)
ax.grid('on')
ax.set_xlabel('lon', fontsize = 15)
ax.set_ylabel('Model - Observed',fontsize = 15)
ax.set_title('Chl', fontsize = 20)
<matplotlib.text.Text at 0x7f6a80756a20>
fig, ax = plt.subplots(figsize = (20,8))
ax.plot(list_of_datetimes, list_of_model_chl - list_of_cs_chl, 'ro', alpha =0.5)
ax.grid('on')
ax.set_xlabel('Date', fontsize = 15)
ax.set_ylabel('Model - Observed',fontsize = 15)
ax.set_title('Chl', fontsize = 20)
<matplotlib.text.Text at 0x7f6a806fadd8>