A Generic Example of Code to Convert Regular Gridded NetCDF with Winds/Pressures to NWS13 NetCDF, and also the structure of the code that does the same for sets of ascii NWS12 OWI WIN/PRE domains

generic_nws13_example.py

import netCDF4
import xarray as xr
import numpy as np
from glob import glob


params = ["U10", "V10", "PSFC", "lon", "lat", "time"]

def xrreaddumb(fpathpattern, ilatdim="lati", ilondim="long", 
               ilatcoord="latitude", iloncoord="longitude",
               iuwnd="uwnd", ivwnd="vwnd", ipres="msl"):
    fpaths = glob(fpathpattern)
    ds = xr.open_mfdataset(fpaths, 
                           combine="nested", concat_dim="time",
                           )
    # Rename input dimension, coordinate and wind/pres variable names
    ds = ds.rename_dims({ 
                               ilatdim: "yi", 
                               ilondim: "xi"})
    ds = ds.rename({iloncoord: "lon",
                    ilatcoord: "lat",
                    iuwnd: "U10",
                    ivwnd: "V10",
                    ipres: "PSFC",
                    })

    ds = ds.set_coords(("time", "lat", "lon"))
    
    ds["lon"].attrs = {
        "units": "degrees_east",
        "standard_name": "longitude",
        "axis": "X",
        "coordinates": "time lat lon",
    }
    ds["lat"].attrs = {
        "units": "degrees_north",
        "standard_name": "latitude",
        "axis": "Y",
        "coordinates": "time lat lon",
    }
    ds["time"].attrs = {"axis": "T", "coordinates": "time"}
    ds["U10"].attrs = {"units": "m s-1", "coordinates": "time lat lon"}
    ds["V10"].attrs = {"units": "m s-1", "coordinates": "time lat lon"}
    ds["PSFC"].attrs = {"units": "mb", "coordinates": "time lat lon"}
    ds.encoding = {}
    return ds[params]


def add_global_nws13(
    output,
    group_order,
    institution="Your Institution",
    conventions=["CF-1.6", "OWI-NWS13"],
):
    with netCDF4.Dataset(output, "a") as nc:
        nc.setncattr("group_order", " ".join(group_order))
        nc.setncattr("institution", institution)
        nc.setncattr("conventions", " ".join(conventions))


def to_nws13_group(output, group, rank, wind=None, pre=None, mode="w"):
    """Write win/pre data out to netcdf adcirc NWS13 formatted group.

    Args:
        output (str): Output netcdf filepath.
        group (str): Group name.
        rank (int): Group rank.
        wind (xarray.Dataset): Dataset with U10, V10 (and possibly PSFC), and coordinates.
        pre (xarray.Dataset): Dataset with PSFC varaible. Defaults to None.
        mode (str): Set to "a" to append to exisiting file, "w" will overwrite. Defaults to "w".

    Returns:
        None

    Examples:
        >>> mode = "w"
        >>> rank = 1
        >>> outgroups = ingroups
        >>> for inwin, inpre, group in zip(inwinfiles, inprefiles, ingroups):
        >>>     wind = xr.concat(
        >>>         list(iter_file_xrdataset(inwin, ftype="win")),
        >>>         dim="time",
        >>>         coords="different",
        >>>         data_vars="different",
        >>>     )
        >>>     to_nws13_group(output, group, rank, wind=wind, mode=mode)
        >>>     mode = "a"
        >>>     del wind
        >>>     pres = xr.concat(
        >>>         list(iter_file_xrdataset(inpre, ftype="pre")),
        >>>         dim="time",
        >>>         coords="different",
        >>>         data_vars="different",
        >>>     )
        >>>     to_nws13_group(output, group, rank, pre=pres, mode=mode)
        >>>     rank += 1
        >>>     del pres

    """
    if wind is not None:
        if "WS" in wind:
            del wind["WS"]
        wind.attrs["rank"] = np.int32(rank)
        wind.time.encoding = {
            "units": "minutes since 1990-01-01 01:00:00",
            "dtype": np.int,
        }
        wind[["U10", "V10"]].to_netcdf(output, group=group, mode=mode)
    if pre is not None:
        pre.time.encoding = {
            "units": "minutes since 1990-01-01 01:00:00",
            "dtype": np.int,
        }
        pre["PSFC"].to_netcdf(output, group=group, mode="a")


def dumb2nws13(inncpatterns, ingroups, output, mode="w", verbose=False):
    if mode == "a":
        data = netCDF4.Dataset(output)
        group_order = data.group_order.split()
        rank = len(group_order) + 1
        outgroups = group_order + ingroups
        data.close()
    else:
        rank = 1
        outgroups = ingroups
    for inwin, group in zip(inncpatterns, ingroups):
        if verbose:
            print(f"Reading {group}, Rank {rank} NetCDF, File Pattern: {inwin}")
        ds = xrreaddumb(inwin)
        if verbose:
            print(f"Writing WIN {group}, Rank {rank} to {output} with mode={mode}")
        to_nws13_group(output, group, rank, 
                       pre=ds, 
                       wind=ds,
                       mode=mode)
        rank += 1
        mode = "a"
    if verbose:
        print(
            "Writing global attributes, group_order='{0}'".format(" ".join(outgroups))
        )
    add_global_nws13(output, outgroups)


