Following the CF Conventions for Dimensionless Vertical Coordinates here: http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.6/apd.html
Here we are using sample global NCOM data from Paul Martin
In [1]:
cd /usgs/data2/rsignell/models/ncom
In [2]:
import netCDF4
In [3]:
# input file
nci = netCDF4.Dataset('t3d_glb8_3b_2013010100.nc')
In [4]:
tvar=nci.variables['Temperature']
In [5]:
fmin=tvar.pack_min
fmax=tvar.pack_max
nbits=tvar.pack_nbits
a=1./float(2**nbits-1)
In [6]:
levs=(nci.variables['level'][:]) - 1
In [7]:
levs=levs.astype(int)
In [8]:
temp= ones(shape(tvar))
print shape(temp)
Lz, ny, nx = shape(temp)
In [9]:
for lev in levs:
temp[lev,:,:] = nci.variables['Temperature'][lev,:,:]*(fmax[lev]-fmin[lev])*a + fmin[lev]
figure(figsize=(12,8))
pcolormesh(temp[0,::5,::5])
colorbar()
Out[9]:
In [10]:
# read vertical grid file
infile='ovgrd_1.D'
f = open(infile,'rb')
# seems to be an extra value at beginning, so read that as 'foo'
foo = fromfile(f,dtype='float32',count=1).byteswap()
al = fromfile(f,dtype='float32',count=1).byteswap()
als = fromfile(f,dtype='float32',count=1).byteswap()
L=int(al)
Ls=int(als)
zw = fromfile(f,dtype='float32',count=L).byteswap()
f.close()
print 'L=%d' % L
print 'Ls=%s' % Ls
In [11]:
depth_c = -zw[Ls-1]
print depth_c
sigma = copy(zw)
sigma[:Ls] = zw[:Ls]/depth_c
sigma[Ls:]= -1
zlev = zw
zlev[:Ls] = zw[Ls-1]
# average levels to layer centers
zlay = 0.5*(zlev[1:]+zlev[:-1])
siglay = 0.5*(sigma[1:]+sigma[:-1])
In [12]:
zlev
Out[12]:
In [13]:
# read horizontal grid file
infile='ohgrd_1.A'
f = open(infile,'rb')
lon = fromfile(f,dtype='float32',count=nx*ny).reshape((nx,ny),order='FORTRAN').byteswap().T
lat = fromfile(f,dtype='float32',count=nx*ny).reshape((nx,ny),order='FORTRAN').byteswap().T
dx = fromfile(f,dtype='float32',count=nx*ny).reshape((nx,ny),order='FORTRAN').byteswap().T
dy = fromfile(f,dtype='float32',count=nx*ny).reshape((nx,ny),order='FORTRAN').byteswap().T
h = fromfile(f,dtype='float32',count=nx*ny).reshape((nx,ny),order='FORTRAN').byteswap().T
ang = fromfile(f,dtype='float32',count=nx*ny).reshape((nx,ny),order='FORTRAN').byteswap().T
f.close()
In [14]:
depth = -h # depth positive
#depth = ma.masked_where(depth==0,depth)
print shape(depth)
figure(figsize=(12,8))
pcolormesh(depth[::5,::5])
colorbar()
Out[14]:
In [15]:
# create output file
nco = netCDF4.Dataset('ncom_sigma_over_z.nc','w', clobber=True)
nco.createDimension('level', Lz)
nco.createDimension('time', 1)
nco.createDimension('lon',nx)
nco.createDimension('lat',ny)
nco.Conventions='CF-1.6'
timeo = nco.createVariable('time', 'f4', ('time'))
# zlev(k) and sigma(k), and depth_c are the controlling variables for "ocean_sigma_z_coordinate"
zlevo = nco.createVariable('zlay', 'f4', ('level'))
sigmao = nco.createVariable('siglay', 'f4', ('level'))
depthc = nco.createVariable('depth_c', 'f4')
nsiglayo = nco.createVariable('nsiglay', 'i4')
lono = nco.createVariable('lon', 'f4', ('lat','lon'))
lato = nco.createVariable('lat', 'f4', ('lat','lon'))
depo = nco.createVariable('depth', 'f4', ('lat','lon'))
ango = nco.createVariable('angle', 'f4', ('lat','lon'))
dxo = nco.createVariable('dx', 'f4', ('lat','lon'))
dyo = nco.createVariable('dy', 'f4', ('lat','lon'))
zetao = nco.createVariable('zeta', 'f4', ('time', 'lat', 'lon'))
tempo = nco.createVariable('temp', 'f4', ('time', 'level', 'lat', 'lon'))
timeo.units = 'days since 2013-01-01 00:00'
timeo.standard_name='time'
tempo.units= 'degC'
tempo.standard_name='sea_water_potential_temperature'
# I think this is wrong, but currently required by Iris
#tempo.coordinates='time siglay lat lon'
tempo.coordinates='siglay zeta depth depth_c nsiglay zlay lat lon'
zetao.units='m'
zetao.standard_name='sea_surface_height'
zetao.coordinates='time lat lon'
depo.units = 'm'
depo.standard_name='depth'
depo.coordinates='lat lon'
lono.units= 'degrees_east'
lono.standard_name='longitude'
lato.units= 'degrees_north'
lato.standard_name='latitude'
sigmao.units='1'
sigmao.standard_name='ocean_sigma_z_coordinate'
sigmao.positive='up'
sigmao.formula_terms = 'sigma: siglay eta: zeta depth: depth depth_c: depth_c nsigma: nsiglay zlev: zlay'
depthc.units='m'
zlevo.units='m'
# write time
timeo[0]=0.0
# write zlev, sigma and depth_c
zlevo[:]=zlay
sigmao[:]=siglay
depthc[:]=depth_c
nsiglay = Ls - 1
nsiglayo[:]=nsiglay
# write depth
depo[:,:]=depth
# write lon,lat
lono[:,:]=lon
lato[:,:]=lat
ango[:,:]=ang
dxo[:,:]=dx
dyo[:,:]=dy
# write all zeros to zeta
zeta=zeros_like(depth)
zetao[0,:,:]=zeta
# write temperature
tempo[0,:,:,:]=temp
nco.close()
In [16]:
# sample point to look at z positions
ii=400
jj=500
In [17]:
print zlay
In [18]:
print siglay
In [19]:
print depth_c
In [20]:
print depth[jj,ii]
In [21]:
# point where depth < depth_c
ii=1780
jj=1279
print depth[jj,ii]
In [22]:
# plot vertical coordinate at a specified point
zij1 = siglay*(min(depth[jj,ii],depth_c))
zij1[nsiglay:] = zlay[nsiglay:]
print zij1
plot(zij1[:nsiglay],'ko-')
grid()
In [23]:
zij1[:nsiglay+1]
print nsiglay
print shape(temp)
In [24]:
print temp[nsiglay-1:nsiglay+1,jj,ii]
print zij1[nsiglay-1:nsiglay+1]
print depth[jj,ii]
In [25]:
# point where depth > depth_c
ii=400
jj=500
print depth[jj,ii]
In [26]:
# plot vertical coordinate at a specified point
zij1 = siglay*(min(depth[jj,ii],depth_c))
zij1[nsiglay:] = zlay[nsiglay:]
ind = where(zij1 > -depth[jj,ii])[0]
print zij1[ind]
plot(zij1[ind],'ko-')
grid()
In [27]:
t = temp[:,jj,ii]
print t[21:23]
print zij1[21:23]
Okay, this looks good: the last non-zero temperature point t[21] has a z value z[21]=-216 that is just shallower than the bottom (depth=235), and the first zero temperature value t[22] has a z value z[22]=-258 that is deeper than the bottom.
Now let's see what we get from Iris:
In [28]:
import iris
In [29]:
cubes = iris.load('http://geoport.whoi.edu/thredds/dodsC/usgs/data2/rsignell/models/ncom/ncom_sigma_over_z.nc')
In [30]:
cube = iris.load('ncom_sigma_over_z.nc')
In [31]:
print cube
In [32]:
slice=cube[3]
In [33]:
print slice
In [34]:
z=slice.coord('sea_surface_height_above_reference_ellipsoid')
In [35]:
shape(z)
Out[35]:
In [36]:
zij=z[0,:,jj,ii].points
print zij[ind]
plot(zij[ind],'o-')
grid()
Uh oh, it doesn't look like the sigma part of the sigma-over-z is working correctly.... we should not have zeros at the surface
Just for reference, here is the original Fortran code I received from Paul Martin at NRL to read the vertical and horizontal grid files
In [37]:
'''
c Dimensions are n = 2048, m = 1280, L = 41, Ls = 20.
real elon(n,m),alat(n,m),dx(n,m),dy(n,m),h(n,m),ang(n,m),zw(L)
c The switch from sigma to z-level occurs at interface = Ls
c L is the total number of interfaces (number of layers + 1).
c read horizontal grid fi
nrecl=n*m*4
open(199,file=infile,form='unformatted',status='old',
& access='direct',recl=nrecl)
read(199,rec=1) elon
read(199,rec=2) alat
read(199,rec=3) dx
read(199,rec=4) dy
read(199,rec=5) h
read(199,rec=6) ang
close(199)
c read vertical grid file.
infile='ovgrd_1.D'
open(199,file=infile,form='unformatted',status='old')
read(199) al,als,zw
close(199)
subroutine zwt_msk_sz(n,m,L,Ls, h,zw, zwt,amsk)
c subroutine to compute 3D array of interface depths (zwt) and
c 3D land-sea mask (amsk) for a sigma/z-level grid.
c n,m = horizontal grid dimensions in x and y.
c L = total number of layers +1.
c Ls = index of interface for switch from sigma to z-layers.
c this is the number of sigma layers +1.
c h = 2D array of static bottom depth + upward in m.
c h >= 0 at land pts, h < 0 at sea pts.
c zw = input 1D array of reference depths, + upward in m.
c these are the static depths from the surface to the
c layer interfaces on the z-level grid, and on the sigma
c grid in water of depth h <= zw(ls).
c zwt = returned 3D array of interface depths + upward in m.
c amsk = returned 3D land-sea mask defined at the grid-cell
c centers: =0 at land pts, =1 at sea pts.
c created 2013-06-11, Paul J Martin, NRL.
implicit none
c
c declare passed variables.
integer n,m,L,Ls
real h(n,m),zw(L),zwt(n,m,L),amsk(n,m,L)
c
c declare local temporary variables.
integer i,j,k
c
c define default values (used at land pts).
do k=1,L
do j=1,m
do i=1,n
zwt(i,j,k)=zw(k)
amsk(i,j,k)=0.0
enddo
enddo
enddo
c
c calculate the 3D array of depths and 3D land-sea mask at sea pts.
do j=1,m
do i=1,n
if (h(i,j) .lt. 0.0) then
if (h(i,j) .gt. zw(Ls)) then
do k=1,Ls
zwt(i,j,k)=h(i,j)*(zw(k)/zw(Ls))
if (k .lt. Ls) amsk(i,j,k)=1.0
enddo
else
do k=1,L-1
if (zw(k) .gt. h(i,j)) amsk(i,j,k)=1.0
enddo
endif
endif
enddo
enddo
c
return
end
''';
In [37]: