Demonstration using the NetCDF4-Python library to compute bottom stress and bottom velocity (1 meter above bottom) from a triangular grid ocean model (FVCOM) via OPeNDAP. The results are stored in a new NetCDF4 file.
NECOFS (Northeastern Coastal Ocean Forecast System) is run by groups at the University of Massachusetts Dartmouth and the Woods Hole Oceanographic Institution, led by Drs. C. Chen, R. C. Beardsley, G. Cowles and B. Rothschild. Funding is provided to run the model by the NOAA-led Integrated Ocean Observing System and the State of Massachusetts.
NECOFS is a coupled numerical model that uses nested weather models, a coastal ocean circulation model, and a wave model. The ocean model is a volume-mesh model with horizontal resolution that is finer in complicated regions. It is layered (not depth-averaged) and includes the effects of tides, winds, and varying water densities caused by temperature and salinity changes.
%matplotlib inline
import matplotlib.tri as Tri
import netCDF4
import datetime as dt
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
# DAP Data URL
#url = 'http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST.nc'
url = 'http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_GOM3_FORECAST.nc'
# input dataset from Open DAP
nci = netCDF4.Dataset(url)
print(nci['u'])
nt,nz,nele = nci['u'].shape
node = len(nci['h'])
time_var = nci['time']
<class 'netCDF4._netCDF4.Variable'> float32 u(time, siglay, nele) long_name: Eastward Water Velocity units: meters s-1 type: data missing_value: -999.0 field: ua, scalar coverage_content_type: modelResult standard_name: eastward_sea_water_velocity coordinates: time siglay latc lonc mesh: fvcom_mesh location: face unlimited dimensions: time current shape = (145, 40, 99137) filling off
# OUTPUT: create NetCDF4 file with deflation on variables
url_out = '/home/jovyan/out.nc'
nco = netCDF4.Dataset(url_out, 'w')
ncov = nco.variables
chunk_lon=512
chunk_time=1
sigdigits=6
nco.createDimension('nele', nele)
nco.createDimension('node', node)
nco.createDimension('three', 3)
nco.createDimension('time', None)
timeo = nco.createVariable('time', 'f4', ('time'))
ho = nco.createVariable('h', 'f4', ('node'))
nvo = nco.createVariable('nv', 'i4', ('three', 'nele'))
lonco = nco.createVariable('lonc', 'f4', ( 'nele'))
latco = nco.createVariable('latc', 'f4', ( 'nele'))
lono = nco.createVariable('lon', 'f4', ( 'node'))
lato = nco.createVariable('lat', 'f4', ( 'node'))
ubot = nco.createVariable('ubot', 'f4', ('time', 'nele'),
zlib=True,least_significant_digit=sigdigits,chunksizes=[chunk_time,chunk_lon])
vbot = nco.createVariable('vbot', 'f4', ('time', 'nele'),
zlib=True,least_significant_digit=sigdigits,chunksizes=[chunk_time,chunk_lon])
bustr = nco.createVariable('bustr', 'f4', ('time', 'nele'),
zlib=True,least_significant_digit=sigdigits,chunksizes=[chunk_time,chunk_lon])
bvstr = nco.createVariable('bvstr', 'f4', ('time', 'nele'),
zlib=True,least_significant_digit=sigdigits,chunksizes=[chunk_time,chunk_lon])
timeo.units=nci['time'].units
ho.units=nci['h'].units
lono.units=nci['lon'].units
lato.units=nci['lat'].units
lonco.units=nci['lonc'].units
latco.units=nci['latc'].units
ubot.units=nci['u'].units
vbot.units=nci['v'].units
bustr.units='N/m2'
bvstr.units='N/m2'
# write data with no time dimension
lonco[:]=nci['lonc'][:]
latco[:]=nci['latc'][:]
lono[:]=nci['lon'][:]
lato[:]=nci['lat'][:]
nvo[:]=nci['nv'][:]
ho[:]=nci['h'][:]
# take a look at the "metadata" for the variable "u"
print(nci['u'])
<class 'netCDF4._netCDF4.Variable'> float32 u(time, siglay, nele) long_name: Eastward Water Velocity units: meters s-1 type: data missing_value: -999.0 field: ua, scalar coverage_content_type: modelResult standard_name: eastward_sea_water_velocity coordinates: time siglay latc lonc mesh: fvcom_mesh location: face unlimited dimensions: time current shape = (145, 40, 99137) filling off
itime=0
dtime = netCDF4.num2date(time_var[itime],time_var.units)
daystr = dtime.strftime('%Y-%b-%d %H:%M')
print(daystr)
2018-Apr-21 00:00
lono[:]
array([-59.810688, -59.752678, -59.701324, ..., -73.95236 , -73.95679 , -73.948135], dtype=float32)
nv = nvo[:,:].T-1
tri = Tri.Triangulation(lono[:],lato[:], triangles=nv)
#get height of velocity points above bottom: zr = (bottom level - bottom layer) * water_depth
a,b=nci['siglay'].shape
f=nci['siglay'].shape
type(f)
tuple
# specify bottom layer, handling case where there is just 1 layer in input file
if nci['siglay'].shape[0]==1:
ilayer = 0
else:
ilayer = -1
def cdtoz0(cd=2.5e-3,zr=1.0):
kappa = 0.4
z0 = zr/(np.exp(kappa*np.ones_like(cd)/np.sqrt(cd)))
return z0
def z0tocd(z0=3.3546e-04,zr=1.0):
kappa = 0.4
cd=(kappa*np.ones_like(zr)/np.log(zr/z0))**2
return cd
# neither z0 or cd is saved in output, so just use canonical kb=0.5 cm
kb=0.005
z0=kb/30.
def w100(w=0.1+0j,z0=3.35e-04,zr=1):
""" Compute the velocity 1 meter above bottom and friction velocity
from velocity measured at height zr above bottom.
Keyword arguments
-----------------
w : east velocity component + j*north velocity component (m/s) [complex]
z0 : roughness height = kb/30 (m)
zr : height above bottom for input velocity "w" (m)
Returns
-------
ustar : friction velocity (m/s) [complex]
w : velocity 1 mab (m/s) [complex]
"""
cd = z0tocd(z0,zr)
ustar = np.sqrt(cd)*w
kappa = 0.4
ur = np.abs(ustar)/kappa*np.log10(np.ones_like(zr)/z0)
wbot = w*ur/(np.abs(w)+1e-16)
return ustar,wbot
#for itime in range(len(time_var)):
for itime in range(12):
print('{}% done'.format(float(itime)/len(time_var)*100.))
# get current at bottom layer
zr = 0.5*(nci['siglay'][ilayer,:]+1)*(nci['h'][:]+nci['zeta'][itime,:])
ur = nci['u'][itime, ilayer, :]
vr = nci['v'][itime, ilayer, :]
# average nodes to get bottom layer thicknesses at faces (velocity points)
zr_face = np.mean(zr[nv],axis=1)
w = ur +1j*vr # create complex velocity from components
ustar,wbot = w100(w,z0,zr_face) # compute bottom stress and velocity at 1 mab
bustr[itime,:]=ustar.real
bvstr[itime,:]=ustar.imag
ubot[itime,:]=wbot.real
vbot[itime,:]=wbot.imag
0.0% done 0.6896551724137931% done 1.3793103448275863% done 2.0689655172413794% done 2.7586206896551726% done 3.4482758620689653% done 4.137931034482759% done 4.827586206896552% done 5.517241379310345% done 6.206896551724138% done 6.896551724137931% done 7.586206896551724% done
#woods hole
levels = np.arange(-30,2,1)
ax = [-70.7, -70.6, 41.48, 41.55]
maxvel = 1.0
subsample = 2
#boston harbor
levels= np.arange(-34,2,1) # depth contours to plot
ax= [-70.97, -70.82, 42.25, 42.35] #
maxvel = 0.5
subsample = 3
# north shore
levels= np.arange(-20,4,1) # depth contours to plot
ax=[ -70.7874, -70.7702, 42.5548, 42.5671]
maxvel = 0.1
subsample = 1
# find velocity points in bounding box
ind = np.argwhere((lonco[:] >= ax[0]) & (lonco[:] <= ax[1]) & (latco[:] >= ax[2]) & (latco[:] <= ax[3]))
ind.shape
(506, 1)
np.random.shuffle(ind)
Nvec = int(len(ind) / subsample)
idv = ind[:Nvec]
ubot
<class 'netCDF4._netCDF4.Variable'> float32 ubot(time, nele) least_significant_digit: 6 units: meters s-1 unlimited dimensions: time current shape = (12, 99137) filling on, default _FillValue of 9.969209968386869e+36 used
itime = 0
idv=idv.flatten()
# tricontourf plot of water depth with vectors on top
plt.figure(figsize=(18,10))
plt.subplot(111,aspect=(1.0/np.cos(np.mean(lato[:])*np.pi/180.0)))
plt.tricontourf(tri, -ho[:],levels=levels,cmap=plt.cm.gist_earth)
plt.axis(ax)
plt.gca().patch.set_facecolor('0.5')
cbar = plt.colorbar()
cbar.set_label('Water Depth (m)', rotation=-90)
Q = plt.quiver(lonco[idv],latco[idv],ubot[itime,idv],vbot[itime,idv],scale=5)
maxstr='%3.1f m/s' % maxvel
qk = plt.quiverkey(Q,0.92,0.08,maxvel,maxstr,labelpos='W')
plt.title('NECOFS near-bottom velocity (1 mab): %s UTC' % ( daystr));
# tricontourf plot of water depth with vectors on top
plt.figure(figsize=(12,8))
plt.subplot(111,aspect=(1.0/np.cos(np.mean(lato[:])*np.pi/180.0)))
plt.tricontourf(tri, -ho[:],levels=levels,cmap=plt.cm.gist_earth,)
plt.axis(ax)
plt.gca().patch.set_facecolor('0.5')
cbar = plt.colorbar()
cbar.set_label('Water Depth (m)', rotation=-90)
Q = plt.quiver(lonco[idv],latco[idv],bustr[itime,idv],bvstr[itime,idv],scale=0.5)
maxstr=('{} m/s'.format(maxvel))
qk = plt.quiverkey(Q,0.92,0.08,maxvel,maxstr,labelpos='W')
plt.title('NECOFS ustart bottom stress (1 mab) {} UTC)'.format(daystr))
Text(0.5,1,'NECOFS ustart bottom stress (1 mab) 2018-Apr-21 00:00 UTC)')
nci.close()
nco.close()