Skip to content

Instantly share code, notes, and snippets.

@acrosby

acrosby/generic_nws13_example.py

Created January 18, 2023 22:18
Show Gist options
  • Save acrosby/b2acea4488ee437cc9564d9b7054a228 to your computer and use it in GitHub Desktop.
Save acrosby/b2acea4488ee437cc9564d9b7054a228 to your computer and use it in GitHub Desktop.
Download ZIP
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
Raw
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)
Raw
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
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment