Notebook

Testing parallel tidal analysis with Utide and Dask Delayed¶

We have a 3D (time,lat,lon) cube of water level data and we want to run a tidal analysis function ("solve" from Utide) at each lon,lat grid point (or every other lat,lon grid point). This is embarrassingly parallel, and we could use other python parallel approaches, but can we use Dask Delayed?

In [1]:

%matplotlib inline

import xarray as xr
import matplotlib.pyplot as plt

from dask.diagnostics import Profiler, ResourceProfiler, CacheProfiler
from dask.diagnostics import visualize

from bokeh.io import output_notebook
output_notebook()

from utide import solve
import numpy as np
import warnings

Loading BokehJS ...

In [3]:

from dask.distributed import Client, progress
client = Client(scheduler_file='/home/rsignell/scheduler.json')
client

Out[3]:

Client

Scheduler: tcp://10.12.0.103:8786
Dashboard: http://10.12.0.103:8787/status

Cluster

Workers: 9
Cores: 90
Memory: 81.00 GB

Open Zarr data with Xarray¶

In [5]:

%%time
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)

CPU times: user 61 ms, sys: 6 ms, total: 67 ms
Wall time: 66.3 ms

In [8]:

ds.info()

xarray.Dataset {
dimensions:
	eta_psi = 323 ;
	eta_rho = 324 ;
	eta_u = 324 ;
	eta_v = 323 ;
	ocean_time = 1441 ;
	s_rho = 7 ;
	s_w = 8 ;
	tracer = 2 ;
	xi_psi = 1541 ;
	xi_rho = 1542 ;
	xi_u = 1541 ;
	xi_v = 1542 ;

variables:
	int32 ntimes() ;
		ntimes:long_name = number of long time-steps ;
	int32 ndtfast() ;
		ndtfast:long_name = number of short time-steps ;
	float64 dt() ;
		dt:long_name = size of long time-steps ;
		dt:units = second ;
	float64 dtfast() ;
		dtfast:long_name = size of short time-steps ;
		dtfast:units = second ;
	float64 dstart() ;
		dstart:long_name = time stamp assigned to model initilization ;
		dstart:units = days since 0001-01-01 00:00:00 ;
		dstart:calendar = julian ;
	int32 nHIS() ;
		nHIS:long_name = number of time-steps between history records ;
	int32 ndefHIS() ;
		ndefHIS:long_name = number of time-steps between the creation of history files ;
	int32 nRST() ;
		nRST:long_name = number of time-steps between restart records ;
		nRST:cycle = only latest two records are maintained ;
	float64 Falpha() ;
		Falpha:long_name = Power-law shape barotropic filter parameter ;
	float64 Fbeta() ;
		Fbeta:long_name = Power-law shape barotropic filter parameter ;
	float64 Fgamma() ;
		Fgamma:long_name = Power-law shape barotropic filter parameter ;
	float64 nl_visc2() ;
		nl_visc2:long_name = nonlinear model Laplacian mixing coefficient for momentum ;
		nl_visc2:units = meter2 second-1 ;
	int32 LuvSponge() ;
		LuvSponge:long_name = horizontal viscosity sponge activation switch ;
		LuvSponge:flag_values = [0 1] ;
		LuvSponge:flag_meanings = .FALSE. .TRUE. ;
	float64 Akt_bak(tracer) ;
		Akt_bak:long_name = background vertical mixing coefficient for tracers ;
		Akt_bak:units = meter2 second-1 ;
	float64 Akv_bak() ;
		Akv_bak:long_name = background vertical mixing coefficient for momentum ;
		Akv_bak:units = meter2 second-1 ;
	float64 Akk_bak() ;
		Akk_bak:long_name = background vertical mixing coefficient for turbulent energy ;
		Akk_bak:units = meter2 second-1 ;
	float64 Akp_bak() ;
		Akp_bak:long_name = background vertical mixing coefficient for length scale ;
		Akp_bak:units = meter2 second-1 ;
	float64 rdrg() ;
		rdrg:long_name = linear drag coefficient ;
		rdrg:units = meter second-1 ;
	float64 rdrg2() ;
		rdrg2:long_name = quadratic drag coefficient ;
	float64 Zob() ;
		Zob:long_name = bottom roughness ;
		Zob:units = meter ;
	float64 Zos() ;
		Zos:long_name = surface roughness ;
		Zos:units = meter ;
	float64 gls_p() ;
		gls_p:long_name = stability exponent ;
	float64 gls_m() ;
		gls_m:long_name = turbulent kinetic energy exponent ;
	float64 gls_n() ;
		gls_n:long_name = turbulent length scale exponent ;
	float64 gls_cmu0() ;
		gls_cmu0:long_name = stability coefficient ;
	float64 gls_c1() ;
		gls_c1:long_name = shear production coefficient ;
	float64 gls_c2() ;
		gls_c2:long_name = dissipation coefficient ;
	float64 gls_c3m() ;
		gls_c3m:long_name = buoyancy production coefficient (minus) ;
	float64 gls_c3p() ;
		gls_c3p:long_name = buoyancy production coefficient (plus) ;
	float64 gls_sigk() ;
		gls_sigk:long_name = constant Schmidt number for TKE ;
	float64 gls_sigp() ;
		gls_sigp:long_name = constant Schmidt number for PSI ;
	float64 gls_Kmin() ;
		gls_Kmin:long_name = minimum value of specific turbulent kinetic energy ;
	float64 gls_Pmin() ;
		gls_Pmin:long_name = minimum Value of dissipation ;
	float64 Charnok_alpha() ;
		Charnok_alpha:long_name = Charnok factor for surface roughness ;
	float64 Zos_hsig_alpha() ;
		Zos_hsig_alpha:long_name = wave amplitude factor for surface roughness ;
	float64 sz_alpha() ;
		sz_alpha:long_name = surface flux from wave dissipation ;
	float64 CrgBan_cw() ;
		CrgBan_cw:long_name = surface flux due to Craig and Banner wave breaking ;
	float64 Znudg() ;
		Znudg:long_name = free-surface nudging/relaxation inverse time scale ;
		Znudg:units = day-1 ;
	float64 M2nudg() ;
		M2nudg:long_name = 2D momentum nudging/relaxation inverse time scale ;
		M2nudg:units = day-1 ;
	float64 M3nudg() ;
		M3nudg:long_name = 3D momentum nudging/relaxation inverse time scale ;
		M3nudg:units = day-1 ;
	float64 Tnudg(tracer) ;
		Tnudg:long_name = Tracers nudging/relaxation inverse time scale ;
		Tnudg:units = day-1 ;
	float64 Tnudg_SSS() ;
		Tnudg_SSS:long_name = SSS nudging/relaxation inverse time scale ;
		Tnudg_SSS:units = day-1 ;
	float64 rho0() ;
		rho0:long_name = mean density used in Boussinesq approximation ;
		rho0:units = kilogram meter-3 ;
	float64 R0() ;
		R0:long_name = background density used in linear equation of state ;
		R0:units = kilogram meter-3 ;
	float64 Tcoef() ;
		Tcoef:long_name = thermal expansion coefficient ;
		Tcoef:units = Celsius-1 ;
	float64 Scoef() ;
		Scoef:long_name = Saline contraction coefficient ;
	float64 gamma2() ;
		gamma2:long_name = slipperiness parameter ;
	int32 LuvSrc() ;
		LuvSrc:long_name = momentum point sources and sink activation switch ;
		LuvSrc:flag_values = [0 1] ;
		LuvSrc:flag_meanings = .FALSE. .TRUE. ;
	int32 LwSrc() ;
		LwSrc:long_name = mass point sources and sink activation switch ;
		LwSrc:flag_values = [0 1] ;
		LwSrc:flag_meanings = .FALSE. .TRUE. ;
	int32 LtracerSrc(tracer) ;
		LtracerSrc:long_name = tracer point sources and sink activation switch ;
		LtracerSrc:flag_values = [0 1] ;
		LtracerSrc:flag_meanings = .FALSE. .TRUE. ;
	int32 LsshCLM() ;
		LsshCLM:long_name = sea surface height climatology processing switch ;
		LsshCLM:flag_values = [0 1] ;
		LsshCLM:flag_meanings = .FALSE. .TRUE. ;
	int32 Lm2CLM() ;
		Lm2CLM:long_name = 2D momentum climatology processing switch ;
		Lm2CLM:flag_values = [0 1] ;
		Lm2CLM:flag_meanings = .FALSE. .TRUE. ;
	int32 Lm3CLM() ;
		Lm3CLM:long_name = 3D momentum climatology processing switch ;
		Lm3CLM:flag_values = [0 1] ;
		Lm3CLM:flag_meanings = .FALSE. .TRUE. ;
	int32 LtracerCLM(tracer) ;
		LtracerCLM:long_name = tracer climatology processing switch ;
		LtracerCLM:flag_values = [0 1] ;
		LtracerCLM:flag_meanings = .FALSE. .TRUE. ;
	int32 LnudgeM2CLM() ;
		LnudgeM2CLM:long_name = 2D momentum climatology nudging activation switch ;
		LnudgeM2CLM:flag_values = [0 1] ;
		LnudgeM2CLM:flag_meanings = .FALSE. .TRUE. ;
	int32 LnudgeM3CLM() ;
		LnudgeM3CLM:long_name = 3D momentum climatology nudging activation switch ;
		LnudgeM3CLM:flag_values = [0 1] ;
		LnudgeM3CLM:flag_meanings = .FALSE. .TRUE. ;
	int32 LnudgeTCLM(tracer) ;
		LnudgeTCLM:long_name = tracer climatology nudging activation switch ;
		LnudgeTCLM:flag_values = [0 1] ;
		LnudgeTCLM:flag_meanings = .FALSE. .TRUE. ;
	int32 spherical() ;
		spherical:long_name = grid type logical switch ;
		spherical:flag_values = [0 1] ;
		spherical:flag_meanings = Cartesian spherical ;
	float64 xl() ;
		xl:long_name = domain length in the XI-direction ;
		xl:units = meter ;
	float64 el() ;
		el:long_name = domain length in the ETA-direction ;
		el:units = meter ;
	int32 Vtransform() ;
		Vtransform:long_name = vertical terrain-following transformation equation ;
	int32 Vstretching() ;
		Vstretching:long_name = vertical terrain-following stretching function ;
	float64 theta_s() ;
		theta_s:long_name = S-coordinate surface control parameter ;
	float64 theta_b() ;
		theta_b:long_name = S-coordinate bottom control parameter ;
	float64 Tcline() ;
		Tcline:long_name = S-coordinate surface/bottom layer width ;
		Tcline:units = meter ;
	float64 hc() ;
		hc:long_name = S-coordinate parameter, critical depth ;
		hc:units = meter ;
	int32 grid() ;
		grid:cf_role = grid_topology ;
		grid:topology_dimension = 2 ;
		grid:node_dimensions = xi_psi eta_psi ;
		grid:face_dimensions = xi_rho: xi_psi (padding: both) eta_rho: eta_psi (padding: both) ;
		grid:edge1_dimensions = xi_u: xi_psi eta_u: eta_psi (padding: both) ;
		grid:edge2_dimensions = xi_v: xi_psi (padding: both) eta_v: eta_psi ;
		grid:node_coordinates = lon_psi lat_psi ;
		grid:face_coordinates = lon_rho lat_rho ;
		grid:edge1_coordinates = lon_u lat_u ;
		grid:edge2_coordinates = lon_v lat_v ;
		grid:vertical_dimensions = s_rho: s_w (padding: none) ;
	float64 s_rho(s_rho) ;
		s_rho:long_name = S-coordinate at RHO-points ;
		s_rho:valid_min = -1.0 ;
		s_rho:valid_max = 0.0 ;
		s_rho:positive = up ;
		s_rho:standard_name = ocean_s_coordinate_g2 ;
		s_rho:formula_terms = s: s_rho C: Cs_r eta: zeta depth: h depth_c: hc ;
		s_rho:field = s_rho, scalar ;
	float64 s_w(s_w) ;
		s_w:long_name = S-coordinate at W-points ;
		s_w:valid_min = -1.0 ;
		s_w:valid_max = 0.0 ;
		s_w:positive = up ;
		s_w:standard_name = ocean_s_coordinate_g2 ;
		s_w:formula_terms = s: s_w C: Cs_w eta: zeta depth: h depth_c: hc ;
		s_w:field = s_w, scalar ;
	float64 Cs_r(s_rho) ;
		Cs_r:long_name = S-coordinate stretching curves at RHO-points ;
		Cs_r:valid_min = -1.0 ;
		Cs_r:valid_max = 0.0 ;
		Cs_r:field = Cs_r, scalar ;
	float64 Cs_w(s_w) ;
		Cs_w:long_name = S-coordinate stretching curves at W-points ;
		Cs_w:valid_min = -1.0 ;
		Cs_w:valid_max = 0.0 ;
		Cs_w:field = Cs_w, scalar ;
	float64 h(eta_rho, xi_rho) ;
		h:long_name = bathymetry at RHO-points ;
		h:units = meter ;
		h:grid = grid ;
		h:location = face ;
		h:field = bath, scalar ;
	float64 f(eta_rho, xi_rho) ;
		f:long_name = Coriolis parameter at RHO-points ;
		f:units = second-1 ;
		f:grid = grid ;
		f:location = face ;
		f:field = coriolis, scalar ;
	float64 pm(eta_rho, xi_rho) ;
		pm:long_name = curvilinear coordinate metric in XI ;
		pm:units = meter-1 ;
		pm:grid = grid ;
		pm:location = face ;
		pm:field = pm, scalar ;
	float64 pn(eta_rho, xi_rho) ;
		pn:long_name = curvilinear coordinate metric in ETA ;
		pn:units = meter-1 ;
		pn:grid = grid ;
		pn:location = face ;
		pn:field = pn, scalar ;
	float64 lon_rho(eta_rho, xi_rho) ;
		lon_rho:long_name = longitude of RHO-points ;
		lon_rho:units = degree_east ;
		lon_rho:standard_name = longitude ;
		lon_rho:field = lon_rho, scalar ;
	float64 lat_rho(eta_rho, xi_rho) ;
		lat_rho:long_name = latitude of RHO-points ;
		lat_rho:units = degree_north ;
		lat_rho:standard_name = latitude ;
		lat_rho:field = lat_rho, scalar ;
	float64 lon_u(eta_u, xi_u) ;
		lon_u:long_name = longitude of U-points ;
		lon_u:units = degree_east ;
		lon_u:standard_name = longitude ;
		lon_u:field = lon_u, scalar ;
	float64 lat_u(eta_u, xi_u) ;
		lat_u:long_name = latitude of U-points ;
		lat_u:units = degree_north ;
		lat_u:standard_name = latitude ;
		lat_u:field = lat_u, scalar ;
	float64 lon_v(eta_v, xi_v) ;
		lon_v:long_name = longitude of V-points ;
		lon_v:units = degree_east ;
		lon_v:standard_name = longitude ;
		lon_v:field = lon_v, scalar ;
	float64 lat_v(eta_v, xi_v) ;
		lat_v:long_name = latitude of V-points ;
		lat_v:units = degree_north ;
		lat_v:standard_name = latitude ;
		lat_v:field = lat_v, scalar ;
	float64 lon_psi(eta_psi, xi_psi) ;
		lon_psi:long_name = longitude of PSI-points ;
		lon_psi:units = degree_east ;
		lon_psi:standard_name = longitude ;
		lon_psi:field = lon_psi, scalar ;
	float64 lat_psi(eta_psi, xi_psi) ;
		lat_psi:long_name = latitude of PSI-points ;
		lat_psi:units = degree_north ;
		lat_psi:standard_name = latitude ;
		lat_psi:field = lat_psi, scalar ;
	float64 mask_rho(eta_rho, xi_rho) ;
		mask_rho:long_name = mask on RHO-points ;
		mask_rho:flag_values = [0. 1.] ;
		mask_rho:flag_meanings = land water ;
		mask_rho:grid = grid ;
		mask_rho:location = face ;
	float64 mask_u(eta_u, xi_u) ;
		mask_u:long_name = mask on U-points ;
		mask_u:flag_values = [0. 1.] ;
		mask_u:flag_meanings = land water ;
		mask_u:grid = grid ;
		mask_u:location = edge1 ;
	float64 mask_v(eta_v, xi_v) ;
		mask_v:long_name = mask on V-points ;
		mask_v:flag_values = [0. 1.] ;
		mask_v:flag_meanings = land water ;
		mask_v:grid = grid ;
		mask_v:location = edge2 ;
	float64 mask_psi(eta_psi, xi_psi) ;
		mask_psi:long_name = mask on psi-points ;
		mask_psi:flag_values = [0. 1.] ;
		mask_psi:flag_meanings = land water ;
		mask_psi:grid = grid ;
		mask_psi:location = node ;
	float64 ocean_time(ocean_time) ;
		ocean_time:long_name = time since initialization ;
		ocean_time:units = seconds since 0001-01-01 00:00:00 ;
		ocean_time:calendar = julian ;
		ocean_time:field = time, scalar, series ;
	float32 wetdry_mask_psi(ocean_time, eta_psi, xi_psi) ;
		wetdry_mask_psi:long_name = wet/dry mask on PSI-points ;
		wetdry_mask_psi:flag_values = [0. 1.] ;
		wetdry_mask_psi:flag_meanings = land water ;
		wetdry_mask_psi:time = ocean_time ;
		wetdry_mask_psi:grid = grid ;
		wetdry_mask_psi:location = node ;
		wetdry_mask_psi:field = wetdry_mask_psi, scalar, series ;
	float32 wetdry_mask_rho(ocean_time, eta_rho, xi_rho) ;
		wetdry_mask_rho:long_name = wet/dry mask on RHO-points ;
		wetdry_mask_rho:flag_values = [0. 1.] ;
		wetdry_mask_rho:flag_meanings = land water ;
		wetdry_mask_rho:time = ocean_time ;
		wetdry_mask_rho:grid = grid ;
		wetdry_mask_rho:location = face ;
		wetdry_mask_rho:field = wetdry_mask_rho, scalar, series ;
	float32 wetdry_mask_u(ocean_time, eta_u, xi_u) ;
		wetdry_mask_u:long_name = wet/dry mask on U-points ;
		wetdry_mask_u:flag_values = [0. 1.] ;
		wetdry_mask_u:flag_meanings = land water ;
		wetdry_mask_u:time = ocean_time ;
		wetdry_mask_u:grid = grid ;
		wetdry_mask_u:location = edge1 ;
		wetdry_mask_u:field = wetdry_mask_u, scalar, series ;
	float32 wetdry_mask_v(ocean_time, eta_v, xi_v) ;
		wetdry_mask_v:long_name = wet/dry mask on V-points ;
		wetdry_mask_v:flag_values = [0. 1.] ;
		wetdry_mask_v:flag_meanings = land water ;
		wetdry_mask_v:time = ocean_time ;
		wetdry_mask_v:grid = grid ;
		wetdry_mask_v:location = edge2 ;
		wetdry_mask_v:field = wetdry_mask_v, scalar, series ;
	float32 zeta(ocean_time, eta_rho, xi_rho) ;
		zeta:long_name = free-surface ;
		zeta:units = meter ;
		zeta:time = ocean_time ;
		zeta:grid = grid ;
		zeta:location = face ;
		zeta:field = free-surface, scalar, series ;
	float32 ubar(ocean_time, eta_u, xi_u) ;
		ubar:long_name = vertically integrated u-momentum component ;
		ubar:units = meter second-1 ;
		ubar:time = ocean_time ;
		ubar:grid = grid ;
		ubar:location = edge1 ;
		ubar:field = ubar-velocity, scalar, series ;
	float32 vbar(ocean_time, eta_v, xi_v) ;
		vbar:long_name = vertically integrated v-momentum component ;
		vbar:units = meter second-1 ;
		vbar:time = ocean_time ;
		vbar:grid = grid ;
		vbar:location = edge2 ;
		vbar:field = vbar-velocity, scalar, series ;
	float32 u(ocean_time, s_rho, eta_u, xi_u) ;
		u:long_name = u-momentum component ;
		u:units = meter second-1 ;
		u:time = ocean_time ;
		u:grid = grid ;
		u:location = edge1 ;
		u:field = u-velocity, scalar, series ;
	float32 v(ocean_time, s_rho, eta_v, xi_v) ;
		v:long_name = v-momentum component ;
		v:units = meter second-1 ;
		v:time = ocean_time ;
		v:grid = grid ;
		v:location = edge2 ;
		v:field = v-velocity, scalar, series ;
	float32 w(ocean_time, s_w, eta_rho, xi_rho) ;
		w:long_name = vertical momentum component ;
		w:units = meter second-1 ;
		w:time = ocean_time ;
		w:standard_name = upward_sea_water_velocity ;
		w:grid = grid ;
		w:location = face ;
		w:field = w-velocity, scalar, series ;
	float32 omega(ocean_time, s_w, eta_rho, xi_rho) ;
		omega:long_name = S-coordinate vertical momentum component ;
		omega:units = meter second-1 ;
		omega:time = ocean_time ;
		omega:grid = grid ;
		omega:location = face ;
		omega:field = omega, scalar, series ;
	float32 temp(ocean_time, s_rho, eta_rho, xi_rho) ;
		temp:long_name = potential temperature ;
		temp:units = Celsius ;
		temp:time = ocean_time ;
		temp:grid = grid ;
		temp:location = face ;
		temp:field = temperature, scalar, series ;
	float32 salt(ocean_time, s_rho, eta_rho, xi_rho) ;
		salt:long_name = salinity ;
		salt:time = ocean_time ;
		salt:grid = grid ;
		salt:location = face ;
		salt:field = salinity, scalar, series ;

// global attributes:
	:file = ocean_his_gsb_tides_55nb.nc ;
	:format = netCDF-3 64bit offset file ;
	:Conventions = CF-1.4 ;
	:type = ROMS/TOMS history file ;
	:title = Great south Bay ;
	:var_info = varinfo.dat ;
	:rst_file = ocean_rst.nc ;
	:his_file = ocean_his_gsb_tides_55nb.nc ;
	:grd_file = ../grids/GSB_55nb.nc ;
	:frc_file_01 = ../forcings/tide_forc_GSB_55.nc ;
	:script_file =  ;
	:NLM_LBC = 
EDGE:  WEST   SOUTH  EAST   NORTH  
zeta:  Cha    Cha    Cha    Cha    
ubar:  Fla    Fla    Fla    Fla    
vbar:  Fla    Fla    Fla    Fla    
u:     Gra    Gra    Gra    Gra    
v:     Gra    Gra    Gra    Gra    
temp:  Gra    Gra    Gra    Gra    
salt:  Gra    Gra    Gra    Gra    
tke:   Gra    Gra    Gra    Gra ;
	:svn_url = https:://myroms.org/svn/src ;
	:svn_rev =  ;
	:code_dir = /cxfs/projects/usgs/hazards/cmgp/woodshole/aaretxabaleta/models/COAWST ;
	:header_dir = /cxfs/projects/usgs/hazards/cmgp/woodshole/aaretxabaleta/projects/GSB_tides_55nb ;
	:header_file = gsb.h ;
	:os = Linux ;
	:cpu = x86_64 ;
	:compiler_system = ifort ;
	:compiler_command = /opt/intel/impi/5.0.1.035/intel64/bin/mpif90 ;
	:compiler_flags = -heap-arrays -fp-model precise -ip -O3 -xW -free ;
	:tiling = 036x010 ;
	:history = ROMS/TOMS, Version 3.7, Monday - February 26, 2018 - 10:23:23 AM ;
	:ana_file = ROMS/Functionals/ana_btflux.h, ROMS/Functionals/ana_fsobc.h, ROMS/Functionals/ana_initial.h, ROMS/Functionals/ana_m2obc.h, ROMS/Functionals/ana_smflux.h, ROMS/Functionals/ana_stflux.h ;
	:CPP_options = GSB, ADD_FSOBC, ADD_M2OBC, ANA_BSFLUX, ANA_BTFLUX, ANA_FSOBC, ANA_INITIAL, ANA_M2OBC, ANA_SMFLUX, ANA_SSFLUX, ANA_STFLUX, ASSUMED_SHAPE, DJ_GRADPS, DOUBLE_PRECISION, GLS_MIXING, KANTHA_CLAYSON, MASKING, MIX_S_UV, MPI, NONLINEAR, !NONLIN_EOS, N2S2_HORAVG, POWER_LAW, PROFILE, K_GSCHEME, RAMP_TIDES, !RST_SINGLE, SOLVE3D, SSH_TIDES, TS_C4HADVECTION, TS_C4VADVECTION, TS_FIXED, UV_ADV, UV_COR, UV_U3HADVECTION, UV_C4VADVECTION, UV_LOGDRAG, UV_TIDES, UV_VIS2, VAR_RHO_2D, WET_DRY ;
}

In [ ]:

#client.get_versions(check=True)

In [9]:

dt, n, m = ds['zeta'].shape
print(dt,n,m)

1441 324 1542

In [10]:

ds['zeta']

Out[10]:

<xarray.DataArray 'zeta' (ocean_time: 1441, eta_rho: 324, xi_rho: 1542)>
dask.array<shape=(1441, 324, 1542), dtype=float32, chunksize=(1441, 324, 40)>
Coordinates:
    lon_rho     (eta_rho, xi_rho) float64 dask.array<shape=(324, 1542), chunksize=(324, 40)>
    lat_rho     (eta_rho, xi_rho) float64 dask.array<shape=(324, 1542), chunksize=(324, 40)>
  * ocean_time  (ocean_time) float64 0.0 1.8e+03 3.6e+03 5.4e+03 7.2e+03 ...
Dimensions without coordinates: eta_rho, xi_rho
Attributes:
    long_name:  free-surface
    units:      meter
    time:       ocean_time
    grid:       grid
    location:   face
    field:      free-surface, scalar, series

In [11]:

nsub=4

In [12]:

%time z = ds['zeta'][:,::nsub,::nsub].load().values

CPU times: user 770 ms, sys: 236 ms, total: 1.01 s
Wall time: 2min 7s

In [13]:

