Notebook to prepare NEMO 3.4 Tide Files based on CONCEPTS 110 Tide Files
In [102]:
%matplotlib inline
from matplotlib import pylab
import netCDF4 as NC
from numpy import sin, cos, pi
from numpy import arange, zeros
Open CONCEPTS Tide Files & Bathymetry File
In [9]:
fT = NC.Dataset('../../input_files/SubDom_tide_elev.nc','r')
fU = NC.Dataset('../../input_files/SubDom_tide_ubar.nc','r')
fV = NC.Dataset('../../input_files/SubDom_tide_vbar.nc','r')
fB = NC.Dataset('../../input_files/SubDom_bathy_meter_NOBCchancomp.nc','r')
Extract Bathymetric Depths
In [144]:
Depth = fB.variables['Bathymetry']
depth = Depth[:]
# masking just causes problems, set depth to 0 on land
depth[depth.mask] = 0
print depth.shape
print depth[180,:]
(345, 398) [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 4.0 8.0 10.0 6.0 5.0 9.0 17.0 25.0 49.0 65.0 66.0 73.0 83.0 98.0 116.0 125.0 137.0 140.0 154.0 173.0 185.0 202.0 212.0 220.0 224.0 225.0 228.0 227.0 232.0 230.0 228.0 224.0 224.0 216.0 212.0 210.0 207.0 204.0 201.0 202.0 206.0 203.0 209.0 212.0 198.0 188.0 169.0 155.0 137.0 130.0 117.0 112.0 100.0 78.0 61.0 46.0 37.0 33.0 41.0 25.0 24.0 4.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 23.0 51.0 53.0 53.0 39.0 4.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.0 25.0 28.0 31.0 36.0 17.0 10.0 0.0 0.0 0.0 0.0 24.0 51.0 74.0 75.0 51.0 26.0 0.0 0.0 0.0 0.0 7.0 15.0 7.0 0.0 0.0 0.0 0.0 4.0 23.0 50.0 72.0 110.0 125.0 174.0 203.0 186.0 172.0 157.0 148.0 142.0 140.0 136.0 134.0 131.0 128.0 127.0 126.0 125.0 125.0 121.0 120.0 121.0 120.0 120.0 120.0 118.0 118.0 117.0 117.0 115.0 115.0 115.0 116.0 117.0 117.0 121.0 128.0 116.0 111.0 108.0 107.0 84.0 94.0 26.0 33.0 35.0 17.0 12.0 8.0 10.0 14.0 16.0 17.0 17.0 21.0 24.0 27.0 29.0 29.0 31.0 31.0 32.0 32.0 32.0 31.0 30.0 20.0 8.0 7.0 6.0 12.0 12.0 7.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
Find the depths on the u-grid and v-grid along the Northern Boundary. Complicated by the fact the Bathymetry array is masked.
In [146]:
# Northern Boundary
iindex_N = 345-1
jlast = 398-1
# u is off-set from T and bathymetry so find values in between
depthu_N = zeros((jlast+1), dtype=float)
depthu_N[0:jlast-1] = 0.5*(depth[iindex_N,0:jlast-1]+depth[iindex_N,1:jlast])
depthu_N[jlast] = depth[iindex_N,jlast]
# plot to check
pylab.plot (depthu_N,'x',depth[iindex_N,:],'r+')
pylab.xlim((175,205))
# now the V points, offset in the other direction
depthv_N = zeros((jlast+1), dtype=float)
depthv_N = 0.5*(depth[iindex_N,:]+depth[iindex_N-1,:])
In [147]:
# Western Boundary
ilast = 345-1
jindex_W = 1-1
# First the u-points
depthu_W = zeros((ilast+1), dtype=float)
depthu_W = 0.5*(depth[:,jindex_W]+depth[:,jindex_W+1])
# now V points
depthv_W = zeros((ilast+1), dtype=float)
depthv_W[0:ilast-1] = 0.5*(depth[0:ilast-1,jindex_W]+depth[1:ilast,jindex_W])
depthv_W[ilast] = depth[ilast,jindex_W]
# plot to check
pylab.plot (depthv_W,'x',depth[:,jindex_W],'r+')
pylab.xlim((250,275))
Out[147]:
(250, 275)
Extract Tides along the Northern and Western Boundaries
In [148]:
etaelev = fT.variables['M2_amp_elev']
umax = fU.variables['M2_amp_u']
vmax = fV.variables['M2_amp_v']
print etaelev.shape, umax.shape, vmax.shape
etaphase = fT.variables['M2_phi_elev']
uphase = fU.variables['M2_phi_u']
vphase = fV.variables['M2_phi_v']
etaelev_N = etaelev[iindex_N,:]
umax_N = umax[iindex_N,:]/depthu_N # divide by depth because CONCEPTS uses flux not current
vmax_N = vmax[iindex_N,:]/depthv_N
etaphase_N = etaphase[iindex_N,:]*pi/180.
uphase_N = uphase[iindex_N,:]*pi/180.
vphase_N = vphase[iindex_N,:]*pi/180.
etaZ1_N = etaelev_N*cos(etaphase_N)
etaZ2_N = etaelev_N*sin(etaphase_N)
uZ1_N = umax_N*cos(uphase_N)
uZ2_N = umax_N*sin(uphase_N)
vZ1_N = vmax_N*cos(vphase_N)
vZ2_N = vmax_N*sin(vphase_N)
print etaZ1_N.size
jindex_W = 1-1
etaelev_W = etaelev[:,jindex_W]
umax_W = umax[:,jindex_W]/depthu_W
vmax_W = vmax[:,jindex_W]/depthv_W
etaphase_W = etaphase[:,jindex_W]*pi/180.
uphase_W = uphase[:,jindex_W]*pi/180.
vphase_W = vphase[:,jindex_W]*pi/180.
etaZ1_W = etaelev_W*cos(etaphase_W)
etaZ2_W = etaelev_W*sin(etaphase_W)
uZ1_W = umax_W*cos(uphase_W)
uZ2_W = umax_W*sin(uphase_W)
vZ1_W = vmax_W*cos(vphase_W)
vZ2_W = vmax_W*sin(vphase_W)
# pad
uZ1_W[185:189+1] = uZ1_W[190]
(345, 398) (345, 398) (345, 398) 398
In [149]:
pylab.plot(umax_N[177:197],'b+',depth[344,177:197]/100.,'rx',umax_N[177:197]/depthu_N[177:197]*20,'g+')
Out[149]:
[<matplotlib.lines.Line2D at 0x108331d50>, <matplotlib.lines.Line2D at 0x108331690>, <matplotlib.lines.Line2D at 0x108331dd0>]
Test that I got the conversion correct
In [150]:
phi = arange(0,6.2,0.1)
eta_C = etaelev_W[200]*cos(phi-etaphase_W[200])
eta_N = etaZ1_W[200]*cos(phi)+etaZ2_W[200]*sin(phi)
pylab.plot(phi,eta_C,'+g',phi,eta_N,'-r')
Out[150]:
[<matplotlib.lines.Line2D at 0x10831f910>, <matplotlib.lines.Line2D at 0x108172e90>]
Create the NEMO 3.4 File
In [151]:
nemo = NC.Dataset('JPP_bdytide_M2_grid_T.nc', 'w', format='NETCDF3_CLASSIC')
nemo.description = 'Tide data from WebTide via JPP'
# dimensions
nemo.createDimension('xb', etaZ1_N.size+etaZ1_W.size)
nemo.createDimension('yb', 1)
# variables
xb = nemo.createVariable('xb', 'int32', ('xb',))
xb.units = 'non dim'
xb.longname = 'counter around boundary'
yb = nemo.createVariable('yb', 'int32', ('yb',))
yb.units = 'non dim'
nbidta = nemo.createVariable('nbidta', 'int32' , ('yb','xb'))
nbidta.units = 'non dim'
nbidta.longname = 'i grid position'
nbjdta = nemo.createVariable('nbjdta', 'int32' , ('yb','xb'))
nbjdta.units = 'non dim'
nbjdta.longname = 'j grid position'
nbrdta = nemo.createVariable('nbrdta', 'int32' , ('yb','xb'))
nbrdta.units = 'non dim'
z1 = nemo.createVariable('z1','float32',('yb','xb'))
z1.units = 'm'
z1.longname = 'tidal elevation: cosine'
z2 = nemo.createVariable('z2','float32',('yb','xb'))
z2.units = 'm'
z2.longname = 'tidal elevation: sine'
# data
total_length = etaZ1_N.size+etaZ1_W.size
xb[:] = arange(1, total_length+1)
yb[0] = 1
etaZ1_W[etaZ1_W.mask] = 0.
etaZ1_N[etaZ1_N.mask] = 0.
etaZ2_W[etaZ2_W.mask] = 0.
etaZ2_N[etaZ2_N.mask] = 0.
# West Boundary
nbidta[0,0:etaZ1_W.size] = 1
nbjdta[0,0:etaZ1_W.size] = arange(1,etaZ1_W.size+1)
z1[0,0:etaZ1_W.size] = etaZ1_W
z2[0,0:etaZ1_W.size] = etaZ2_W
# North Boundary
nbidta[0,etaZ1_W.size:] = arange(1,etaZ1_N.size+1)
nbjdta[0,etaZ1_W.size:] = iindex_N+1
z1[0,etaZ1_W.size:] = etaZ1_N
z2[0,etaZ2_W.size:] = etaZ2_N
# and both
nbrdta[:] = 1
pylab.plot(nbidta[0,:etaZ1_W.size],nbjdta[0,:etaZ1_W.size],'bo',
nbidta[0,etaZ1_W.size:],nbjdta[0,etaZ1_W.size:],'gx')
nemo.close()
In [152]:
u = NC.Dataset('JPP_bdytide_M2_grid_U.nc', 'w', format='NETCDF3_CLASSIC')
u.description = 'Tide data from WebTide via JPP'
#dimensions
u.createDimension('xb', etaZ1_N.size+etaZ1_W.size)
u.createDimension('yb', 1)
# variables
xb = u.createVariable('xb', 'int32', ('xb',))
xb.units = 'non dim'
xb.longname = 'counter around boundary'
yb = u.createVariable('yb', 'int32', ('yb',))
yb.units = 'non dim'
nbidta = u.createVariable('nbidta', 'int32' , ('yb','xb'))
nbidta.units = 'non dim'
nbidta.longname = 'i grid position'
nbjdta = u.createVariable('nbjdta', 'int32' , ('yb','xb'))
nbjdta.units = 'non dim'
nbjdta.longname = 'j grid position'
nbrdta = u.createVariable('nbrdta', 'int32' , ('yb','xb'))
nbrdta.units = 'non dim'
u1 = u.createVariable('u1','float32',('yb','xb'))
u1.units = 'm'
u1.longname = 'tidal x-velocity: cosine'
u2 = u.createVariable('u2','float32',('yb','xb'))
u2.units = 'm'
u2.longname = 'tidal x-velocity: sine'
# data
total_length = etaZ1_N.size+etaZ1_W.size
xb[:] = arange(1, total_length+1)
yb[0] = 1
uZ1_W[uZ1_W.mask] = 0.
uZ1_N[uZ1_N.mask] = 0.
uZ2_W[uZ2_W.mask] = 0.
uZ2_N[uZ2_N.mask] = 0.
# West Boundary
nbidta[0,0:etaZ1_W.size] = 1
nbjdta[0,0:etaZ1_W.size] = arange(1,etaZ1_W.size+1)
u1[0,0:etaZ1_W.size] = uZ1_W
u2[0,0:etaZ1_W.size] = uZ2_W
# North Boundary
nbidta[0,etaZ1_W.size:] = arange(1,etaZ1_N.size+1)
nbjdta[0,etaZ1_W.size:] = iindex_N+1
u1[0,etaZ1_W.size:] = uZ1_N
u2[0,etaZ2_W.size:] = uZ2_N
# and both
nbrdta[:] = 1
pylab.plot(nbidta[0,:etaZ1_W.size],nbjdta[0,:etaZ1_W.size],'bo',
nbidta[0,etaZ1_W.size:],nbjdta[0,etaZ1_W.size:],'gx')
u.close()
In [158]:
v = NC.Dataset('JPP_bdytide_M2_grid_V.nc', 'w', format='NETCDF3_CLASSIC')
v.description = 'Tide data from WebTide via JPP'
# dimensions
v.createDimension('xb', etaZ1_N.size+etaZ1_W.size)
v.createDimension('yb', 1)
# variables
xb = v.createVariable('xb', 'int32', ('xb',))
xb.units = 'non dim'
xb.longname = 'counter around boundary'
yb = v.createVariable('yb', 'int32', ('yb',))
yb.units = 'non dim'
nbidta = v.createVariable('nbidta', 'int32' , ('yb','xb'))
nbidta.units = 'non dim'
nbidta.longname = 'i grid position'
nbjdta = v.createVariable('nbjdta', 'int32' , ('yb','xb'))
nbjdta.units = 'non dim'
nbjdta.longname = 'j grid position'
nbrdta = v.createVariable('nbrdta', 'int32' , ('yb','xb'))
nbrdta.units = 'non dim'
v1 = v.createVariable('v1','float32',('yb','xb'))
v1.units = 'm'
v1.longname = 'tidal y-velocity: cosine'
v2 = v.createVariable('v2','float32',('yb','xb'))
v2.units = 'm'
v2.longname = 'tidal y-velocity: sine'
# data
total_length = etaZ1_N.size+etaZ1_W.size
xb[:] = arange(1, total_length+1)
yb[0] = 1
vZ1_W[vZ1_W.mask] = 0.
vZ1_N[vZ1_N.mask] = 0.
vZ2_W[vZ2_W.mask] = 0.
vZ2_N[vZ2_N.mask] = 0.
# West Boundary
nbidta[0,0:etaZ1_W.size] = 1
nbjdta[0,0:etaZ1_W.size] = arange(1,etaZ1_W.size+1)
v1[0,0:etaZ1_W.size] = vZ1_W
v2[0,0:etaZ1_W.size] = vZ2_W
# North Boundary
nbidta[0,etaZ1_W.size:] = arange(1,etaZ1_N.size+1)
nbjdta[0,etaZ1_W.size:] = iindex_N+1
v1[0,etaZ1_W.size:] = vZ1_N
v2[0,etaZ2_W.size:] = vZ2_N
# and both
nbrdta[:] = 1
pylab.plot(nbidta[0,:etaZ1_W.size],nbjdta[0,:etaZ1_W.size],'bo',
nbidta[0,etaZ1_W.size:],nbjdta[0,etaZ1_W.size:],'gx')
v.close()
In [142]:
In [ ]: