miχpods for LES simulations#

import ast
import os

import dcpy
import intake
import matplotlib as mpl
import matplotlib.pyplot as plt
import xarray as xr
import xgcm
from datatree import DataTree
import hvplot.xarray

import pump
from pump import mixpods

plt.rcParams["figure.dpi"] = 140

import holoviews as hv

hv.notebook_extension("bokeh")

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:20: FutureWarning: tmpfile is deprecated and will be removed in a future release. Please use dask.utils.tmpfile instead.
  from distributed.utils import tmpfile

Read data#

LES#

les_catalog = intake.open_esm_datastore(
    "../catalogs/pump-les-catalog.json",
    read_csv_kwargs={"converters": {"variables": ast.literal_eval}},
)
les_catalog.df
length kind longitude latitude month path variables
0 5-day average -120 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
1 5-day average -165 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
2 5-day average -115 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
3 5-day average -110 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
4 5-day average -105 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
5 5-day average -100 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
6 5-day average -160 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
7 5-day average -155 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
8 5-day average -150 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
9 5-day average -145 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
10 5-day average -140 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
11 5-day average -135 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
12 5-day average -130 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
13 5-day average -125 0.0 may /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
14 month average -140 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
15 month mooring -140 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [u, v, w, temp, salt, nu_sgs, kappa_sgs, alpha...
16 month average -140 1.5 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
17 month mooring -140 1.5 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [u, v, w, temp, salt, nu_sgs, kappa_sgs, alpha...
18 month average -140 -1.5 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
19 month mooring -140 -1.5 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [u, v, w, temp, salt, nu_sgs, kappa_sgs, alpha...
20 5-day average -120 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
21 5-day average -165 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
22 5-day average -115 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
23 5-day average -110 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
24 5-day average -105 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
25 5-day average -100 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
26 5-day average -160 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
27 5-day average -155 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
28 5-day average -150 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
29 5-day average -145 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
30 5-day average -135 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
31 5-day average -130 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
32 5-day average -125 0.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
33 month average -140 3.0 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
34 month average -140 4.5 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [ume, vme, tempme, saltme, urms, vrms, wrms, t...
35 month mooring -140 4.5 oct /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_r... [u, v, w, temp, salt, nu_sgs, kappa_sgs, alpha...

Only month-long “moorings”

mooring_datasets = les_catalog.search(kind="mooring", length="month").to_dataset_dict(
    preprocess=pump.les.preprocess_les_dataset
)
moorings = DataTree.from_dict(mooring_datasets).squeeze()

--> The keys in the returned dictionary of datasets are constructed as follows:
	'latitude.longitude.month.kind.length'
100.00% [4/4 00:00<00:00] 
avg_datasets = catalog.search(kind="average", length="month").to_dataset_dict(
    preprocess=pump.les.preprocess_les_dataset
)
avgs = DataTree.from_dict(avg_datasets).squeeze()
avgs = avgs.rename_vars({"ume": "u", "vme": "v"})les_

--> The keys in the returned dictionary of datasets are constructed as follows:
	'latitude.longitude.month.kind.length'
100.00% [5/5 00:10<00:00] 
def clear(dataset):
    new = dataset.copy()
    for var in dataset.variables:
        del new[var]
    return new


def clear_root(tree):
    new = tree.copy()
    for var in tree.ds.variables:
        del new.ds[var]
    return new


def extract(tree, varnames):
    return tree.map_over_subtree(lambda ds: ds[varnames])


def to_dataset(tree, dim):
    return xr.concat(
        [
            child.ds.expand_dims({dim: [name]} if dim not in child.ds else dim)
            for name, child in tree.children.items()
        ],
        dim=dim,
    )


def add_ancillary_mixpod_variables(tree):
    from datatree import DataTree

    tree = clear_root(tree)

    # grid = xgcm.Grid(
    #    avgs["0.0.-140.oct.average.month"].ds,
    #    coords={"Z": {"center": "z", "inner": "zc"}},
    #    metrics={("Z",): "dz"},
    # )

    tree.map_over_subtree_inplace(mixpods.prepare)
    tree = tree.assign({"n2s2pdf": mixpods.pdf_N2S2})
    tree["n2s2pdf"] = (
        to_dataset(extract(tree, "n2s2pdf"), dim="latitude")
        .sortby("latitude")
        .to_array()
        .squeeze("variable")
        .load()
    )
    return tree

moorings = add_ancillary_mixpod_variables(moorings)
avgs = add_ancillary_mixpod_variables(avgs)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))

Read TAO#

tao_gridded = xr.open_dataset(
    os.path.expanduser("~/work/pump/zarrs/tao-gridded-ancillary.zarr"),
    chunks="auto",
    engine="zarr",
).sel(longitude=-140, time=slice("2005-Jun", "2015"))
tao_gridded["depth"].attrs["axis"] = "Z"
# eucmax exists
tao_gridded.coords["eucmax"] = pump.calc.get_euc_max(
    tao_gridded.u.reset_coords(drop=True), kind="data"
)
# pump.calc.calc_reduced_shear(tao_gridded)
tao_gridded.coords["enso_transition"] = pump.obs.make_enso_transition_mask().reindex(
    time=tao_gridded.time, method="nearest"
)

tao_gridded = tao_gridded.update(
    {
        "n2s2pdf": mixpods.pdf_N2S2(
            tao_gridded[["S2", "N2T"]]
            .drop_vars(["shallowest", "zeuc"])
            .rename_vars({"N2T": "N2"})
        ).load()
    }
)

ds = tao_gridded[["N2", "S2"]].drop_vars(["shallowest", "zeuc"])

import flox.xarray
import numpy as np

tao_gridded["n2s2pdf_monthly"] = mixpods.to_density(
    flox.xarray.xarray_reduce(
        ds.S2,
        np.log10(4 * ds.N2),
        np.log10(ds.S2),
        ds.time.dt.month,
        func="count",
        expected_groups=(np.linspace(-5, -2, 30), np.linspace(-5, -2, 30), None),
        isbin=(True, True, False),
    ).load()
)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:248: UserWarning: The specified Dask chunks separate the stored chunks along dimension "depth" starting at index 58. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:248: UserWarning: The specified Dask chunks separate the stored chunks along dimension "time" starting at index 139586. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:248: UserWarning: The specified Dask chunks separate the stored chunks along dimension "longitude" starting at index 2. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))

fg = tao_gridded.n2s2pdf_monthly.plot(col="month", col_wrap=4, robust=True)
fg.map(dcpy.plots.line45)

<xarray.plot.facetgrid.FacetGrid at 0x2b2714aeaa40>
_images/19d05e0ed91ca353df8ce0bc4e4cd8b699195bdcf21d3ba1a9f35dc912afa49b.png

Read microstructure#

ls ~/work/datasets/microstructure/osu/

adcp_eq08_30sec.mat        T_0_10W_monthly.mat
adcp_eq08.mat              T_0_110W_monthly.mat
chipod/                    T_0_140W_monthly.mat
chipods_0_10W_hourly.mat   T_0_23W_monthly.mat
chipods_0_110W_hourly.mat  tao0N140W_make_summary_all_deployments.m
chipods_0_110W.nc          th84_timeseries_2009.mat
chipods_0_140W_hourly.mat  tiwe.nc
chipods_0_140W.nc          tropicheat.nc
chipods_0_23W_hourly.mat   tw91_2009.mat
eq08_EUC.mat               tw91_sum.mat
eq08_sum_deglitched.mat    tw91_velocity.mat
equix/                     vel_0_10W_monthly.mat
equix.nc                   vel_0_110W_monthly.mat
mfiles/                    vel_0_140W_monthly.mat
notebooks/                 vel_0_23W_monthly.mat
osu_apl_eps.png

micro = mixpods.load_microstructure()
micro

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:771: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)

<xarray.Dataset>
Dimensions:  ()
Data variables:
    *empty*

Stability Diagram#

  1. Contours enclose 50%, 75% of the data

  2. Colors are for LES

  3. White contours are 50%, 75% contours for TAO 0, 140

  4. thnk about LES resampling?

Domain averages vs moorings#

Averaging over the domain (red) tightens the scatter plot along the 1:1 line

f, ax = plt.subplots(1, 4, sharex=True, sharey=True, constrained_layout=True)
for lat, axx in zip(moorings.ds["latitude"].data, ax):
    mixpods.plot_n2s2pdf(
        moorings.ds.n2s2pdf.sel(enso_transition_phase="none", latitude=lat),
        ax=axx,
        pcolor=False,
        colors="k",
    )
    mixpods.plot_n2s2pdf(
        avgs.ds.n2s2pdf.sel(enso_transition_phase="none", latitude=lat),
        ax=axx,
        pcolor=False,
        colors="r",
    )
    axx.set_title(f"latitude={lat}")
dcpy.plots.clean_axes(ax)
f.set_size_inches((8, 6))
_images/45f4acfe223f62cecdb7f396ca3ae4aafe6b5051c8aa50ab802fe0584dad250b.png

vs TAO#

f, ax = plt.subplots(1, 4, sharex=True, sharey=True, constrained_layout=True)
for lat, axx in zip(moorings.ds["latitude"].data, ax):
    mixpods.plot_n2s2pdf(
        moorings.ds.n2s2pdf.sel(latitude=lat),
        ax=axx,
        vmin=0,
        vmax=0.8,
        add_colorbar=False,
        cmap="viridis",
    )

    mixpods.plot_n2s2pdf(
        tao_gridded.n2s2pdf.sel(enso_transition_phase="none"),
        ax=axx,
        pcolor=False,
        colors="w",
    )
    axx.set_title(f"latitude={lat}")
