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.
In [8]:
from pylab import *
import matplotlib.tri as Tri
import netCDF4
import datetime as dt
In [9]:
# DAP Data URL
url = 'http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST.nc'
# Open DAP
nci = netCDF4.Dataset(url)
nciv =nc.variables
nciv.keys()
Out[9]:
In [11]:
print ncv['u']
nt,nz,nele = shape(nciv['u'])
node = len(nciv['h'])
In [20]:
# OUTPUT: create NetCDF4 file with deflation on variables
url_out = '/usgs/data2/rsignell/models/fvcom/massbay_forecast_wbot.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=nciv['time'].units
ho.units=nciv['h'].units
lono.units=nciv['lon'].units
lato.units=nciv['lat'].units
lonco.units=nciv['lonc'].units
latco.units=nciv['latc'].units
ubot.units=nciv['u'].units
vbot.units=nciv['v'].units
bustr.units='N/m2'
bvstr.units='N/m2'
# write data with no time dimension
lonco[:]=nciv['lonc'][:]
latco[:]=nciv['latc'][:]
lono[:]=nciv['lon'][:]
lato[:]=nciv['lat'][:]
nvo[:]=nciv['nv'][:]
ho[:]=nciv['h'][:]
In [139]:
# take a look at the "metadata" for the variable "u"
print nc['u']
In [142]:
dtime = netCDF4.num2date(time_var[itime],time_var.units)
daystr = dtime.strftime('%Y-%b-%d %H:%M')
print daystr
In [143]:
tri = Tri.Triangulation(lon,lat, triangles=nv)
In [144]:
#get height of velocity points above bottom: zr = (bottom level - bottom layer) * water_depth
a,b=shape(nc['siglay'])
f=shape(nc['siglay'])
type(f)
Out[144]:
In [145]:
# specify bottom layer, handling case where there is just 1 layer in input file
if shape(nc['siglay'])[0]==1:
ilayer = 0
else:
ilayer = -1
In [146]:
def cdtoz0(cd=2.5e-3,zr=1.0):
kappa = 0.4
z0 = zr/(exp(kappa*ones_like(cd)/sqrt(cd)))
return z0
In [147]:
def z0tocd(z0=3.3546e-04,zr=1.0):
kappa = 0.4
cd=(kappa*ones_like(zr)/log(zr/z0))**2
return cd
In [148]:
# neither z0 or cd is saved in output, so just use canonical kb=0.5 cm
kb=0.005
z0=kb/30.
In [149]:
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 = sqrt(cd)*w
kappa = 0.4
ur = abs(ustar)/kappa*log(ones_like(zr)/z0)
wbot = w*ur/(abs(w)+1e-16)
return ustar,wbot
In [150]:
for itime in range(len(time_var)):
print '%4.1f percent done' % (float(itime)/len(time_var)*100.)
# get current at layer [0 = surface, -1 = bottom]
zr = (nc['siglay'][ilayer,:]-nc['siglev'][ilayer,:])*(nc['h'][:]+nc['zeta'][itime,:])
u = nc['u'][itime, ilayer, :]
v = nc['v'][itime, ilayer, :]
# average nodes to get bottom layer thicknesses at faces (velocity points)
zr_face = mean(zr[nv],axis=1)
w=u+1j*v # create complex velocity from components
ustar,wbot = w100(w,z0,zr_face) # compute bottom stress and velocity at 1 mab
nc_out['u'][itime,0,:]=ustar.real
nc_out['v'][itime,0,:]=ustar.imag
In [151]:
ncd.close()
ncd_out.close()
In [45]:
abs(array(w100()))
Out[45]:
In [46]:
#woods hole
levels=arange(-30,2,1)
ax = [-70.7, -70.6, 41.48, 41.55]
maxvel = 1.0
subsample = 2
In [59]:
#boston harbor
levels=arange(-34,2,1) # depth contours to plot
ax= [-70.97, -70.82, 42.25, 42.35] #
maxvel = 0.5
subsample = 3
In [60]:
# north shore
levels=arange(-20,4,1) # depth contours to plot
ax=[ -70.7874, -70.7702, 42.5548, 42.5671]
maxvel = 0.1
subsample = 1
In [61]:
# find velocity points in bounding box
ind = argwhere((lonc >= ax[0]) & (lonc <= ax[1]) & (latc >= ax[2]) & (latc <= ax[3]))
In [62]:
np.random.shuffle(ind)
Nvec = int(len(ind) / subsample)
idv = ind[:Nvec]
In [80]:
len(idv)
Out[80]:
In [77]:
ubot=wbot.real
vbot=wbot.imag
In [98]:
# tricontourf plot of water depth with vectors on top
figure(figsize=(18,10))
subplot(111,aspect=(1.0/cos(mean(latc)*pi/180.0)))
tricontourf(tri, -h,levels=levels,shading='faceted',cmap=plt.cm.gist_earth)
axis(ax)
gca().patch.set_facecolor('0.5')
cbar=colorbar()
cbar.set_label('Water Depth (m)', rotation=-90)
Q = quiver(lonc[idv],latc[idv],ubot[idv].flatten(),vbot[idv].flatten(),scale=1)
maxstr='%3.1f m/s' % maxvel
qk = quiverkey(Q,0.92,0.08,maxvel,maxstr,labelpos='W')
title('NECOFS near-bottom velocity (1 mab): %s UTC' % ( daystr));
In [100]:
# tricontourf plot of water depth with vectors on top
figure(figsize=(18,10))
subplot(111,aspect=(1.0/cos(mean(lat)*pi/180.0)))
tricontourf(tri, -h,levels=levels,shading='faceted',cmap=plt.cm.gist_earth)
axis(ax)
gca().patch.set_facecolor('0.5')
cbar=colorbar()
cbar.set_label('Water Depth (m)', rotation=-90)
Q = quiver(lonc[idv],latc[idv],ustar[idv].real.flatten(),ustar[idv].imag.flatten(),scale=0.05)
maxstr='%3.1f m/s' % maxvel
qk = quiverkey(Q,0.92,0.08,maxvel,maxstr,labelpos='W')
title('NECOFS ustart bottom stress (1 mab): %s UTC' % ( daystr));
In [23]:
# now here is the funny thing:
zeta=nc['zeta'][1,:]
h=nc['h'][:]
htot=h+zeta
In [24]:
# minimum of h (water depth) is -2:
h.min()
Out[24]:
In [25]:
# and the zeta values at these locations is 2.05
zeta[h==-2]
Out[25]:
In [26]:
# oh, I guess these are dry points waiting to be wetted, aren't they