NetCDF Converter
The page below provides code for converting MITgcm binary output to netCDF format. This works for most cases, but will break if you give it data files that are surface-only fields, like surface ice thickness and ice area fraction. The code below has been extended with a command line interface and ability to convert 2D data fields as part of the gcmpy package, in particular gcmpy.netCDFbinary
. See the documentation here. This package adds another python package requirement: f90nml
, which is available on the package index (pip install f90nml
)
MITgcm dat to NetCDF converter
The code below provides a python script to convert MITgcm .data
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 install
ed, 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 (or before if the script complains) 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 13-15
# 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-10 and comment-out lines 13-15
# 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()