dcpy.plots.clean_axes(ax)
f.set_size_inches((8, 6))
_images/2ff518fcb63b7a0c2afa47ffbbb0c1daa133540d4e9fe8037ca5ccf252b61262.png

vs TAO in October#

Comparing October contours only in TAO makes it look better

f, ax = plt.subplots(1, 4, sharex=True, sharey=True, constrained_layout=True)
for lat, axx in zip(moorings.ds.latitude.data, ax):
    mixpods.plot_n2s2pdf(
        moorings.ds.n2s2pdf.sel(latitude=lat),
        ax=axx,
        vmin=0,
        vmax=0.8,
        add_colorbar=False,
        cmap="viridis",
    )

    mixpods.plot_n2s2pdf(
        tao_gridded.n2s2pdf_monthly.sel(month=10), ax=axx, pcolor=False, colors="w"
    )
    axx.set_title(f"latitude={lat}")
dcpy.plots.clean_axes(ax)
f.set_size_inches((8, 6))
_images/3b38794325dc3a3fe223d8a389c971bfeeddb49011a0ccb2a15af62367144067.png

Here is TAO 0, 140 by itself

mixpods.plot_n2s2pdf(
    tao_gridded.n2s2pdf.sel(enso_transition_phase="none"),
)

<matplotlib.contour.QuadContourSet at 0x2af9d92f38b0>
_images/4184133bb68a088e2e1f068873e89b9928a20a57c1d1e4cd660f2731617d377d.png

vs microstructure#

def integrate(da, dim):

    import pandas as pd

    bin_ = da[dim][0].data.item()
    result = da.sum(dim) * bin_.length
    for other_dim in result.dims:
        maybe_bin = da[other_dim][0].data.item()
        if isinstance(maybe_bin, pd.Interval):
            result[other_dim] = pd.IntervalIndex(result[other_dim].data).mid.to_numpy()
    return result
Hide code cell source 
to_plot = {
    "/les moor": moorings["0.0.-140.oct.mooring.month"].ds.load(),
    "/les avg": avgs["0.0.-140.oct.average.month"].ds.load(),
    "/tao oct": tao_gridded[["n2s2pdf_monthly"]]
    .sel(month=10)
    .rename({"n2s2pdf_monthly": "n2s2pdf"}),
}
to_plot.update(micro.to_dict())

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in true_divide
  return func(*(_execute_task(a, cache) for a in args))
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
Hide code cell source 
colors = {
    "les avg": "magenta",
    "les moor": "red",
    "tao oct": "black",
    "equix": "green",
    "tropicheat": "darkgray",
    "tiwe": "blue",
}

hdls = []
s2hdls = []
n2hdls = []
for name, node in to_plot.items():
    if node.data_vars:
        kwargs = dict(line_width=2, label=name[1:])
        hdls.append(
            mixpods.hvplot_n2s2pdf(
                node["n2s2pdf"],
                targets=(0.5,),
                pcolor=False,
                cmap=[colors[name[1:]]],
                **kwargs
            )
        )

        s2hdls.append(
            integrate(node.n2s2pdf.squeeze(), "N2_bins").hvplot.line(
                **kwargs, color=[colors[name[1:]]]
            )
        )
        n2hdls.append(
            integrate(node.n2s2pdf.squeeze(), "S2_bins").hvplot.line(
                **kwargs, color=[colors[name[1:]]]
            )
        )

from functools import reduce
import operator

joint = (
    reduce(operator.mul, hdls)
    # .opts(frame_width=400, frame_height=400)
    # .opts(xlabel="log S2", ylabel="log 4N2")
)
s2 = reduce(operator.mul, s2hdls)
n2 = reduce(operator.mul, n2hdls)

(
    joint
    << n2.opts(show_legend=False)  # frame_width=200, frame_height=400)
    << s2.opts(show_legend=False)  ##.opts(frame_width=400, frame_height=200)
).opts(width=500)

(n2.opts(show_legend=True) + s2.opts(show_legend=True)).cols(1)

TODO#

  1. Combined 5-day distributions: all longitudes

  2. hour average for microstructure; maybe vertical scales too

  3. LES Ri PDF with different subsampling z; match with Smyth & Moum (2013)

  4. Add domain depth to catalog

TODO: Turbulence patterns#

moorings = moorings.assign(
    {
        "epsilon": lambda ds: (15 / 4 * ds.nu_sgs * ds.S2)
        if "nu_sgs" in ds
        else xr.DataArray([0])
    }
)
del moorings.ds["epsilon"]
moorings

<xarray.Dataset>
Dimensions:                (latitude: 4, enso_transition_phase: 1, N2_bins: 29,
                            S2_bins: 29)
Coordinates:
    dz                     float64 0.5
    longitude              int64 140
  * latitude               (latitude) float64 -1.5 0.0 1.5 4.5
  * enso_transition_phase  (enso_transition_phase) object 'none'
  * N2_bins                (N2_bins) object (-5.0, -4.896551724137931] ... (-...
  * S2_bins                (S2_bins) object (-5.0, -4.896551724137931] ... (-...
    bin_areas              (N2_bins, S2_bins) float64 0.0107 0.0107 ... 0.0107
    variable               <U7 'n2s2pdf'
Data variables:
    n2s2pdf                (latitude, enso_transition_phase, N2_bins, S2_bins) float64 ...
moorings.ds["epsilon"] = (
    to_dataset(extract(clear_root(moorings), "epsilon"), dim="latitude")
    .to_array()
    .squeeze()
)

dz
longitude
latitude
enso_transition_phase
N2_bins
S2_bins
bin_areas
variable
n2s2pdf

avgs.ds["epsilon"] = (
    to_dataset(extract(clear_root(avgs), "epsilon"), dim="latitude")
    .to_array()
    .squeeze()
)

dz
longitude
latitude
enso_transition_phase
N2_bins
S2_bins
bin_areas
variable
n2s2pdf

Build Catalog#

Using ecgtools

Dan:

(100 W is grid point “1402” and 165 W is grid point “102”; 140W is “602” first date “25MAY1985” indicates first simulated date; second date 24NOV2021 is the nominal date I started the runs):

import glob
import pathlib

import ecgtools
import intake
import intake_esm
import xarray as xr

root = (
    "/glade/p/cgd/oce/people/dwhitt/TPOS/"
    "tpos_LES_runs_setup_scripts/tpos-DIABLO/diablo_2.0/post_process/diablo_analysis"
)

def parse_les_file(path):
    path = pathlib.Path(path)
    split = path.stem.split("_")

    # print(split)
    # print(path)

    if split[3] == "242":
        length = "5-day"
        kind = "average"
        longitude = int(-165 + (int(split[4]) - 102) * 0.05)
        latitude = 0
    else:
        length = "month"
        kind = "mooring" if "mooring" in split[-1] else "average"
        longitude = -1 * int(split[3][-4:-1])
        latitude = float(split[3][:-5]) * (-1 if "S" in split[3] else 1)

    with xr.open_dataset(path) as ds:
        info = {
            "length": length,
            "kind": kind,
            "longitude": longitude,
            "latitude": latitude,
            "month": "may" if "MAY" in split[-2] else "oct",
            "path": path,
            "variables": [k for k in ds],
        }
    return info


parse_les_file(f"{root}/ROMS_PSH_6HRLIN_0N140W_360x360x216_22OCT2020.nc")

{'length': 'month',
 'kind': 'average',
 'longitude': -140,
 'latitude': 0.0,
 'month': 'oct',
 'path': PosixPath('/glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_runs_setup_scripts/tpos-DIABLO/diablo_2.0/post_process/diablo_analysis/ROMS_PSH_6HRLIN_0N140W_360x360x216_22OCT2020.nc'),
 'variables': ['ume',
  'vme',
  'tempme',
  'saltme',
  'urms',
  'vrms',
  'wrms',
  'temprms',
  'saltrms',
  'uw',
  'vw',
  'saltw',
  'tempw',
  'nududz',
  'nudvdz',
  'kappadsdz',
  'kappadtdz',
  'dTdtSOLAR',
  'dTdtRESTORE',
  'dTdtFORCE',
  'dUdtRESTORE',
  'dUdtFORCE',
  'dVdtRESTORE',
  'dVdtFORCE',
  'epsilon',
  'N2',
  'S2',
  'RIG',
  'nududztop',
  'nudvdztop',
  'kappadsdztop',
  'kappadtdztop',
  'nududzbot',
  'nudvdzbot',
  'kappadsdzbot',
  'kappadtdzbot',
  'alpha',
  'beta',
  'T0',
  'S0',
  'rho0']}

builder = ecgtools.Builder(
    root,
    exclude_patterns=[
        "*test*",
        "*irene_*",
        "*spectra*",
        "*fixedeps*",
        "*timeavg*",
        "*0N140W*5OCT202*",
        "*0N140W*20OCT202*",
        "*0N140W*29OCT202*",
        "*6mavg*",
        "*54x54*",
    ],
    njobs=-1,
)
builder.build(parse_les_file)

[Parallel(n_jobs=-1)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=-1)]: Done   1 out of   1 | elapsed:    0.9s finished
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=-1)]: Done  14 tasks      | elapsed:    3.2s
[Parallel(n_jobs=-1)]: Done  36 out of  36 | elapsed:   15.7s finished

Builder(root_path=PosixPath('/glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_runs_setup_scripts/tpos-DIABLO/diablo_2.0/post_process/diablo_analysis'), extension='.nc', depth=0, exclude_patterns=['*test*', '*irene_*', '*spectra*', '*fixedeps*', '*timeavg*', '*0N140W*5OCT202*', '*0N140W*20OCT202*', '*0N140W*29OCT202*', '*6mavg*', '*54x54*'], njobs=-1)

Summarize catalog#

builder.df.drop("variables", axis=1).groupby(
    ["length", "kind", "latitude", "longitude", "month"]
).count()
path
length kind latitude longitude month
5-day average 0.0 -165 may 1
oct 1
-160 may 1
oct 1
-155 may 1
oct 1
-150 may 1
oct 1
-145 may 1
oct 1
-140 may 1
-135 may 1
oct 1
-130 may 1
oct 1
-125 may 1
oct 1
-120 may 1
oct 1
-115 may 1
oct 1
-110 may 1
oct 1
-105 may 1
oct 1
-100 may 1
oct 1
month average -1.5 -140 oct 1
0.0 -140 oct 1
1.5 -140 oct 1
3.0 -140 oct 1
4.5 -140 oct 1
mooring -1.5 -140 oct 1
0.0 -140 oct 1
1.5 -140 oct 1
4.5 -140 oct 1 
builder.save(
    "../catalogs/pump-les-catalog.csv",
    path_column_name="path",
    variable_column_name="variables",
    data_format="netcdf",
    groupby_attrs=["latitude", "longitude", "month", "kind", "length"],
    aggregations=[
        {"type": "union", "attribute_name": "variables"},
        {"type": "join_new", "attribute_name": "latitude"},
        {"type": "join_new", "attribute_name": "longitude"},
    ],
)

Saved catalog location: ../catalogs/pump-les-catalog.json and ../catalogs/pump-les-catalog.csv

import ast
from datatree import DataTree

catalog = intake.open_esm_datastore(
    "../catalogs/pump-les-catalog.json",
    csv_kwargs={"converters": {"variables": ast.literal_eval}},
)
catalog.df

dataset_dict = catalog.to_dataset_dict(cdf_kwargs={})

--> The keys in the returned dictionary of datasets are constructed as follows:
	'latitude.longitude.month.kind.length'
100.00% [36/36 00:28<00:00] 
tree = DataTree.from_dict(dataset_dict)

tree["-140.may"].ds

<xarray.Dataset>
Dimensions:       (z: 288, time: 1313, longitude: 1, month: 1)
Coordinates:
  * z             (z) float64 -143.5 -143.0 -142.5 ... -1.0 -0.5 9.969e+36
  * time          (time) datetime64[ns] 1985-05-25T03:00:00 ... 1985-05-30T00...
  * month         (month) <U3 'may'
  * longitude     (longitude) int64 -140
Data variables: (12/41)
    ume           (longitude, month, time, z) float32 dask.array<chunksize=(1, 1, 1313, 288), meta=np.ndarray>
    vme           (longitude, month, time, z) float32 dask.array<chunksize=(1, 1, 1313, 288), meta=np.ndarray>
    tempme        (longitude, month, time, z) float32 dask.array<chunksize=(1, 1, 1313, 288), meta=np.ndarray>
    saltme        (longitude, month, time, z) float32 dask.array<chunksize=(1, 1, 1313, 288), meta=np.ndarray>
    urms          (longitude, month, time, z) float32 dask.array<chunksize=(1, 1, 1313, 288), meta=np.ndarray>
    vrms          (longitude, month, time, z) float32 dask.array<chunksize=(1, 1, 1313, 288), meta=np.ndarray>
    ...            ...
    kappadtdzbot  (longitude, month, time) float32 dask.array<chunksize=(1, 1, 1313), meta=np.ndarray>
    alpha         (longitude, month) float64 0.0002976
    beta          (longitude, month) float64 -0.0007386
    T0            (longitude, month) float64 25.0
    S0            (longitude, month) float64 35.25
    rho0          (longitude, month) float64 1.024e+03
Attributes:
    type:                    DIABLO LES, processed means
    title:                   ROMS_PSH_6HRLIN_0N140W_360x360x288_5OCT2021
    history:                 Tue Dec 28 09:09:16 2021: ncatted -O -a units,ti...
    NCO:                     netCDF Operators version 4.9.5 (Homepage = http:...
    intake_esm_varname:      ('ume', 'vme', 'tempme', 'saltme', 'urms', 'vrms...
    intake_esm_dataset_key:  -140.may