Bulk Ri analysis#

Consolidating code for bulk Ri paper

%load_ext watermark

import dcpy
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import facetgrid
import pump

mpl.rcParams["savefig.dpi"] = 300
mpl.rcParams["savefig.bbox"] = "tight"
mpl.rcParams["figure.dpi"] = 140

xr.set_options(keep_attrs=True)

%watermark -iv

pump      : 0.1
xarray    : 0.17.1.dev3+g48378c4b1
matplotlib: 3.3.4
dcpy      : 0.1
numpy     : 1.20.1

import dask
import distributed
import ncar_jobqueue

cluster = ncar_jobqueue.NCARCluster(project="ncgd0011")
cluster.scale(4)
client = distributed.Client(cluster)
client

/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/ncar_jobqueue/cluster.py:29: UserWarning: Unable to determine which NCAR cluster you are running on...Using a local cluster via `distributed.LocalCluster`.
  warn(

Client

Cluster

  • Workers: 0
  • Cores: 0
  • Memory: 0 B

Read datasets#

mimoc = dcpy.oceans.read_mimoc("/home/deepak/datasets/mimoc/")
mimoc["longitude"] = mimoc.longitude - 360
mimocmld = (
    mimoc.DEPTH_MIXED_LAYER.coarsen(longitude=5, boundary="pad")
    .mean()
    .sel(longitude=slice(-220, -80))
    .sel(latitude=0, method="nearest")
    .max("time")
)

mimocmld.plot()

/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide
  x = np.divide(x1, x2, out)

[<matplotlib.lines.Line2D at 0x7f9ffe6faa60>]
_images/4c2afa56774bc95ca458633566c8960ba3f998ab4fdcd7e15c1418578bea2c16.png 
johnson = pump.obs.read_johnson("~/datasets/johnson.cdf").sel(latitude=0)
johnson["eucmax"] = pump.get_euc_max(johnson.u)
johnson = pump.calc.estimate_euc_depth_terms(johnson)
johnson.attrs["name"] = "Johnson"
johnson["mld"] = mimocmld.interp(longitude=johnson.longitude.values).load()
johnson.load()

/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide
  x = np.divide(x1, x2, out)

<xarray.Dataset>
Dimensions:    (depth: 50, lon_edges: 11, longitude: 10)
Coordinates:
  * longitude  (longitude) float64 -217.0 -204.0 -195.0 ... -125.0 -110.0 -95.0
  * lon_edges  (lon_edges) float64 136.5 149.5 160.5 172.5 ... 242.5 257.5 272.5
    latitude   float64 0.0
  * depth      (depth) float64 -5.0 -15.0 -25.0 -35.0 ... -475.0 -485.0 -495.0
Data variables: (12/18)
    temp       (depth, longitude) float32 29.32 29.32 28.99 ... 8.132 8.024
    salt       (depth, longitude) float32 34.19 34.45 34.77 ... 34.63 34.63
    dens       (depth, longitude) float32 1.021e+03 1.022e+03 ... 1.027e+03
    u          (depth, longitude) float32 0.08696 0.1127 0.01524 ... nan nan nan
    TSPTS      (depth, longitude) float32 20.0 34.0 28.0 46.0 ... 47.0 38.0 24.0
    UPTS       (depth, longitude) float32 20.0 34.0 27.0 45.0 ... nan 2.362 nan
    ...         ...
    b          (depth, longitude) float32 -9.766 -9.774 -9.776 ... -9.822 -9.822
    bs         (longitude) float32 -9.767 -9.778 -9.777 ... -9.794 -9.791 -9.79
    beuc       (longitude) float32 -9.81 -9.81 -9.81 -9.81 ... -9.81 -9.81 -9.81
    db         (longitude) float32 0.04333 0.03237 0.03319 ... 0.01902 0.01996
    Rib        (longitude) float64 56.18 20.04 36.27 12.32 ... 1.185 1.453 1.692
    mld        (longitude) float64 29.57 45.73 38.86 37.31 ... 25.44 17.44 17.51
Attributes:
    history:  FERRET V5.41    1-Oct-02
    name:     Johnson
tao = pump.obs.read_tao_zarr("ancillary").rename({"Ri": "Rig"})
tao

<xarray.Dataset>
Dimensions:             (depth: 61, longitude: 5, time: 287335)
Coordinates:
    deepest             (time, longitude) float64 dask.array<chunksize=(100000, 1), meta=np.ndarray>
  * depth               (depth) float64 -300.0 -295.0 -290.0 ... -10.0 -5.0 0.0
    eucmax              (time, longitude) float64 dask.array<chunksize=(287335, 4), meta=np.ndarray>
    latitude            float32 ...
  * longitude           (longitude) float64 -204.0 -195.0 -170.0 -140.0 -110.0
    mld                 (time, longitude) float64 dask.array<chunksize=(100000, 1), meta=np.ndarray>
    reference_pressure  int64 ...
    shallowest          (time, longitude) float64 dask.array<chunksize=(100000, 1), meta=np.ndarray>
  * time                (time) datetime64[ns] 1988-05-15T18:00:00 ... 2021-02-24
    zeuc                (depth, time, longitude) float64 dask.array<chunksize=(61, 287335, 4), meta=np.ndarray>
Data variables:
    N2                  (time, longitude, depth) float64 dask.array<chunksize=(100000, 1, 61), meta=np.ndarray>
    N2T                 (time, longitude, depth) float64 dask.array<chunksize=(100000, 1, 61), meta=np.ndarray>
    Rig                 (time, longitude, depth) float64 dask.array<chunksize=(100000, 1, 61), meta=np.ndarray>
    Rig_T               (time, longitude, depth) float64 dask.array<chunksize=(100000, 1, 61), meta=np.ndarray>
    S                   (time, longitude, depth) float64 dask.array<chunksize=(100000, 1, 61), meta=np.ndarray>
    S2                  (time, depth, longitude) float32 dask.array<chunksize=(100000, 61, 1), meta=np.ndarray>
    T                   (time, longitude, depth) float64 dask.array<chunksize=(100000, 1, 61), meta=np.ndarray>
    dens                (time, longitude, depth) float64 dask.array<chunksize=(100000, 1, 61), meta=np.ndarray>
    densT               (time, longitude, depth) float64 dask.array<chunksize=(100000, 1, 61), meta=np.ndarray>
    u                   (time, depth, longitude) float32 dask.array<chunksize=(100000, 61, 1), meta=np.ndarray>
    v                   (time, depth, longitude) float32 dask.array<chunksize=(100000, 61, 1), meta=np.ndarray>
median_Rig_z = (
    tao[["Rig", "Rig_T"]]
    .groupby("time.season")
    .apply(lambda x: x.chunk({"time": -1}).quantile(q=0.5, dim="time"))
    .reindex(season=["DJF", "MAM", "JJA", "SON"])
)
median_Rig_z.load()

/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/numpy/lib/nanfunctions.py:1389: RuntimeWarning: All-NaN slice encountered
  result = np.apply_along_axis(_nanquantile_1d, axis, a, q,

<xarray.Dataset>
Dimensions:    (depth: 61, longitude: 5, season: 4)
Coordinates:
  * season     (season) <U3 'DJF' 'MAM' 'JJA' 'SON'
  * depth      (depth) float64 -300.0 -295.0 -290.0 -285.0 ... -10.0 -5.0 0.0
  * longitude  (longitude) float64 -204.0 -195.0 -170.0 -140.0 -110.0
    quantile   float64 0.5
Data variables:
    Rig        (season, longitude, depth) float64 nan nan nan ... nan nan nan
    Rig_T      (season, longitude, depth) float64 nan nan 0.582 ... nan nan nan
tao_seasonal = (
    tao.groupby("time.season").mean().reindex(season=["DJF", "MAM", "JJA", "SON"])
)
tao_seasonal["eucmax"] = pump.calc.get_euc_max(tao_seasonal.u)
tao_seasonal = pump.calc.estimate_euc_depth_terms(tao_seasonal)
tao_seasonal.attrs["name"] = "TAO"
tao_seasonal.load()

# is this useful?
# tao_clim["Rig"] = pump.calc.calc_tao_ri(
#    tao_clim[["u", "v"]].expand_dims(time=1),
#    tao_clim["T"].expand_dims(time=1)
# )

/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide
  x = np.divide(x1, x2, out)
/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide
  x = np.divide(x1, x2, out)

<xarray.Dataset>
Dimensions:    (depth: 101, longitude: 5, season: 4)
Coordinates:
  * season     (season) object 'DJF' 'MAM' 'JJA' 'SON'
  * depth      (depth) float64 -500.0 -495.0 -490.0 -485.0 ... -10.0 -5.0 0.0
    latitude   float32 0.0
  * longitude  (longitude) float64 -204.0 -195.0 -170.0 -140.0 -110.0
Data variables: (12/16)
    Rig        (season, longitude, depth) float64 nan nan nan ... nan nan nan
    T          (season, longitude, depth) float64 8.184 8.241 ... 23.07 23.04
    dens       (season, longitude, depth) float64 1.025e+03 ... 1.024e+03
    u          (season, depth, longitude) float32 nan nan nan ... nan nan
    v          (season, depth, longitude) float32 nan nan nan ... nan nan
    eucmax     (season, longitude) float64 -190.0 -190.0 -170.0 ... -115.0 -85.0
    ...         ...
    dens_euc   (season, longitude) float64 1.024e+03 1.024e+03 ... 1.025e+03
    b          (season, longitude, depth) float64 -9.815 -9.816 ... -9.799 -9.8
    bs         (season, longitude) float64 -9.784 -9.785 -9.794 ... -9.801 -9.8
    beuc       (season, longitude) float64 -9.81 -9.81 -9.81 ... -9.81 -9.81
    db         (season, longitude) float64 0.02576 0.02522 ... 0.009429 0.01011
    Rib        (season, longitude) float64 31.47 10.08 2.628 ... 0.7851 1.134
Attributes:
    name:     TAO
tao_clim = tao.mean("time")
tao_clim["eucmax"] = pump.calc.get_euc_max(tao_clim.u)
tao_clim = pump.calc.estimate_euc_depth_terms(tao_clim)
tao_clim.attrs["name"] = "TAO"
tao_clim.load()

# is this useful?
tao_clim["Rig"] = pump.calc.calc_tao_ri(
    tao_clim[["u", "v"]].expand_dims(time=1), tao_clim["T"].expand_dims(time=1)
)

/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide
  x = np.divide(x1, x2, out)
/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide
  x = np.divide(x1, x2, out)

tao_zeuc = xr.open_zarr("tao-zeuc.zarr", consolidated=True)
tao_zeuc

<xarray.Dataset>
Dimensions:     (longitude: 5, time: 310315, zeuc: 59)
Coordinates:
    deepest     (time, longitude) float64 dask.array<chunksize=(5000, 2), meta=np.ndarray>
    eucmax      (time, longitude) float64 dask.array<chunksize=(5000, 4), meta=np.ndarray>
    latitude    float32 ...
  * longitude   (longitude) float64 -204.0 -195.0 -170.0 -140.0 -110.0
    mld         (time, longitude) float64 dask.array<chunksize=(5000, 2), meta=np.ndarray>
    shallowest  (time, longitude) float64 dask.array<chunksize=(5000, 2), meta=np.ndarray>
  * time        (time) datetime64[ns] 1985-10-01T06:00:00 ... 2021-02-24
  * zeuc        (zeuc) float64 -47.5 -42.5 -37.5 -32.5 ... 232.5 237.5 242.5
Data variables:
    Rig         (time, longitude, zeuc) float64 dask.array<chunksize=(5000, 2, 59), meta=np.ndarray>
    T           (time, longitude, zeuc) float64 dask.array<chunksize=(5000, 2, 59), meta=np.ndarray>
    dens        (time, longitude, zeuc) float64 dask.array<chunksize=(5000, 2, 59), meta=np.ndarray>
    u           (time, longitude, zeuc) float64 dask.array<chunksize=(5000, 2, 59), meta=np.ndarray>
    v           (time, longitude, zeuc) float64 dask.array<chunksize=(5000, 2, 59), meta=np.ndarray>
Ri_q = xr.load_dataarray("tao-hourly-Ri-seasonal-percentiles.nc")
Ri_q.attrs["long_name"] = "$Ri_g$"
Ri_q

<xarray.DataArray 'Rig' (season: 4, longitude: 5, zeuc: 59, quantile: 3)>
array([[[[8.65716985e-01, 1.84699736e+00, 4.07929186e+00],
         [9.68607747e-01, 2.07586669e+00, 4.69831821e+00],
         [1.12714689e+00, 2.45510051e+00, 5.63565819e+00],
         ...,
         [8.25933339e-01, 1.21812980e+00, 2.22266901e+00],
         [3.31773526e-02, 2.15271859e-01, 6.91555336e-01],
         [6.30251313e-01, 9.59637420e-01, 1.28902353e+00]],

        [[7.71597015e-01, 1.43807574e+00, 3.14489151e+00],
         [8.06376970e-01, 1.53085922e+00, 3.27291715e+00],
         [8.95512070e-01, 1.70564919e+00, 3.66061253e+00],
         ...,
         [           nan,            nan,            nan],
         [           nan,            nan,            nan],
         [           nan,            nan,            nan]],

        [[1.07121354e+00, 1.87813410e+00, 3.69584450e+00],
         [1.11844147e+00, 2.00944921e+00, 4.03503153e+00],
         [1.19428670e+00, 2.16651316e+00, 4.50967620e+00],
         ...,
...
         ...,
         [           nan,            nan,            nan],
         [           nan,            nan,            nan],
         [           nan,            nan,            nan]],

        [[6.15421945e-01, 1.04513786e+00, 2.04579183e+00],
         [6.42674890e-01, 1.10834780e+00, 2.11288515e+00],
         [6.63607716e-01, 1.15214874e+00, 2.20694345e+00],
         ...,
         [           nan,            nan,            nan],
         [           nan,            nan,            nan],
         [           nan,            nan,            nan]],

        [[4.20117100e-01, 7.23407176e-01, 1.42003331e+00],
         [4.19239467e-01, 7.26062245e-01, 1.43285230e+00],
         [4.16259694e-01, 7.33555903e-01, 1.46662031e+00],
         ...,
         [           nan,            nan,            nan],
         [           nan,            nan,            nan],
         [           nan,            nan,            nan]]]])
Coordinates:
  * season     (season) object 'DJF' 'MAM' 'JJA' 'SON'
  * longitude  (longitude) float64 -204.0 -195.0 -170.0 -140.0 -110.0
  * zeuc       (zeuc) float64 -47.5 -42.5 -37.5 -32.5 ... 232.5 237.5 242.5
  * quantile   (quantile) float64 0.25 0.5 0.75
    latitude   float32 0.0
    num_obs    (season, longitude, zeuc) int64 2956 3450 3685 4162 ... 0 0 0 0
Attributes:
    long_name:  $Ri_g$
median_Rig = Ri_q.sel(quantile=0.5, zeuc=slice(0, None))
median_Rig.coords["z_above_euc"] = median_Rig.zeuc + johnson.eucmax.interp(
    longitude=median_Rig.longitude.values
)
median_Rig.z_above_euc.attrs["long_name"] = "Height above mean EUC max"
median_Rig.z_above_euc.attrs["units"] = "m"
median_Rig

<xarray.DataArray 'Rig' (season: 4, longitude: 5, zeuc: 49)>
array([[[1.21434122e+01, 5.20354724e+00, 3.50501919e+00, 3.30736311e+00,
         3.33031631e+00, 3.43964060e+00, 2.92148757e+00, 2.48129720e+00,
         2.11243240e+00, 1.83912229e+00, 1.69382899e+00, 1.61765319e+00,
         1.57829252e+00, 1.65197799e+00, 1.67660646e+00, 1.70361004e+00,
         1.75811871e+00, 1.77633452e+00, 1.78274006e+00, 1.77130445e+00,
         1.76453557e+00, 1.59982041e+00, 1.53944503e+00, 1.46579473e+00,
         1.38052897e+00, 1.36555182e+00, 1.22801891e+00, 1.11480848e+00,
         9.81393217e-01, 8.09430894e-01, 7.11033315e-01, 6.64250979e-01,
         6.77348794e-01, 6.59767973e-01, 6.68314872e-01, 6.45685835e-01,
         6.25755328e-01, 6.47728194e-01, 6.16354965e-01, 7.43610355e-01,
         7.47971748e-01, 7.51553840e-01, 6.64788050e-01, 6.37038637e-01,
         7.09365153e-01, 1.08899569e+00, 1.21812980e+00, 2.15271859e-01,
         9.59637420e-01],
        [1.19743455e+01, 4.94836731e+00, 3.32542799e+00, 2.83437588e+00,
         2.52501184e+00, 2.21110672e+00, 1.92064009e+00, 1.70944355e+00,
         1.57992060e+00, 1.47963011e+00, 1.50067241e+00, 1.45032254e+00,
         1.51007357e+00, 1.50672529e+00, 1.58891969e+00, 1.53438176e+00,
         1.58523445e+00, 1.45635478e+00, 1.41217789e+00, 1.23123271e+00,
         1.14912456e+00, 9.52980219e-01, 8.41401552e-01, 6.95025110e-01,
         5.87664619e-01, 5.10455407e-01, 4.66768745e-01, 4.36939883e-01,
...
         2.38389300e-01, 2.20526505e-01, 1.79629081e-01, 1.63788227e-01,
         1.88565053e-01, 1.81312113e-01, 3.56486993e-01, 4.76187062e-02,
                    nan,            nan,            nan,            nan,
                    nan,            nan,            nan,            nan,
                    nan,            nan,            nan,            nan,
                    nan,            nan,            nan,            nan,
                    nan],
        [5.60872189e+00, 1.44782782e+00, 6.76429622e-01, 4.91971220e-01,
         4.57825586e-01, 4.37887622e-01, 4.16151139e-01, 3.88965627e-01,
         3.52532334e-01, 3.19898671e-01, 2.89160931e-01, 2.66395005e-01,
         2.49164377e-01, 2.44682953e-01, 2.48068307e-01, 2.36800210e-01,
         2.37940240e-01, 2.26757634e-01, 2.49925464e-01, 2.65475665e-01,
         2.90705739e-01, 3.15365763e-01, 3.18039032e-01, 4.66446360e-01,
         3.98508291e-01, 4.66065807e-01,            nan,            nan,
                    nan,            nan,            nan,            nan,
                    nan,            nan,            nan,            nan,
                    nan,            nan,            nan,            nan,
                    nan,            nan,            nan,            nan,
                    nan,            nan,            nan,            nan,
                    nan]]])
Coordinates:
  * season       (season) object 'DJF' 'MAM' 'JJA' 'SON'
  * longitude    (longitude) float64 -204.0 -195.0 -170.0 -140.0 -110.0
  * zeuc         (zeuc) float64 2.5 7.5 12.5 17.5 ... 227.5 232.5 237.5 242.5
    quantile     float64 0.5
    latitude     float32 0.0
    num_obs      (season, longitude, zeuc) int64 5189 5238 5231 5230 ... 0 0 0 0
    z_above_euc  (zeuc, longitude) float64 -182.5 -182.5 -152.5 ... 127.5 167.5
Attributes:
    long_name:  $Ri_g$

Test linear velocity profile?#

If the velocity profile is linear then bulk Ri is related to gradient Ri only in marginally stable zone.

Looks like the linear velocity profile is a decent assumption east of 195W, though 170W & 155W are just OK.

f, ax = plt.subplots(1, 2, sharey=True)

(johnson.u + (johnson.longitude / 50)).plot.line(
    hue="longitude", y="depth", add_legend=False, ax=ax[0]
)
(johnson.u.differentiate("depth") + (johnson.longitude / 1000)).plot.line(
    hue="longitude", y="depth", add_legend=False, ax=ax[1]
)
plt.gcf().legend(
    list(map(str, johnson.longitude.values)), title="longitude", loc="right"
)
plt.gcf().set_size_inches((8, 4))
plt.tight_layout()
ax[0].set_title("u (displaced)")
ax[1].set_title("$u_z$ (displaced)")

Text(0.5, 1.0, '$u_z$ (displaced)')
_images/65de8f659c8dd4228f5f2f86ced4bba48fa6763f3fe4feca98be5d5748349b08.png 
mld_interp = -1 * mimocmld.interp(longitude=johnson.longitude).load()
mld_interp.loc[{"longitude": -195}] = -70
mld_interp.loc[{"longitude": -180}] = -60


def plot_mld_eucmax(fg):
    import warnings

    for loc, ax in zip(fg.name_dicts.flat, fg.axes.flat):
        mld = mld_interp.sel(loc)
        eucmax = johnson.eucmax.sel(loc)
        if loc is None:
            continue
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            dcpy.plots.liney([mld, eucmax], ax=ax)
            ax.plot(
                fg.data.sel(loc).interp(depth=[mld, eucmax]).data,
                [mld.data, eucmax.data],
                color="k",
                ls="--",
            )


subset = johnson.sel(longitude=slice(-200, -95), depth=slice(-200))
fg = subset.u.plot.line(col="longitude", y="depth", col_wrap=4)
fg.map(lambda: plt.axvline(0))
plot_mld_eucmax(fg)


fg = subset.dens.plot.line(col="longitude", y="depth", col_wrap=4)
plot_mld_eucmax(fg)

/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide
  x = np.divide(x1, x2, out)
_images/23ff4c643ae33f1634a4417d8dce1506932bd6dd12b5424de8926e8ab7368cab.png _images/6ab33f7758ce4264ffb2823ebd3827dca9197fc0730fe9db5886daa04a4b4964.png

median Ri_g profiles in depth & euc space#

with mpl.rc_context(pump.plot.sm13_cycler):
    fg = facetgrid.facetgrid(
        ["depth", "euc"],
        median_Rig_z.longitude,
        func=xr.plot.line,
        plot_kwargs=dict(
            ylim=[-150, 20],
            xlim=[0.1, 3.5],
            xscale="log",
            _labels=False,
            hue="season",
            add_legend=False,
        ),
    )

fg.map_row("depth", median_Rig_z, y="depth")
fg.map_row("euc", median_Rig.where(median_Rig.num_obs > 30 * 24), y="z_above_euc")

fg.set_xlabels("$Ri_g$")
fg.set_ylabels(["depth [m]", "Height above $z_{EUC}$ [m]"])
fg.set_col_labels(["156°E", "165°E", "170°W", "140°W", "110°W"])
fg.axes[-1, -1].legend(fg.handles["euc"][-1], median_Rig.season.values)

[dcpy.plots.linex(0.25, lw=0.5, ax=ax) for ax in fg.axes.flat]
[dcpy.plots.liney(0, lw=0.5, ax=ax) for ax in fg.axes.flat]

fg.fig.set_size_inches((dcpy.plots.pub_fig_width("jpo", "two column"), 5))
dcpy.plots.label_subplots(fg.axes.flat, fontsize="small")
# fg.fig.savefig("../images/rig-profiles-depth-zeuc.png")

[Text(0.05, 0.9, '(a) '),
 Text(0.05, 0.9, '(b) '),
 Text(0.05, 0.9, '(c) '),
 Text(0.05, 0.9, '(d) '),
 Text(0.05, 0.9, '(e) '),
 Text(0.05, 0.9, '(f) '),
 Text(0.05, 0.9, '(g) '),
 Text(0.05, 0.9, '(h) '),
 Text(0.05, 0.9, '(i) '),
 Text(0.05, 0.9, '(j) ')]
_images/b1e3535892c14fd531a4f7f6485137e1b69a8e7e4f42070fee12533b08a18ca9.png 
with mpl.rc_context(sm13_colors):
    tao_seasonal.Rib.plot.line(hue="season", yscale="log", marker="o")
dcpy.plots.liney(1)
_images/7f24575884fe4d82e8e696aaa0f42d71225da613a3e56325604a825ff4e917d6.png

TAO seasonal variations#

with mpl.rc_context(sm13_colors):
    tao_seasonal.eucmax.plot.line(hue="season", marker="o")
_images/d027ed5e029044507733d5dea9d952c4f733f7c610eb7d4bbaee43e9b12c4b81.png 
f, ax = pump.plot.plot_bulk_Ri_diagnosis(
    tao_seasonal.sel(season="JJA"), ls="none", marker="s", **kwargs
)
_images/6703b21c5d038b11c5b09a2e479f2f4a2ad707ea6c1123001f6d1160de2ae3b4.png 
f, ax = pump.plot.plot_bulk_Ri_diagnosis(
    tao_seasonal.sel(season="MAM"), ls="none", marker="s", **kwargs
)
_images/24f9c256068b44ad929a88f259299c7a64177fe7ade7c7b5e42947317ace2fbb.png

fractional contributions#

tao_zeuc = tao_zeuc.where(tao_zeuc.count("time") > 30 * 24)
Ri_q_annual = (
    tao_zeuc.Rig.sel(zeuc=slice(0, None))
    .chunk({"time": -1, "zeuc": -1})
    .quantile(q=[0.25, 0.5, 0.75], dim=["time", "zeuc"])
)
Ri_q_annual.load()

<xarray.DataArray 'Rig' (quantile: 3, longitude: 5)>
array([[0.8068002 , 0.55212802, 0.25348482, 0.20803469, 0.28020857],
       [1.95963135, 1.54025941, 0.68420086, 0.39974787, 0.59899564],
       [5.42141063, 4.57124483, 2.60617434, 1.1664548 , 1.8581921 ]])
Coordinates:
  * longitude  (longitude) float64 -204.0 -195.0 -170.0 -140.0 -110.0
  * quantile   (quantile) float64 0.25 0.5 0.75
kwargs = {"buoy": False}

f, ax = pump.plot.plot_bulk_Ri_diagnosis(
    tao_clim.isel(longitude=slice(1, None)), ls="none", marker="s", **kwargs
)
# pump.plot.plot_bulk_Ri_diagnosis(tao_daily, ls="none", marker= 'o', f=f, ax=ax)

# take out 2 westernmost points
johnson_proc = johnson.sel(longitude=slice(-200, None))
johnson_proc["longitude"] = johnson_proc.longitude + 2.05

pump.plot.plot_bulk_Ri_diagnosis(
    johnson_proc, ls="none", marker="o", f=f, ax=ax, **kwargs
)
# ri.sel(zeuc=slice(0,None)).mean("zeuc").plot(ax=ax["Rib"], color='k')

# dcpy.plots.fill_between(
#    Ri_q_annual.sel(quantile=[0.25, 0.75]),
#    axis="y",
#    y="quantile",
#    x="longitude",
#    alpha=0.1,
#    ax=ax["Rib"],
#    color="k",
# )
# Ri_q_annual.sel(quantile=0.5).plot(color="k", ax=ax["Rib"], _labels=False)

ax["contrib"].set_xlim((-200, None))
ax["Rib"].set_title("")

ax["Rib"].legend(["TAO", "Johnson"], ncol=2)
f.set_size_inches(dcpy.plots.pub_fig_width("jpo", "two column"), 7)
# f.savefig("../images/bulk-ri-fractional-contrib.png")
_images/23573527d289ab3319cf895b4e27aea2bb564e44abac6e06c46c18de5e3ecb69.png

Old version#

Ri profiles#

with mpl.rc_context(sm13_colors):
    fg = median_Rig_z.squeeze().plot.line(
        col="longitude",
        hue="season",
        y="depth",
        ylim=[-150, 0],
        xlim=[0.1, 3.5],
        xscale="log",
    )
fg.map(lambda: dcpy.plots.linex([0.25], lw=0.75))
plt.gcf().set_size_inches((8, 4))
_images/aa87e965e43962be6a081f29afdcbee7c91cc8388dba320b889ba1b1714f8342.png 
with mpl.rc_context(sm13_colors):
    fg = median_Rig.where(median_Rig.num_obs > 30 * 24).plot.line(
        col="longitude",
        hue="season",
        y="z_above_euc",
        xlim=[1e-1, 3.5],
        ylim=[-150, 5],
        xscale="log",
    )

fg.map(lambda: dcpy.plots.linex(0.25, lw=0.75))
fg.fig.set_size_inches((10, 4))
_images/ff4430562da0002dbec79e717e9d451f9dce99f1b30f6c4025eb9a15bd2b1f12.png

Test MLD vs 25m for surface#

Ideally we would only see errors in \(h\), that seems to be the case

TAO#

Also not bad!

tao_clim = pump.calc.estimate_bulk_Ri_terms(tao_clim)
f, ax = pump.plot.plot_bulk_Ri_diagnosis(tao_clim, ls="none", marker="^")

tao_clim = pump.calc.estimate_bulk_Ri_terms(tao_clim, use_mld=False)
pump.plot.plot_bulk_Ri_diagnosis(tao_clim, ls="none", marker="o", f=f, ax=ax);
_images/1d67a1c6af749f938627c6b4549d2552de30ad17395ab346abb4a8ea9af6d1ff.png

Johnson#

Not bad once I take out the western most two points where the linear approx is bad. Before that:

  1. With MLD, it looked like Ri_b was (high, low, high) so dRi/dx was positive.

  2. Fixed by throwing out the westernmost two points.

Triangles and circles agree pretty closely so this is not bad

johnson = pump.calc.estimate_bulk_Ri_terms(johnson.isel(longitude=slice(2, None)))
f, ax = pump.plot.plot_bulk_Ri_diagnosis(johnson, ls="none", marker="^")

johnson = pump.calc.estimate_bulk_Ri_terms(johnson, use_mld=False)
pump.plot.plot_bulk_Ri_diagnosis(johnson, ls="none", marker="o", f=f, ax=ax);
_images/398713181b0dcd3946e082edc06ee0eb53af1a95456e8d00864c6c857aff4fa4.png