if __name__ == "__main__":
    patterns = [
               "NWS13/Wind_ExtendedSmoothT.nc",
               ]
    groups = ["Main"]
    output = "NWS13/test4_nws13.nc"
    dumb2nws13(patterns, groups, output, mode="w", verbose=True)

partial_winpre2nws13.py

"""usage: winpre_to_nws13.py [-h] -w WIN [WIN ...] -p PRE [PRE ...] -o OUTPUT
                          [-g GROUP [GROUP ...]] [-a] [-v]

Convert a series of OWI WIN/PRE files to an ADCIRC NWS13 Netcdf file.

optional arguments:
  -h, --help            show this help message and exit
  -w WIN [WIN ...], --win WIN [WIN ...]
                        Input WIN file(s)
  -p PRE [PRE ...], --pre PRE [PRE ...]
                        Input PRE file(s)
  -o OUTPUT, --output OUTPUT
                        Output file path
  -g GROUP [GROUP ...], --group GROUP [GROUP ...]
                        Input group name(s)
  -a, --append          If set, do not clobber the existing file, but append groups to it.
  -v, --verbose         Verbose log output

"""
import owiWinPre
import xarray as xr
import numpy as np
import os
import netCDF4
import argparse


# Fancy command line argument parsing setup
parser = argparse.ArgumentParser(
    description="Convert a series of OWI WIN/PRE files to an ADCIRC NWS13 Netcdf file.",
    formatter_class=argparse.RawTextHelpFormatter,
)
parser.add_argument(
    "-w", "--win", type=str, nargs="+", help="Input WIN file(s)", required=True
)
parser.add_argument(
    "-p", "--pre", type=str, nargs="+", help="Input PRE file(s)", required=True
)
parser.add_argument("-o", "--output", type=str, help="Output file path", required=True)
parser.add_argument(
    "-g",
    "--group",
    type=str,
    nargs="+",
    help="Input group name(s)",
    required=False,
    default=["Main"],
)
parser.add_argument(
    "-a",
    "--append",
    action="store_true",
    help="If set, do not clobber the existing file, but append groups to it.",
    required=False,
    default=False,
)
parser.add_argument(
    "-v",
    "--verbose",
    action="store_true",
    help="Verbose log output",
    required=False,
    default=False,
)


def to_nws13_group(output, group, rank, wind=None, pre=None, mode="w"):
    """Write win/pre data out to netcdf adcirc NWS13 formatted group.

    Args:
        output (str): Output netcdf filepath.
        group (str): Group name.
        rank (int): Group rank.
        wind (xarray.Dataset): Dataset with U10, V10 (and possibly PSFC), and coordinates.
        pre (xarray.Dataset): Dataset with PSFC varaible. Defaults to None.
        mode (str): Set to "a" to append to exisiting file, "w" will overwrite. Defaults to "w".

    Returns:
        None

    Examples:
        >>> mode = "w"
        >>> rank = 1
        >>> outgroups = ingroups
        >>> for inwin, inpre, group in zip(inwinfiles, inprefiles, ingroups):
        >>>     wind = xr.concat(
        >>>         list(iter_file_xrdataset(inwin, ftype="win")),
        >>>         dim="time",
        >>>         coords="different",
        >>>         data_vars="different",
        >>>     )
        >>>     to_nws13_group(output, group, rank, wind=wind, mode=mode)
        >>>     mode = "a"
        >>>     del wind
        >>>     pres = xr.concat(
        >>>         list(iter_file_xrdataset(inpre, ftype="pre")),
        >>>         dim="time",
        >>>         coords="different",
        >>>         data_vars="different",
        >>>     )
        >>>     to_nws13_group(output, group, rank, pre=pres, mode=mode)
        >>>     rank += 1
        >>>     del pres

    """
    if wind is not None:
        if "WS" in wind:
            del wind["WS"]
        wind.attrs["rank"] = np.int32(rank)
        wind.time.encoding = {
            "units": "minutes since 1990-01-01 01:00:00",
            "dtype": np.int,
        }
        wind.to_netcdf(output, group=group, mode=mode)
    if pre is not None:
        pre.time.encoding = {
            "units": "minutes since 1990-01-01 01:00:00",
            "dtype": np.int,
        }
        pre["PSFC"].to_netcdf(output, group=group, mode="a")


def iter_file_xrdataset(fpath, ftype=None):
    """Given a WIN/PRE filepath, return generator/iterable for timesteps that
    outputs xarray.Datasets.

    Args:
        fpath (str): Path to single WIN/PRE file.
        ftype (str): Force the filetype to either "win" or "pre". Defaults to None.

    Yields:
        xarray.Dataset: For timestep

    Examples:
        >>> wind = xr.concat(
        >>>     list(iter_file_xrdataset(inwin, ftype="win")),
        >>>     dim="time",
        >>>     coords="different",
        >>>     data_vars="different",
        >>> )
        >>> pres = xr.concat(
        >>>     list(iter_file_xrdataset(inpre, ftype="pre")),
        >>>     dim="time",
        >>>     coords="different",
        >>>     data_vars="different",
        >>> )

    """
    dims = ("time", "yi", "xi")
    for data in owiWinPre.iter_file(fpath, ftype):
        ftype = data[0]
        time = data[1]["time"]
        lat, lon = owiWinPre.getlatlon_fromhead(data[1])
        ds = xr.Dataset()
        ds["lat"] = xr.DataArray(lat.astype(np.float32), dims=("yi", "xi"))
        ds["lon"] = xr.DataArray(lon.astype(np.float32), dims=("yi", "xi"))
        ds["lon"].attrs = {
            "units": "degrees_east",
            "standard_name": "longitude",
            "axis": "X",
            "coordinates": "time lat lon",
        }
        ds["lat"].attrs = {
            "units": "degrees_north",
            "standard_name": "latitude",
            "axis": "Y",
            "coordinates": "time lat lon",
        }
        ds["time"] = xr.DataArray([pd.to_datetime(time)], dims="time")
        ds["time"].attrs = {"axis": "T", "coordinates": "time"}
        if ftype == "win":
            # ds["WS"] = xr.DataArray(
            #    data[2].reshape((1, *ds.lat.shape)), dims=dims, type=np.float32
            # )
            ds["U10"] = xr.DataArray(
                data[3].reshape((1, *ds.lat.shape)).astype(np.float32), dims=dims
            )
            ds["V10"] = xr.DataArray(
                data[4].reshape((1, *ds.lat.shape)).astype(np.float32), dims=dims
            )
            # ds["WS"].attrs = {"units": "m s-1", "coordinates": "time lat lon"}
            ds["U10"].attrs = {"units": "m s-1", "coordinates": "time lat lon"}
            ds["V10"].attrs = {"units": "m s-1", "coordinates": "time lat lon"}
        else:
            ds["PSFC"] = xr.DataArray(
                data[2].reshape((1, *ds.lat.shape)).astype(np.float32), dims=dims
            )
            ds["PSFC"].attrs = {"units": "mb", "coordinates": "time lat lon"}
        yield ds


def add_global_nws13(
    output,
    group_order,
    institution="Oceanweather Inc. (OWI)",
    conventions=["CF-1.6", "OWI-NWS13"],
):
    with netCDF4.Dataset(output, "a") as nc:
        nc.setncattr("group_order", " ".join(group_order))
        nc.setncattr("institution", institution)
        nc.setncattr("conventions", " ".join(conventions))


def winpre2nws13(inwinfiles, inprefiles, ingroups, output, mode="w", verbose=False):
    if mode == "a":
        data = netCDF4.Dataset(output)
        group_order = data.group_order.split()
        rank = len(group_order) + 1
        outgroups = group_order + ingroups
        data.close()
    else:
        rank = 1
        outgroups = ingroups
    for inwin, inpre, group in zip(inwinfiles, inprefiles, ingroups):
        if verbose:
            print(f"Reading {group}, Rank {rank} WIN")
            print(f"Win {inwin}, Pre {inpre}")
        wind = xr.concat(
            list(iter_file_xrdataset(inwin, ftype="win")),
            dim="time",
            coords="different",
            data_vars="different",
        )
        # for i, iwin in enumerate(iter_file_xrdataset(inwin)):
        #     if i > 0:
        #         wind = xr.concat((wind, iwin), dim="time")
        #     else:
        #         wind = iwin
        print("Wind Dataset", wind)
        if verbose:
            print(f"Writing WIN {group}, Rank {rank} to {output} with mode={mode}")
        to_nws13_group(output, group, rank, wind=wind, mode=mode)
        mode = "a"
        del wind
        if verbose:
            print(f"Reading {group}, Rank {rank} PRE")
        pres = xr.concat(
            list(iter_file_xrdataset(inpre, ftype="pre")),
            dim="time",
            coords="different",
            data_vars="different",
        )
        # for i, iwin in enumerate(iter_file_xrdataset(inpre)):
        #     if i > 0:
        #         pres = xr.concat((wind, iwin), dim="time")
        #     else:
        #         pres = iwin
        if verbose:
            print("Pressure Dataset", pres)
            print(f"Writing PRE {group}, Rank {rank} to {output} with mode={mode}")
        to_nws13_group(output, group, rank, pre=pres, mode=mode)
        rank += 1
        del pres
    if verbose:
        print(
            "Writing global attributes, group_order='{0}'".format(" ".join(outgroups))
        )
    add_global_nws13(output, outgroups)


if __name__ == "__main__":
    inputs = parser.parse_args()
    if inputs.verbose:
        print(inputs)
    inwinfiles = inputs.win
    inprefiles = inputs.pre
    ingroups = inputs.group
    output = inputs.output
    if inputs.append:
        mode = "a"
    else:
        mode = "w"
    # TODO check that the file inputs exist
    winpre2nws13(
        inwinfiles, inprefiles, ingroups, output, mode=mode, verbose=inputs.verbose
    )