In [103]:
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
from salishsea_tools import geo_tools, tidetools, viz_tools, loadDataFRP
import matplotlib.cm as cm
%matplotlib inline
In [2]:
base=nc.Dataset('/ocean/vdo/MEOPAR/completed-runs/threemonthbase/testD/SalishSea_1h_20170406_20170415_grid_T_20170410-20170410.nc')
testa=nc.Dataset('/ocean/vdo/MEOPAR/completed-runs/threemonthsa/test44/SalishSea_1h_20170406_20170415_grid_T_20170410-20170410.nc')
testb=nc.Dataset('/ocean/vdo/MEOPAR/completed-runs/threemonthsb/testd/SalishSea_1h_20170406_20170415_grid_T_20170410-20170410.nc')
In [62]:
Bathymetry = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
bathy, X, Y = tidetools.get_bathy_data(Bathymetry)
In [65]:
mesh = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/mesh_mask201702.nc')
tmask = mesh.variables['tmask'][:]
In [4]:
stationdata, casts = loadDataFRP.loadDataFRP_SSGrid()
In [27]:
stationdata[['Station', 'Date_UTC', 'Time_UTC_hhmmss']]
Out[27]:
In [60]:
max_dsdz = np.array([])
depth_of_max_dsdz = np.array([])
for n in range(1,19):
if n == 14:
depths = casts[14.1].dCast['depth_m']
dz = casts[14.1].dCast['depth_m'].values[1:] - casts[14.1].dCast['depth_m'].values[:-1]
ds = casts[14.1].dCast['gsw_srA0'].values[1:] - casts[14.1].dCast['gsw_srA0'].values[:-1]
depth_of_max_dsdz = np.append(depth_of_max_dsdz,
(depths[np.nanargmax(ds / dz)]
+ depths[np.nanargmax(ds / dz) + 1])/2)
max_dsdz = np.append(max_dsdz, np.nanmax(ds / dz))
depths = casts[14.2].uCast['depth_m']
dz = casts[14.2].uCast['depth_m'].values[1:] - casts[14.2].uCast['depth_m'].values[:-1]
ds = casts[14.2].uCast['gsw_srA0'].values[1:] - casts[14.2].uCast['gsw_srA0'].values[:-1]
depth_of_max_dsdz = np.append(depth_of_max_dsdz,
(depths[np.nanargmax(ds / dz)]
+ depths[np.nanargmax(ds / dz) + 1])/2)
max_dsdz = np.append(max_dsdz, np.nanmax(ds / dz))
else:
depths = casts[n].dCast['depth_m']
dz = casts[n].dCast['depth_m'].values[1:] - casts[n].dCast['depth_m'].values[:-1]
ds = casts[n].dCast['gsw_srA0'].values[1:] - casts[n].dCast['gsw_srA0'].values[:-1]
depth_of_max_dsdz = np.append(depth_of_max_dsdz,
(depths[np.nanargmax(ds / dz)]
+ depths[np.nanargmax(ds / dz) + 1])/2)
max_dsdz = np.append(max_dsdz, np.nanmax(ds / dz))
depths = casts[n].uCast['depth_m']
dz = casts[n].uCast['depth_m'].values[1:] - casts[n].uCast['depth_m'].values[:-1]
ds = casts[n].uCast['gsw_srA0'].values[1:] - casts[n].uCast['gsw_srA0'].values[:-1]
depth_of_max_dsdz = np.append(depth_of_max_dsdz,
(depths[np.nanargmax(ds / dz)]
+ depths[np.nanargmax(ds / dz) + 1])/2)
max_dsdz = np.append(max_dsdz, np.nanmax(ds / dz))
print(max_dsdz)
print(depth_of_max_dsdz)
In [151]:
depths_base = np.array([])
max_dsdz_base = np.array([])
depths_a = np.array([])
max_dsdz_a = np.array([])
depths_b = np.array([])
max_dsdz_b = np.array([])
for n in range(19):
station = stationdata.iloc[[n]]
cast = casts[station['Station'].values[0]]
Yind, Xind = geo_tools.find_closest_model_point(station['LonDecDeg'].values[0],
station['LatDecDeg'].values[0],
X, Y, land_mask = bathy.mask)
if n == 13:
shape_depth = cast.dCast['depth_m'].values.shape[0]
else:
shape_depth = cast.uCast['depth_m'].values.shape[0]
if n < 9:
base=nc.Dataset(
'/ocean/vdo/MEOPAR/completed-runs/threemonthbase/testI/SalishSea_1h_20170526_20170604_grid_T_20170531-20170531.nc')
testa=nc.Dataset(
'/ocean/vdo/MEOPAR/completed-runs/threemonthsa/test99/SalishSea_1h_20170526_20170604_grid_T_20170531-20170531.nc')
testb=nc.Dataset(
'/ocean/vdo/MEOPAR/completed-runs/threemonthsb/testi/SalishSea_1h_20170526_20170604_grid_T_20170531-20170531.nc')
else:
base=nc.Dataset(
'/ocean/vdo/MEOPAR/completed-runs/threemonthbase/testD/SalishSea_1h_20170406_20170415_grid_T_20170410-20170410.nc')
testa=nc.Dataset(
'/ocean/vdo/MEOPAR/completed-runs/threemonthsa/test44/SalishSea_1h_20170406_20170415_grid_T_20170410-20170410.nc')
testb=nc.Dataset(
'/ocean/vdo/MEOPAR/completed-runs/threemonthsb/testd/SalishSea_1h_20170406_20170415_grid_T_20170410-20170410.nc')
if int(station['Time_UTC_hhmmss'].values[0][3:5]) < 30:
delta = (30 + int(station['Time_UTC_hhmmss'].values[0][3:5])) / 60
before = int(station['Time_UTC_hhmmss'].values[0][:2]) - 1
else:
delta = (int(station['Time_UTC_hhmmss'].values[0][3:5]) - 30) / 60
before = int(station['Time_UTC_hhmmss'].values[0][:2])
pt_mask = tmask[0,:shape_depth,Yind,Xind]
deptht = base.variables['deptht'][:shape_depth]
base_sal = np.ma.masked_array(delta * base.variables['vosaline'][before,:shape_depth,Yind, Xind]
+ (1-delta)*base.variables['vosaline'][before+1,:shape_depth,Yind, Xind],
mask = 1-pt_mask)
a_sal = np.ma.masked_array(delta*testa.variables['vosaline'][before,:shape_depth,Yind, Xind]
+(1-delta)*testa.variables['vosaline'][before+1,:shape_depth,Yind,Xind],
mask = 1-pt_mask)
b_sal = np.ma.masked_array(delta*testb.variables['vosaline'][before,:shape_depth,Yind, Xind]
+(1-delta)*testb.variables['vosaline'][before+1,:shape_depth,Yind,Xind],
mask = 1-pt_mask)
dz = deptht[1:] - deptht[:-1]
ds_base = base_sal[1:] - base_sal[:-1]
ds_a = a_sal[1:] - a_sal[:-1]
ds_b = b_sal[1:] - b_sal[:-1]
if n == 13:
depths_base = np.append(depths_base, (deptht[np.nanargmax(ds_base / dz)]
+ deptht[np.nanargmax(ds_base / dz) + 1])/2)
max_dsdz_base = np.append(max_dsdz_base, np.nanmax(ds_base / dz))
depths_a = np.append(depths_a, (deptht[np.nanargmax(ds_a / dz)]
+ deptht[np.nanargmax(ds_a / dz) + 1])/2)
max_dsdz_a = np.append(max_dsdz_a, np.nanmax(ds_a / dz))
depths_b = np.append(depths_b, (deptht[np.nanargmax(ds_b / dz)]
+ deptht[np.nanargmax(ds_b / dz) + 1])/2)
max_dsdz_b = np.append(max_dsdz_b, np.nanmax(ds_b / dz))
elif n == 14:
depths_base = np.append(depths_base, (deptht[np.nanargmax(ds_base / dz)]
+ deptht[np.nanargmax(ds_base / dz) + 1])/2)
max_dsdz_base = np.append(max_dsdz_base, np.nanmax(ds_base / dz))
depths_a = np.append(depths_a, (deptht[np.nanargmax(ds_a / dz)]
+ deptht[np.nanargmax(ds_a / dz) + 1])/2)
max_dsdz_a = np.append(max_dsdz_a, np.nanmax(ds_a / dz))
depths_b = np.append(depths_b, (deptht[np.nanargmax(ds_b / dz)]
+ deptht[np.nanargmax(ds_b / dz) + 1])/2)
max_dsdz_b = np.append(max_dsdz_b, np.nanmax(ds_b / dz))
else:
depths_base = np.append(depths_base, (deptht[np.nanargmax(ds_base / dz)]
+ deptht[np.nanargmax(ds_base / dz) + 1])/2)
max_dsdz_base = np.append(max_dsdz_base, np.nanmax(ds_base / dz))
depths_base = np.append(depths_base, (deptht[np.nanargmax(ds_base / dz)]
+ deptht[np.nanargmax(ds_base / dz) + 1])/2)
max_dsdz_base = np.append(max_dsdz_base, np.nanmax(ds_base / dz))
depths_a = np.append(depths_a, (deptht[np.nanargmax(ds_a / dz)]
+ deptht[np.nanargmax(ds_a / dz) + 1])/2)
max_dsdz_a = np.append(max_dsdz_a, np.nanmax(ds_a / dz))
depths_a = np.append(depths_a, (deptht[np.nanargmax(ds_a / dz)]
+ deptht[np.nanargmax(ds_a / dz) + 1])/2)
max_dsdz_a = np.append(max_dsdz_a, np.nanmax(ds_a / dz))
depths_b = np.append(depths_b, (deptht[np.nanargmax(ds_b / dz)]
+ deptht[np.nanargmax(ds_b / dz) + 1])/2)
max_dsdz_b = np.append(max_dsdz_b, np.nanmax(ds_b / dz))
depths_b = np.append(depths_b, (deptht[np.nanargmax(ds_b / dz)]
+ deptht[np.nanargmax(ds_b / dz) + 1])/2)
max_dsdz_b = np.append(max_dsdz_b, np.nanmax(ds_b / dz))
In [152]:
fig, ax = plt.subplots(figsize = (20, 5))
ax.plot(max_dsdz, 'o', label = 'CTD')
ax.plot(max_dsdz_base, 'ro', label = 'base case')
ax.plot(max_dsdz_a, 'yo', label = 'Test A')
ax.plot(max_dsdz_b, 'o', color = 'gold', label = 'Test B')
ax.set_xlabel('Station')
labels = [n//2 + 1 for n in range(36)]
plt.xticks(np.arange(36), labels);
ax.grid('on')
ax.set_title('Maximum ds / dz')
ax.legend();
print('base case bias = ' + str(-np.mean(max_dsdz) + np.mean(max_dsdz_base)))
print('test a bias = ' + str(-np.mean(max_dsdz) + np.mean(max_dsdz_a)))
print('test b bias = ' + str(-np.mean(max_dsdz) + np.mean(max_dsdz_b)))
In [153]:
fig, ax = plt.subplots(figsize = (20, 5))
ax.plot(depth_of_max_dsdz, 'o', label = 'CTD')
ax.plot(depths_base, 'ro', label = 'base case')
ax.plot(depths_a, 'yo', label = 'Test A')
ax.plot(depths_b, 'o', color = 'gold', label = 'Test B')
ax.set_xlabel('Station')
labels = [n//2 + 1 for n in range(36)]
plt.xticks(np.arange(36), labels);
ax.grid('on')
ax.set_title('Depth of maximum ds / dz')
ax.legend();
print('base case bias = ' + str(-np.mean(depth_of_max_dsdz) + np.mean(depths_base)))
print('test a bias = ' + str(-np.mean(depth_of_max_dsdz) + np.mean(depths_a)))
print('test b bias = ' + str(-np.mean(depth_of_max_dsdz) + np.mean(depths_b)))
In [126]:
colours = cm.tab20(np.linspace(0,1,19))
fig, ax = plt.subplots(figsize = (8,8))
for n in range(19):
station = stationdata.iloc[[n]]
cast = casts[station['Station'].values[0]]
ax.plot(station['LonDecDeg'].values[0], station['LatDecDeg'].values[0], 'o',
color = colours[n], label = 'Station ' + str(station['Station'].values[0]))
viz_tools.plot_coastline(ax, Bathymetry, coords = 'map')
viz_tools.set_aspect(ax)
ax.legend(loc = 'lower left')
ax.set_xlim(-124, -123)
ax.set_ylim(48.8, 49.5)
ax.set_title('Location of Stations');
In [154]:
lons = np.array([])
for n in range(19):
if n == 13:
lons = np.append(lons, stationdata['LonDecDeg'][n])
elif n == 14:
lons = np.append(lons, stationdata['LonDecDeg'][n])
else:
lons = np.append(lons, stationdata['LonDecDeg'][n])
lons = np.append(lons, stationdata['LonDecDeg'][n])
print(lons.shape)
lons
Out[154]:
In [157]:
fig, ax = plt.subplots(figsize = (20, 5))
ax.plot(lons[:18], max_dsdz[:18], 'o', label = 'CTD', color = 'darkblue', alpha = 0.3)
ax.plot(lons[:18], max_dsdz_base[:18], 'ro', label = 'base case', alpha = 0.3)
ax.plot(lons[:18], max_dsdz_a[:18], 'o' ,color = 'maroon', label = 'Test A', alpha = 0.3)
ax.plot(lons[:18], max_dsdz_b[:18], 'o', color = 'deeppink', label = 'Test B', alpha = 0.3)
ax.set_xlabel('Longitude')
ax.grid('on')
ax.set_title('Maximum ds / dz - April 10')
ax.legend();
print('base case bias = ' + str(-np.mean(max_dsdz[:18]) + np.mean(max_dsdz_base[:18])))
print('test a bias = ' + str(-np.mean(max_dsdz[:18]) + np.mean(max_dsdz_a[:18])))
print('test b bias = ' + str(-np.mean(max_dsdz[:18]) + np.mean(max_dsdz_b[:18])))
In [158]:
fig, ax = plt.subplots(figsize = (20, 5))
ax.plot(lons[18:], max_dsdz[18:], 'o', label = 'CTD', color = 'darkblue', alpha = 0.3)
ax.plot(lons[18:], max_dsdz_base[18:], 'ro', label = 'base case', alpha = 0.3)
ax.plot(lons[18:], max_dsdz_a[18:], 'o' ,color = 'maroon', label = 'Test A', alpha = 0.3)
ax.plot(lons[18:], max_dsdz_b[18:], 'o', color = 'deeppink', label = 'Test B', alpha = 0.3)
ax.set_xlabel('Longitude')
ax.grid('on')
ax.set_title('Maximum ds / dz - May 31')
ax.legend();
print('base case bias = ' + str(-np.mean(max_dsdz[18:]) + np.mean(max_dsdz_base[18:])))
print('test a bias = ' + str(-np.mean(max_dsdz[18:]) + np.mean(max_dsdz_a[18:])))
print('test b bias = ' + str(-np.mean(max_dsdz[18:]) + np.mean(max_dsdz_b[18:])))
In [159]:
fig, ax = plt.subplots(figsize = (20, 5))
ax.plot(lons[:18], depth_of_max_dsdz[:18], 'o',color = 'darkblue', label = 'CTD', alpha = 0.3)
ax.plot(lons[:18], depths_base[:18], 'ro', label = 'base case', alpha = 0.3)
ax.plot(lons[:18], depths_a[:18], 'o', color = 'maroon', label = 'Test A', alpha = 0.3)
ax.plot(lons[:18], depths_b[:18], 'o', color = 'hotpink', label = 'Test B', alpha = 0.3)
ax.set_xlabel('Longitude')
ax.grid('on')
ax.set_title('Depth of maximum ds / dz - April 10')
ax.legend();
print('base case bias = ' + str(-np.mean(depth_of_max_dsdz[:18]) + np.mean(depths_base[:18])))
print('test a bias = ' + str(-np.mean(depth_of_max_dsdz[:18]) + np.mean(depths_a[:18])))
print('test b bias = ' + str(-np.mean(depth_of_max_dsdz[:18]) + np.mean(depths_b[:18])))
In [160]:
fig, ax = plt.subplots(figsize = (20, 5))
ax.plot(lons[18:], depth_of_max_dsdz[18:], 'o',color = 'darkblue', label = 'CTD', alpha = 0.3)
ax.plot(lons[18:], depths_base[18:], 'ro', label = 'base case', alpha = 0.3)
ax.plot(lons[18:], depths_a[18:], 'o', color = 'maroon', label = 'Test A', alpha = 0.3)
ax.plot(lons[18:], depths_b[18:], 'o', color = 'hotpink', label = 'Test B', alpha = 0.3)
ax.set_xlabel('Longitude')
ax.grid('on')
ax.set_title('Depth of maximum ds / dz - May 31')
ax.legend();
print('base case bias = ' + str(-np.mean(depth_of_max_dsdz[18:]) + np.mean(depths_base[18:])))
print('test a bias = ' + str(-np.mean(depth_of_max_dsdz[18:]) + np.mean(depths_a[18:])))
print('test b bias = ' + str(-np.mean(depth_of_max_dsdz[18:]) + np.mean(depths_b[18:])))
In [ ]: