NetCDF Converter

MITgcm dat to NetCDF converter

The code below provides a python script to convert MITgcm .dat files to NetCDF .nc files.

Note that using the script requires having the following python packages installed:

• MITgcmutils
• numpy
• netCDF4
• glob
• sys

'numpy' and 'netCDF4' can be pip installed, while 'glob' and 'sys' are likely provided with your python installation.

'MITgcmutils' is hidden within your MITgcm directory (MITgcm/utils/python/MITgcmutils). You will need to run the setup script (python setup.py install). Afterwards you will need to add it to your Python path. This can be done by adding the following line to your ~/.bashrc

export PYTHONPATH="${PYTHONPATH}:/path-to-MITgcm/utils/python/MITgcmutils/"

Usage

Converting specific output times

Comment-out lines 9-10 and uncomment lines 14-16

# Extract for a specific list of times
fields  = ['Rho', 'T']
indices = ['0000017040']

fields specifies which fields (file prefixes) should be converted, while indices specifies which output times should be converted.

Converting all output times

Uncomment lines 9-11 and comment-out lines 14-16

# Extract for all times
fields  = ['T', 'Rho']
names = glob.glob(fields[0] + '.*.data')
indices = [nm[len(fields[0])+1:-5] for nm in names]

fields specifies which fields (file prefixes) should be converted, while the remaining lines automatically detect which outputs exist for those fields.

Code

import MITgcmutils as mgu
import numpy as np
from netCDF4 import Dataset
import glob, sys

# Step 1: Load data

# Extract for a specific list of times
#fields  = ['Rho', 'T']
#indices = ['0000017040']

# Extract for all times
fields  = ['T', 'Rho']
names = glob.glob(fields[0] + '.*.data')
indices = [nm[len(fields[0])+1:-5] for nm in names]

print("Preparing to convert {0:d} outputs to netcdf.".format(len(indices)))

# Manually extract the grids
Rho = mgu.rdmds('Rho.0000000000')

Nx = Rho.shape[2]
Ny = Rho.shape[1]
Nz = Rho.shape[0]

print("(Nx, Ny, Nz) = ({0:d}, {1:d}, {2:d})".format(Nx, Ny, Nz))

XC = mgu.rdmds('XG')
YC = mgu.rdmds('YG')
RC = mgu.rdmds('RC')

depth = mgu.rdmds('Depth')

myX = np.tile(XC.reshape(1, Ny, Nx), (Nz, 1, 1))
myY = np.tile(YC.reshape(1, Ny, Nx), (Nz, 1, 1))
myZ = np.tile(RC.reshape(Nz, 1, 1), [1, Ny, Nx])

myDepth = np.tile(depth.reshape(1, Ny, Nx), (Nz, 1, 1))

myInds = np.zeros(depth.shape)

# Handle the topography
cntP = 0
cntN = 0
for jj in range(Ny):
    for ii in range(Nx):
        zs = myZ[:,jj,ii]
        myInds[jj,ii] = np.sum(zs >= -depth[jj,ii])

        # Any z points below the topography are mapped to the lowest
        #  above-topography point
        if int(myInds[jj,ii]) < Nz:
            zs[int(myInds[jj,ii]):] = zs[int(myInds[jj,ii])]
        if int(myInds[jj,ii]) == 0:
            cntP += 1
        if int(myInds[jj,ii]) == Nz:
            cntN += 1

        myZ[:,jj,ii] = zs

print("{0:d} of {1:d} points ({2:.3g}%) are purely topographic".format(cntP, Nx*Ny, 100.*cntP/(Nx*Ny)))
print("{0:d} of {1:d} points ({2:.3g}%) have no topography".format(cntN, Nx*Ny, 100.*cntN/(Nx*Ny)))

# Step 2: Create netcdf

for index in indices:
    print("Processing index {0}".format(index))

    # Create file
    fp = Dataset('output_{0}.nc'.format(index), 'w', format='NETCDF4')
    
    # Create dimension objects
    x_dim = fp.createDimension('x',Nx)
    y_dim = fp.createDimension('y',Ny)
    z_dim = fp.createDimension('z',Nz)
    t_dim = fp.createDimension('time',1)

    # Create and write grids
    x_grid  = fp.createVariable('x',    np.float64, ('x',))
    y_grid  = fp.createVariable('y',    np.float64, ('y',))
    z_grid  = fp.createVariable('z',    np.float64, ('z',))
    zc_grid = fp.createVariable('zc',   np.float64, ('z','y','x'))
    t_grid  = fp.createVariable('time', np.float64, ('time',))

    x_grid[:]      = myX[0,0,:]
    y_grid[:]      = myY[0,:,0]
    z_grid[:]      = np.arange(Nz)
    zc_grid[:,:,:] = myZ
    t_grid[0]      = int(index)

    # Create and write fields
    field_objs = []
    for f in fields:
        f_var = fp.createVariable(f, np.float64, ('z','y','x'))
        data = mgu.rdmds(f + '.' + index)

        dat_min = 0*np.min(data[myZ >= -myDepth])

        # Screw with the data to make topography work nicely
        for jj in range(Ny):
            for ii in range(Nx):
                if (int(myInds[jj,ii]) == 0):
                    data[:,jj,ii] = dat_min
                elif int(myInds[jj,ii]) < Nz:
                    data[int(myInds[jj,ii]):,jj,ii] = data[int(myInds[jj,ii]),jj,ii]

        f_var[:,:,:] = data

    fp.close()