Create CF-Compliant model data in "Sigma-over-Z" coordinates

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


/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)


(40, 1280, 2048)

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]:
<matplotlib.colorbar.Colorbar instance at 0x3d09878>

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


L=41
Ls=20

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])


137.869

In [12]:
zlev


Out[12]:
array([ -137.8687439 ,  -137.8687439 ,  -137.8687439 ,  -137.8687439 ,
        -137.8687439 ,  -137.8687439 ,  -137.8687439 ,  -137.8687439 ,
        -137.8687439 ,  -137.8687439 ,  -137.8687439 ,  -137.8687439 ,
        -137.8687439 ,  -137.8687439 ,  -137.8687439 ,  -137.8687439 ,
        -137.8687439 ,  -137.8687439 ,  -137.8687439 ,  -137.8687439 ,
        -165.03752136,  -197.3631897 ,  -235.82455444,  -281.58624268,
        -336.03393555,  -400.8163147 ,  -477.89498901,  -569.60394287,
        -678.72009277,  -808.54748535,  -963.01733398, -1146.80700684,
       -1365.48168945, -1625.6628418 , -1935.22888184, -2303.5534668 ,
       -2741.78955078, -3263.20678711, -3883.59375   , -4621.73632812,
       -5500.        ], dtype=float32)

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()


(1280, 2048)
Out[14]:
<matplotlib.colorbar.Colorbar instance at 0x3927758>

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


[ -137.8687439   -137.8687439   -137.8687439   -137.8687439   -137.8687439
  -137.8687439   -137.8687439   -137.8687439   -137.8687439   -137.8687439
  -137.8687439   -137.8687439   -137.8687439   -137.8687439   -137.8687439
  -137.8687439   -137.8687439   -137.8687439   -137.8687439   -151.453125
  -181.2003479   -216.59387207  -258.7053833   -308.81008911  -368.42510986
  -439.35565186  -523.74945068  -624.1619873   -743.63378906  -885.78240967
 -1054.91210938 -1256.14428711 -1495.57226562 -1780.44580078 -2119.39111328
 -2522.67138672 -3002.49804688 -3573.40039062 -4252.66503906 -5060.86816406]

In [18]:
print siglay


[-0.00362664 -0.01156828 -0.02101733 -0.03225989 -0.04563639 -0.06155188
 -0.08048827 -0.10301898 -0.12982622 -0.16172171 -0.19967128 -0.24482402
 -0.29854721 -0.3624675  -0.43852049 -0.52900904 -0.63667321 -0.76477301
 -0.91718733 -1.         -1.         -1.         -1.         -1.         -1.
 -1.         -1.         -1.         -1.         -1.         -1.         -1.
 -1.         -1.         -1.         -1.         -1.         -1.         -1.
 -1.        ]

In [19]:
print depth_c


137.869

In [20]:
print depth[jj,ii]


235.825

In [21]:
# point where depth < depth_c 
ii=1780
jj=1279
print depth[jj,ii]


5.0

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()


[ -1.81331895e-02  -5.78414202e-02  -1.05086632e-01  -1.61299422e-01
  -2.28181943e-01  -3.07759374e-01  -4.02441353e-01  -5.15094876e-01
  -6.49131060e-01  -8.08608532e-01  -9.98356402e-01  -1.22412014e+00
  -1.49273610e+00  -1.81233752e+00  -2.19260240e+00  -2.64504528e+00
  -3.18336606e+00  -3.82386494e+00  -4.58593655e+00  -1.51453125e+02
  -1.81200348e+02  -2.16593872e+02  -2.58705383e+02  -3.08810089e+02
  -3.68425110e+02  -4.39355652e+02  -5.23749451e+02  -6.24161987e+02
  -7.43633789e+02  -8.85782410e+02  -1.05491211e+03  -1.25614429e+03
  -1.49557227e+03  -1.78044580e+03  -2.11939111e+03  -2.52267139e+03
  -3.00249805e+03  -3.57340039e+03  -4.25266504e+03  -5.06086816e+03]

In [23]:
zij1[:nsiglay+1]
print nsiglay
print shape(temp)


19
(40, 1280, 2048)

In [24]:
print temp[nsiglay-1:nsiglay+1,jj,ii]
print zij1[nsiglay-1:nsiglay+1]
print depth[jj,ii]


[ 1.14973044  0.        ]
[  -4.58593655 -151.453125  ]
5.0

In [25]:
# point where depth > depth_c 
ii=400
jj=500
print depth[jj,ii]


235.825

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()


[  -0.5          -1.59490478   -2.89763236   -4.44762993   -6.29183149
   -8.48608017  -11.09681702  -14.20309734  -17.89897728  -22.29636765
  -27.52842903  -33.753582    -41.16032791  -49.97293854  -60.45827103
  -72.933815    -87.77733612 -105.43829346 -126.45146179 -151.453125
 -181.2003479  -216.59387207]

In [27]:
t = temp[:,jj,ii]
print t[21:23]
print zij1[21:23]


[ 20.81197548   0.        ]
[-216.59387207 -258.7053833 ]

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


0: angle / (unknown)                   (*ANONYMOUS*: 1280; *ANONYMOUS*: 2048)
1: dx / (unknown)                      (*ANONYMOUS*: 1280; *ANONYMOUS*: 2048)
2: dy / (unknown)                      (*ANONYMOUS*: 1280; *ANONYMOUS*: 2048)
3: sea_water_potential_temperature / (degC) (time: 1; *ANONYMOUS*: 40; *ANONYMOUS*: 1280; *ANONYMOUS*: 2048)

In [32]:
slice=cube[3]

In [33]:
print slice


sea_water_potential_temperature                  (time: 1; *ANONYMOUS*: 40; *ANONYMOUS*: 1280; *ANONYMOUS*: 2048)
     Dimension coordinates:
          time                                             x               -                -                  -
     Auxiliary coordinates:
          sea_surface_height                               x               -                x                  x
          siglay                                           -               x                -                  -
          zlay                                             -               x                -                  -
          depth                                            -               -                x                  x
          latitude                                         -               -                x                  x
          longitude                                        -               -                x                  x
     Derived coordinates:
          sea_surface_height_above_reference_ellipsoid     x               x                x                  x
     Scalar coordinates:
          depth_c: 137.869 m
          nsiglay: 19
     Attributes:
          Conventions: CF-1.6

In [34]:
z=slice.coord('sea_surface_height_above_reference_ellipsoid')

In [35]:
shape(z)


Out[35]:
(1, 40, 1280, 2048)

In [36]:
zij=z[0,:,jj,ii].points
print zij[ind]
plot(zij[ind],'o-')
grid()


[  -0.5          -1.59490478   -2.89763236   -4.44762993   -6.29183149
   -8.48608017  -11.09681702  -14.20309734  -17.89897728  -22.29636765
  -27.52842903  -33.753582    -41.16032791  -49.97293854  -60.45827103
  -72.933815    -87.77733612 -105.43829346 -126.45146179 -151.453125
 -181.2003479  -216.59387207]

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]: