In [153]:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import netCDF4 as nc
import seaborn as sns
import matplotlib.colors as mcolors
import glob
import os
import xarray as xr
import datetime
from salishsea_tools import viz_tools, tidetools, geo_tools, gsw_calls
import ONC_patrols as onc
%matplotlib inline
In [2]:
nowcast = xr.open_dataset('https://salishsea.eos.ubc.ca/erddap/griddap/ubcSSg3DTracerFields1hV17-02')
In [3]:
nowcast_salinity = nowcast.salinity
nowcast_temperature = nowcast.temperature
In [4]:
grid_B_new=nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
mesh_mask_new=nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/mesh_mask201702.nc')
badQC=[0,3,4,9]
In [5]:
bathy, X, Y = tidetools.get_bathy_data(grid_B_new)
In [ ]:
In [6]:
Victoria_file = "/ocean/nsoontie/MEOPAR/ONC/Patrols/Victoria_Patrol9_CTD_20150220T182104Z_20151001T195734Z-Corrected.csv"
In [7]:
V_data = onc.load_patrol_csv(Victoria_file)
V_data = onc.exclude_bad(V_data,['Practical Salinity Corrected QC Flag '], badQC)
V_data = onc.divide_into_casts(V_data)
In [8]:
V_data.keys()
Out[8]:
In [9]:
lon = V_data['Longitude Corrected (deg)'].values
lat = V_data['Latitude Corrected (deg)'].values
time = V_data['time'].values
salinity = V_data['Practical Salinity Corrected (psu)'].values
depth = V_data['Depth Corrected (m)'].values
temperature = V_data['Temperature Corrected (C)'].values
In [10]:
BS_file = "/ocean/nsoontie/MEOPAR/ONC/Patrols/BaynesSound_Patrol4_CTD_20150218T180311Z_20160301T214950Z-Corrected.csv"
BS_data = onc.load_patrol_csv(BS_file)
BS_data = onc.exclude_bad(BS_data,['Practical Salinity Corrected QC Flag '], badQC)
BS_data = onc.divide_into_casts(BS_data)
GI_file = '/ocean/nsoontie/MEOPAR/ONC/Patrols/GalianoIsland_Patrol9_CTD_20160303T211549Z_20160304T010923Z-Corrected.csv'
GI_data = onc.load_patrol_csv(GI_file)
GI_data = onc.exclude_bad(GI_data,['Practical Salinity Corrected QC Flag '], badQC)
GI_data = onc.divide_into_casts(GI_data)
NQ_file = '/ocean/nsoontie/MEOPAR/ONC/Patrols/Nanaimo_Qualicum_Patrol5_CTD_20150220T163040Z_20151005T204619Z-Corrected.csv'
NQ_data = onc.load_patrol_csv(NQ_file)
NQ_data = onc.exclude_bad(NQ_data,['Practical Salinity Corrected QC Flag '], badQC)
NQ_data = onc.divide_into_casts(NQ_data)
S_file = '/ocean/nsoontie/MEOPAR/ONC/Patrols/Steveston_Patrol8_CTD_20150525T171143Z_20151005T222329Z-Corrected.csv'
S_data = onc.load_patrol_csv(S_file)
S_data = onc.exclude_bad(S_data,['Practical Salinity Corrected QC Flag '], badQC)
S_data = onc.divide_into_casts(S_data)
In [11]:
for data in [BS_data, GI_data, NQ_data, S_data]:
lon = np.append(lon, data['Longitude Corrected (deg)'].values)
lat = np.append(lat, data['Latitude Corrected (deg)'].values)
time = np.append(time, data['time'].values)
salinity = np.append(salinity, data['Practical Salinity Corrected (psu)'].values)
depth = np.append(depth, data['Depth Corrected (m)'].values)
temperature = np.append(temperature, data['Temperature Corrected (C)'].values)
In [12]:
print(lon.shape)
print(depth.shape)
In [74]:
deptht = nowcast.depth.values
tmask = mesh_mask_new.variables['tmask'][:]
In [76]:
list_of_sals2 = np.array([])
list_of_temps2 = np.array([])
land_index = np.array([])
for n in range(39443):
Yind, Xind = geo_tools.find_closest_model_point(lon[n], lat[n], X, Y, land_mask = bathy.mask)
date = pd.Timestamp(time[n]).to_pydatetime()
land = np.abs((deptht - depth[n])).argmin()
if tmask[0,land,Yind,Xind] == 1:
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_val = (delta * ((nowcast_salinity
.sel(time = before, depth = depth[n], method = 'nearest')
.isel(gridY = Yind, gridX = Xind)).values) +
(1- delta)*((nowcast.salinity
.sel(time = before + datetime.timedelta(hours=1),
depth = depth[n], method = 'nearest')
.isel(gridY = Yind, gridX = Xind)).values))
t_val = (delta * ((nowcast_temperature
.sel(time = before, depth = depth[n], method = 'nearest')
.isel(gridY = Yind, gridX = Xind)).values) +
(1- delta)*((nowcast_temperature
.sel(time = before + datetime.timedelta(hours=1),
depth = depth[n], method = 'nearest')
.isel(gridY = Yind, gridX = Xind)).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_val = (delta * ((nowcast_salinity
.sel(time = before, depth = depth[n], method = 'nearest')
.isel(gridY = Yind, gridX = Xind)).values) +
(1- delta)*((nowcast.salinity
.sel(time = before + datetime.timedelta(hours=1),
depth = depth[n], method = 'nearest')
.isel(gridY = Yind, gridX = Xind)).values))
t_val = (delta * ((nowcast_temperature
.sel(time = before, depth = depth[n], method = 'nearest')
.isel(gridY = Yind, gridX = Xind)).values) +
(1- delta)*((nowcast_temperature
.sel(time = before + datetime.timedelta(hours=1),
depth = depth[n], method = 'nearest')
.isel(gridY = Yind, gridX = Xind)).values))
list_of_temps2 = np.append(list_of_temps2, t_val)
list_of_sals2 = np.append(list_of_sals2, s_val)
else:
land_index = np.append(land_index, n)
In [78]:
salinity2 = np.delete(salinity, land_index)
temperature2 = np.delete(temperature, land_index)
In [79]:
print(salinity2.shape)
print(list_of_sals2.shape)
In [155]:
sr_obs = gsw_calls.generic_gsw_caller('gsw_SR_from_SP.m', [salinity2])
In [156]:
fig, ax = plt.subplots(figsize = (8,8))
ax.plot(sr_obs, list_of_sals2, 'b.')
ax.plot(np.arange(0,35), np.arange(0,35), 'r-')
ax.set_ylabel('Nowcast-green (g/kg)')
ax.set_xlabel('Citizen Science Data (g/kg)')
ax.set_title('Modeled Vs Observed Salinity')
ax.grid('on')
In [157]:
print('bias = ' + str(np.mean(sr_obs) - np.mean(list_of_sals2)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_sals2 - sr_obs)**2) / len(list_of_sals2))))
xbar = np.mean(sr_obs)
print('Willmott = ' + str(1-(np.sum((list_of_sals2 - sr_obs)**2) /
np.sum((np.abs(list_of_sals2 - xbar) + np.abs(sr_obs - xbar))**2))))
In [159]:
ct_obs = gsw_calls.generic_gsw_caller('gsw_CT_from_pt.m', [sr_obs, temperature2])
In [162]:
fig, ax = plt.subplots(figsize = (8,8))
ax.plot(ct_obs, list_of_temps2, 'b.')
ax.plot(np.arange(7,22), np.arange(7,22), 'b-')
ax.set_ylabel('Nowcast-green (deg C)')
ax.set_xlabel('Citizen Science Data (deg C)')
ax.set_title('Model Vs Observed Temperature')
ax.grid('on')
In [161]:
print('bias = ' + str(np.mean(ct_obs) - np.mean(list_of_temps2)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_temps2 - ct_obs)**2) / len(list_of_temps2))))
xbar = np.mean(ct_obs)
print('Willmott = ' + str(1-(np.sum((list_of_temps2 - ct_obs)**2) /
np.sum((np.abs(list_of_temps2 - xbar) + np.abs(ct_obs - xbar))**2))))
In [ ]:
In [ ]:
In [ ]: