ESMF Mesh file for unstructured grid

README.md

To run the script just type following:

python mesh_create.py --ifile data.nc

The script requires numpy, xarray, dask and also datetime modules. The required data file can be found in the following repository

https://github.com/uturuncoglu/data/blob/main/mesh/data.nc

mesh_create.py

import os, sys
import argparse
try:
    import numpy as np
except ImportError as ie:
    sys.stderr.write("""You are missing numpy module!""")
    sys.stderr.write(str(ie))
    exit(2)
try:
    import xarray as xr
except ImportError as ie:
    sys.stderr.write("""You are missing xarray module!""")
    sys.stderr.write(str(ie))
    exit(2)
try:
    import dask.array as da
except ImportError as ie:
    sys.stderr.write("""You are missing dask module!""")
    sys.stderr.write(str(ie))
    exit(2)
try:
    from datetime import datetime
except ImportError as ie:
    sys.stderr.write("""You are missing datetime module!""")
    sys.stderr.write(str(ie))
    exit(2)

def main(argv):
    """
    """
    # set defaults for command line arguments
    ifile = ''

    # read command line arguments
    parser = argparse.ArgumentParser()
    parser.add_argument('--ifile', help='Input file', required=True)
    args = parser.parse_args()

    if args.ifile:
        ifile = args.ifile

    # print out configuration
    print("Configuration:")
    print("ifile     = {}".format(ifile))

    # open file
    if os.path.isfile(ifile):
        ds = xr.open_dataset(ifile, mask_and_scale=False, decode_times=False)
    else:
        print('Input file could not find!')
        sys.exit(2)

    # populate required variables to write ESMF mesh file
    lon = ds['longitude']
    lat = ds['latitude'] 
    node_coords = da.stack([lon, lat], axis=1)
    elem_conn = ds['tri']
    elem_coords = np.array([da.mean(node_coords[v.values-1,:], axis=0) for v in elem_conn])

    # write data to file
    out = xr.Dataset()
    out['origGridDims'] = xr.DataArray(np.array(lon.shape, dtype=np.int32), 
                                       dims=('origGridRank'))
    out['nodeCoords']   = xr.DataArray(node_coords.compute_chunk_sizes(),
                                       dims=('nodeCount', 'coordDim'),
                                       attrs={'units': 'degrees'})
    out['elementConn']  = xr.DataArray(elem_conn, 
                                       dims=('elementCount', 'maxNodePElement'), 
                                       attrs={'long_name': 'Node indices that define the element connectivity'})
    out.elementConn.encoding = {'dtype': np.int32}
    out['numElementConn'] = xr.DataArray(3*np.ones(elem_conn.shape[0], dtype=np.int32),
                                         dims=('elementCount'),
                                         attrs={'long_name': 'Number of nodes per element'})
    out['centerCoords'] = xr.DataArray(elem_coords,
                                       dims=('elementCount', 'coordDim'),
                                       attrs={'units': 'degrees'})

    # force no '_FillValue' if not specified
    for v in out.variables:
        if '_FillValue' not in out[v].encoding:
            out[v].encoding['_FillValue'] = None

    # add global attributes
    out.attrs = {'title': 'ESMF unstructured grid file with {} dimension'.format('x'.join(list(map(str,lat.shape)))),
                 'created_by': os.path.basename(__file__),
                 'date_created': '{}'.format(datetime.now()),
                 'conventions': 'ESMFMESH',
                }

    # write output file
    ofile = 'mesh.nc'
    if ofile is not None:
        print('Writing {} ...'.format(ofile))
        out.to_netcdf(ofile)

if __name__== "__main__":
    main(sys.argv[1:])