# convert time to days
t = ds['ocean_time'].values/(3600.*24)

In [14]:

# nominal latitude for tide calcs
lat = 40.7

In [15]:

%%time
# analyze tides at single cell as a test
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    acoef = solve(t=t, u=z[:,20,20], v=None, lat=lat, 
               trend=False, nodal=False, Rayleigh_min=0.95, method='ols',
               conf_int='linear', verbose=False)

CPU times: user 523 ms, sys: 79 ms, total: 602 ms
Wall time: 916 ms

In [16]:

# M2 amplitude
acoef['A'][0]

Out[16]:

0.6128407747097121

In [17]:

import dask.array as da
from dask import delayed

In [18]:

usolve = delayed(solve)

In [19]:

kk,jj,ii = z.shape

In [20]:

%%time 
# set up the Dask delayed task list
coefs = [usolve(t=t, u=z[:,j,i], v=None, lat=lat, 
   trend=False, nodal=False, verbose=False, Rayleigh_min=0.95, method='ols',
   conf_int='linear') for j in range(jj) for i in range(ii)]

CPU times: user 32.2 s, sys: 591 ms, total: 32.8 s
Wall time: 32.3 s

In [21]:

# compute tidal analysis (parallel)
%time total = delayed(coefs).compute()

CPU times: user 33.7 s, sys: 2.66 s, total: 36.4 s
Wall time: 7min 37s

%%time # compute tidal analysis in a regular loop (serial) with warnings.catch_warnings(): warnings.simplefilter("ignore") bcoef = [solve(t=t, u=z[:,j,i], v=None, lat=lat, trend=False, nodal=False, Rayleigh_min=0.95, method='ols', conf_int='linear', verbose=False) for j in range(jj) for i in range(ii)]

In [22]:

m2amp = [f['A'][0] for f in total]

In [23]:

len(m2amp)

Out[23]:

In [24]:

m2amp = np.array(m2amp).reshape((jj,ii))

In [25]:

%matplotlib inline
import matplotlib.pyplot as plt

In [26]:

plt.figure(figsize=(12,8))
plt.pcolormesh(m2amp)
plt.colorbar()
plt.title('M2 Elevation Amplitude');

In [ ]: