from salishsea_tools import grid_tools, nc_tools, timeseries_tools, viz_tools
from mpl_toolkits.basemap import Basemap
from matplotlib import gridspec, patches, dates
from dateutil import parser
import numpy as np
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import pickle
import netCDF4 as nc
%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>''')
f0 = nc.Dataset('/data/vdo/MEOPAR/completed-runs/stats-runs/21sep14/SalishSea_1h_20140921_20140927_ptrc_T.nc')
f1 = nc.Dataset('/data/vdo/MEOPAR/completed-runs/stats-runs/15oct14/SalishSea_1h_20141015_20141025_ptrc_T.nc')
f2 = nc.Dataset('/data/vdo/MEOPAR/completed-runs/stats-runs/03dec14/SalishSea_1h_20141203_20141211_ptrc_T.nc')
f3 = nc.Dataset('/data/vdo/MEOPAR/completed-runs/stats-runs/29jun15/SalishSea_1h_20150629_20150706_ptrc_T.nc')
f4 = nc.Dataset('/data/vdo/MEOPAR/completed-runs/stats-runs/13jul15/SalishSea_1h_20150713_20150722_ptrc_T.nc')
f5 = nc.Dataset('/data/vdo/MEOPAR/completed-runs/stats-runs/27aug15/SalishSea_1h_20150827_20150903_ptrc_T.nc')
f6 = nc.Dataset('/data/vdo/MEOPAR/completed-runs/stats-runs/24jul16/SalishSea_1h_20160724_20160802_ptrc_T.nc')
f7 = nc.Dataset('/data/vdo/MEOPAR/completed-runs/stats-runs/02aug16/SalishSea_1h_20160802_20160828_ptrc_T.nc')
f8 = nc.Dataset('/data/vdo/MEOPAR/completed-runs/stats-runs/03apr17/SalishSea_1h_20170403_20170407_ptrc_T.nc')
f9 = nc.Dataset('/data/vdo/MEOPAR/completed-runs/stats-runs/10jun17/SalishSea_1h_20170610_20170616_ptrc_T.nc')
f10 = nc.Dataset('/data/vdo/MEOPAR/completed-runs/stats-runs/23jul17/SalishSea_1h_20170723_20170728_ptrc_T.nc')
f11 = nc.Dataset('/data/vdo/MEOPAR/completed-runs/stats-runs/12sep17/SalishSea_1h_20170912_20170916_ptrc_T.nc')
f0surface = f0.variables['mytracer3'][:,0,:,:]
f1surface = f1.variables['mytracer3'][:,0,...]
f2surface = f2.variables['mytracer3'][:,0,...]
f3surface = f3.variables['mytracer3'][:,0,...]
f4surface = f4.variables['mytracer3'][:,0,...]
f5surface = f5.variables['mytracer3'][:,0,...]
f6surface = f6.variables['mytracer3'][:,0,...]
f7surface = f7.variables['mytracer3'][:,0,...]
f8surface = f8.variables['mytracer3'][:,0,...]
f9surface = f9.variables['mytracer3'][:,0,...]
f10surface = f10.variables['mytracer3'][:,0,...]
f11surface = f11.variables['mytracer3'][:,0,...]
f9surface.shape
(168, 898, 398)
together = np.append(f0surface, f1surface, axis = 0)
for f in ([f2surface,f3surface,f4surface,f5surface,f6surface,f7surface,f8surface,
f9surface,f10surface,f11surface]):
together = np.append(together, f, axis = 0)
f.shape
(120, 898, 398)
together.shape
(2712, 898, 398)
mesh = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/mesh_mask201702.nc')
reshaped_f = np.ma.reshape(np.ma.masked_array(together,
mask = (1-mesh.variables['tmask'][:,0,...])
*np.ones((2712,1,1))), ((2712, -1)))
Z = np.reshape(reshaped_f.compressed(), ((2712, -1)))
A_primeb, sqrtLb, E_Tb = np.linalg.svd(Z - Z.mean(axis=0)[np.newaxis, :],
full_matrices=False)
Ab = A_primeb.dot(np.diag(sqrtLb))
PercentVarb = sqrtLb**2/(sqrtLb**2).sum()
Z_0b = Ab[:, 0, np.newaxis].dot(E_Tb[0, np.newaxis, :])
Z_1b = Ab[:, 1, np.newaxis].dot(E_Tb[1, np.newaxis, :])
fig, ax = plt.subplots(figsize = ((15,7)))
ax.plot(Ab[:,0].T)
ax.grid('on')
n = 0
for f in ([f0surface, f1surface, f2surface, f3surface, f4surface, f5surface,
f6surface, f7surface, f8surface, f9surface, f10surface, f11surface]):
j = f.shape[0]/2
ax.plot((j+n, j+n), (-50,30), '--')
print(j+n)
n = n + f.shape[0]
84.0 300.0 540.0 744.0 960.0 1176.0 1392.0 1836.0 2220.0 2364.0 2520.0 2652.0
from salishsea_tools import timeseries_tools, geo_tools, tidetools
grid = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
bathy, X, Y = tidetools.get_bathy_data(grid)
lats = grid.variables['nav_lat'][:]
lons = grid.variables['nav_lon'][:]
compressed_lats0 = np.ma.masked_array(lats,
mask = 1 -
mesh.variables['tmask'][0,0,...]).compressed()
compressed_lons0 = np.ma.masked_array(lons,
mask = 1
- mesh.variables['tmask'][0,0,...]).compressed()
Yinds = np.array([])
Xinds = np.array([])
for lon, lat in zip(compressed_lons0, compressed_lats0):
Yind, Xind = geo_tools.find_closest_model_point(lon, lat, X,Y, land_mask = bathy.mask)
Yinds = np.append(Yinds, Yind)
Xinds = np.append(Xinds, Xind)
import cmocean
PercentVarb[0]
0.4422418
for n in [84,300,540,744,960,1176,1392,1836,2220,2364,2520,2652]:
print(Z_0b[n,:].max(), Z_0b[n,:].min())
0.092359535 -0.014878162 0.02227008 -0.0035874788 0.27118152 -0.04368453 0.08791757 -0.0141626075 0.06470277 -0.010422944 0.20773636 -0.03346417 0.0086134765 -0.05347009 0.041771963 -0.25930884 0.25677288 -0.041363448 0.063548006 -0.0102369245 0.12326104 -0.019856075 0.046929296 -0.0075598224
for n,p in zip([84,300,540,744,960,1176,1392,1836,2220,2364,2520,2652],
[f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11]):
gridded = np.zeros((898,398))
for Yind, Xind, data in zip(Yinds, Xinds, Z_0b[n,:]):
gridded[int(Yind), int(Xind)] = data
fig, ax = plt.subplots(figsize = ((20,12)))
z = ax.pcolormesh(np.ma.masked_array(gridded,
mask = 1 - mesh.variables['tmask'][0,0,...]),
vmin = -0.25, vmax = 0.25, cmap = cmocean.cm.balance)
fig.colorbar(z, ax=ax)
middle = p.variables['time_counter'].shape[0]/2
ax.set_title('PC mode 1 ' + str(nc.num2date(p.variables['time_counter'][middle],
p.variables['time_counter'].units)))
viz_tools.set_aspect(ax);