%matplotlib inline
import xarray as xr
import matplotlib.pyplot as plt
# from dask.dot import dot_graph
from dask.diagnostics import Profiler, ResourceProfiler, CacheProfiler
from dask.diagnostics import visualize
from bokeh.io import output_notebook
output_notebook()
from dask.distributed import Client
client = Client(scheduler_file='scheduler.json')
client
Client
|
Cluster
|
ds = xr.open_dataset('/cxfs/projects/usgs/hazards/cmgp/woodshole/aaretxabaleta/projects/GSB_tides_55nb/ocean_his_gsb_tides_55nb.nc',
engine='netcdf4', chunks={'eta_u':40, 'xi_rho':40}, decode_times=False)
import netCDF4
nc = netCDF4.Dataset('/cxfs/projects/usgs/hazards/cmgp/woodshole/aaretxabaleta/projects/GSB_tides_55nb/ocean_his_gsb_tides_55nb.nc')
nc.variables.keys()
odict_keys(['ntimes', 'ndtfast', 'dt', 'dtfast', 'dstart', 'nHIS', 'ndefHIS', 'nRST', 'Falpha', 'Fbeta', 'Fgamma', 'nl_visc2', 'LuvSponge', 'Akt_bak', 'Akv_bak', 'Akk_bak', 'Akp_bak', 'rdrg', 'rdrg2', 'Zob', 'Zos', 'gls_p', 'gls_m', 'gls_n', 'gls_cmu0', 'gls_c1', 'gls_c2', 'gls_c3m', 'gls_c3p', 'gls_sigk', 'gls_sigp', 'gls_Kmin', 'gls_Pmin', 'Charnok_alpha', 'Zos_hsig_alpha', 'sz_alpha', 'CrgBan_cw', 'Znudg', 'M2nudg', 'M3nudg', 'Tnudg', 'Tnudg_SSS', 'rho0', 'R0', 'Tcoef', 'Scoef', 'gamma2', 'LuvSrc', 'LwSrc', 'LtracerSrc', 'LsshCLM', 'Lm2CLM', 'Lm3CLM', 'LtracerCLM', 'LnudgeM2CLM', 'LnudgeM3CLM', 'LnudgeTCLM', 'spherical', 'xl', 'el', 'Vtransform', 'Vstretching', 'theta_s', 'theta_b', 'Tcline', 'hc', 'grid', 's_rho', 's_w', 'Cs_r', 'Cs_w', 'h', 'f', 'pm', 'pn', 'lon_rho', 'lat_rho', 'lon_u', 'lat_u', 'lon_v', 'lat_v', 'lon_psi', 'lat_psi', 'mask_rho', 'mask_u', 'mask_v', 'mask_psi', 'ocean_time', 'wetdry_mask_psi', 'wetdry_mask_rho', 'wetdry_mask_u', 'wetdry_mask_v', 'zeta', 'ubar', 'vbar', 'u', 'v', 'w', 'omega', 'temp', 'salt'])
dt, n, m = nc['zeta'].shape
print(dt,n,m)
1441 324 1542
t = nc['ocean_time'][:]/(3600*24)
t[:10]
array([0. , 0.02083333, 0.04166667, 0.0625 , 0.08333333, 0.10416667, 0.125 , 0.14583333, 0.16666667, 0.1875 ])
from utide import solve
from matplotlib.dates import date2num
import numpy as np
lat = 35.0147
#time = date2num(ds['ocean_time'].to_pydatetime())
%%time
acoef = solve(t, nc['zeta'][:,10,10], 0*t, lat,
trend=False, nodal=False, Rayleigh_min=0.95, method='ols',
conf_int='linear')
val = acoef['Lsmaj']
solve: matrix prep ... solution ... diagnostics ... done. CPU times: user 76 ms, sys: 15 ms, total: 91 ms Wall time: 94.7 ms
/home/rsignell/miniconda3/envs/pangeo_mpich/lib/python3.6/site-packages/utide/_solve.py:289: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions. To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`. m = np.linalg.lstsq(B, xraw)[0] # Model coefficients. /home/rsignell/miniconda3/envs/pangeo_mpich/lib/python3.6/site-packages/utide/confidence.py:237: RuntimeWarning: invalid value encountered in double_scalars varXv = Pvv[c] * varXv / den /home/rsignell/miniconda3/envs/pangeo_mpich/lib/python3.6/site-packages/utide/confidence.py:238: RuntimeWarning: invalid value encountered in double_scalars varYv = Pvv[c] * varYv / den
acoef.keys()
dict_keys(['name', 'aux', 'nR', 'nNR', 'nI', 'weights', 'Lsmaj', 'Lsmin', 'theta', 'g', 'umean', 'vmean', 'g_ci', 'Lsmaj_ci', 'Lsmin_ci', 'theta_ci', 'diagn'])
acoef['name']
array(['M2', 'N2', 'S2', 'K1', 'O1', 'ETA2', 'M6', 'M4', 'Q1', 'UPS1', 'OO1', 'J1', 'NO1', '2Q1', 'MO3', 'M3', 'SK3', 'MK3', 'MN4', 'MSF', 'MS4', 'S4', '2MK5', '2SK5', '3MK7', '2MN6', '2MS6', '2SM6', 'M8'], dtype=object)
import dask.array as da
from dask import delayed
usolve = delayed(solve)
z = nc['zeta']
%%time
nsub = 2
coefs = [usolve(t, z[:,j*nsub,i*nsub], t*0.0, lat,
trend=False, nodal=False, Rayleigh_min=0.95, method='ols',
conf_int='linear') for j in range(int(n/nsub)) for i in range(int(m/nsub))]
CPU times: user 3min 9s, sys: 4min 49s, total: 7min 59s Wall time: 8min 21s
arrays = [da.from_delayed(coef['Lsmaj'], # Construct a small Dask array
dtype=val.dtype, # for every lazy value
shape=val.shape)
for coef in coefs]
stack = da.stack(arrays, axis=0) # Stack all small Dask arrays into one
stack
dask.array<stack, shape=(124902, 29), dtype=float64, chunksize=(1, 29)>
%%time
amps = stack.compute()
CPU times: user 3min 20s, sys: 8.01 s, total: 3min 28s Wall time: 25min 29s
m2amp = amps[:,0].reshape((int(n/nsub),int(m/nsub)))
%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(figsize=(12,8))
plt.pcolormesh(m2amp)
plt.colorbar()
plt.title('M2 Elevation Amplitude');