BoB ASIRI 2013 CTD χpod#

This notebook will attempt to infer eddy diffusivity \(K_e\) from the basin-wide transect of \(χ, K_T\)

(1)#\[\begin{equation} K_e = \frac{⟨\widetilde{χ}⟩/2 - ⟨K_T \, ∂_z\widetilde{θ}⟩ \; ∂_zθ_m}{|∇θ_m|²} \end{equation}\]
import eddydiff as ed

/home/deepak/anaconda3/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88
  return f(*args, **kwds)

Transect TS plots#

tr1 = xr.open_dataset("../datasets/bob-ctd-chipod/transect_1.nc", autoclose=True)
tr2 = xr.open_dataset("../datasets/bob-ctd-chipod/transect_2.nc", autoclose=True)
tr3 = xr.open_dataset("../datasets/bob-ctd-chipod/transect_3.nc", autoclose=True)

_, bins = pd.cut(tr1.rho.values.ravel(), 8, retbins=True)

dcpy.oceans.TSplot(tr1.S, tr1["T"], tr1.pres, rho_levels=bins - 1000)
plt.title("Transect 1")
plt.savefig("../images/bob-TS-transect-1.png", bbox_inches="tight")
dcpy.oceans.TSplot(tr2.S, tr2["T"], tr2.pres, rho_levels=bins - 1000)
plt.title("Transect 2")
plt.savefig("../images/bob-TS-transect-2.png", bbox_inches="tight")
dcpy.oceans.TSplot(tr3.S, tr3["T"], tr3.pres, rho_levels=bins - 1000)
plt.title("Transect 3")
plt.savefig("../images/bob-TS-transect-3.png", bbox_inches="tight")
_images/791d0798fa98b9c0271467ee6f2b6da39238f45c788922e0e66675333f537153.png _images/527d5a557e4605371b450d968f1dce17916f865267892fa03f84ae9394b5a61e.png _images/218ec9d71567187754bd81e48c9dc4b1b26822a57c5fc2fed50f961ba317ee34.png

Read and plot single transect#

transect = xr.open_dataset(
    "../datasets/bob-ctd-chipod/transect_1.nc", autoclose=True
).sel(cast=slice(10, 36))

transect["KtTz"] = transect["KT"] * transect["dTdz"]
transect["Jq"] = 1025 * 4200 * transect["KtTz"]

transect["KtTz"].values[np.abs(transect["Jq"].values) > 500] = np.nan
transect["Jq"].values[np.abs(transect["Jq"].values) > 500] = np.nan

np.log10(transect["KT"]).plot(robust=True, yincrease=False)

plt.figure()
transect["dTdz"].plot(robust=True, yincrease=False)

plt.figure()
np.log10(transect["KtTz"]).plot(robust=True, yincrease=False)

plt.figure()
transect["T"].sel(cast=slice(10, 36)).plot.contour(
    levels=40, colors="k", robust=True, yincrease=False
)
transect["rho"].sel(cast=slice(10, 36)).plot.contour(
    levels=20, colors="r", robust=True, yincrease=False
)
plt.title("black=T, red=ρ")

# transect = transect.drop(30, dim='cast').sel(P=slice(None, 110))

/home/deepak/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in greater
  import sys
/home/deepak/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in greater
  

Text(0.5,1,'black=T, red=ρ')
_images/3d5362712c0a555a397777c0efac6fa9632b8445c1a8ff277deda4912dc9971e.png _images/1ca9204d8c45e2fc1b3c30556e1dedb578b5cab46d5b73a434acef77c8d3d33f.png _images/0ce610e25bec61b5492daf2eead860e40979600d805031c06eeef1a9afac02ef.png _images/9d8597dfa6864458f0bdcef0f45b1433934617579f16a54ec0d795b68b3c2eaf.png

Make sure \(J_q\) isn’t crazy

np.log10(np.abs(transect["Jq"])).plot.hist()

(array([ 30., 179., 503., 586., 491., 301., 138.,  67.,  36.,  25.]),
 array([-3.48594309, -2.87431314, -2.2626832 , -1.65105325, -1.03942331,
        -0.42779336,  0.18383659,  0.79546653,  1.40709648,  2.01872642,
         2.63035637]),
 <a list of 10 Patch objects>)
_images/9fc0ba833ce5e8262498cec2c19b88232211929b4f5fc340939b029d74b720dc.png

Process transect data#

transect = xr.open_dataset(
    "../datasets/bob-ctd-chipod/transect_1.nc", autoclose=True
).sel(cast=slice(10, 36))

transect["KtTz"] = transect["KT"] * transect["dTdz"]
transect["Jq"] = 1025 * 4200 * transect["KtTz"]
transect["KtTz"].values[np.abs(transect["Jq"].values) > 1000] = np.nan
transect["Jq"].values[np.abs(transect["Jq"].values) > 1000] = np.nan

transect = transect.sel(P=slice(None, 210))

# def process_transect
# takes transect in z-space as input.

trmean_rho = ed.transect_to_density_space(transect)

/home/deepak/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in greater
  
/home/deepak/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in greater
  import sys

Read in other datasets#

cole = ed.read_cole()

argograd = xr.open_dataset(
    "../datasets/argo_dec_iso_gradients.nc", decode_times=False, autoclose=True
).load()

eccograd = xr.open_dataset(
    "../datasets/ecco_monthly_iso_gradient.nc", decode_times=False, autoclose=True
).load()

Using climatological gradients#

eccoKe = ed.process_transect_1d(transect, eccograd, "ECCO transect1")
argoKe = ed.process_transect_1d(transect, argograd, "ARGO transect1")

Ke = pd.DataFrame()
Ke["ecco"] = eccoKe.Ke
Ke["argo"] = argoKe.Ke
# Ke['rho'] = Ke.ecco.index.mid.astype('float')
eccoKe
chi KtTz dTdz dTmdz dTiso KT Ke
rho
(1019.0, 1020.5] 1.087531e-07 3.332968e-06 0.017668 -0.019127 2.065550e-06 0.000894 -2196.788029
(1020.5, 1021.0] 4.242331e-08 1.787747e-06 0.023994 -0.015216 3.913078e-06 0.003574 -391.230758
(1021.0, 1022.4] 1.383193e-07 7.649035e-07 0.099101 0.010196 1.521268e-06 0.000010 26514.376375
(1022.4, 1023.5] 7.467831e-08 3.012975e-07 0.131284 0.106301 1.086169e-06 0.000003 4501.655168
(1023.5, 1024.2] 1.822328e-07 7.604694e-07 0.121188 0.130226 7.190793e-07 0.000006 -15309.981990
(1024.2, 1024.9] 2.110117e-07 1.077232e-06 0.107393 0.113445 4.345696e-07 0.000011 -88435.950761
(1024.9, 1025.4] 1.015220e-07 5.481218e-07 0.093691 0.091905 2.953978e-07 0.000006 4423.677911
(1025.4, 1025.8] 1.114811e-07 6.400163e-07 0.079396 0.082158 2.854238e-07 0.000007 38766.088873
(1025.8, 1026.1] 1.349486e-07 8.839158e-07 0.060302 0.077620 2.923084e-07 0.000105 -13281.205312
(1026.1, 1026.5] 2.656529e-08 4.195098e-07 0.041820 0.076793 2.593703e-07 0.000235 -281430.091207 
/home/deepak/anaconda3/lib/python3.6/site-packages/_pytest/fixtures.py:847: DeprecationWarning: The `convert` argument is deprecated in favor of `converter`.  It will be removed after 2019/01.
  params = attr.ib(convert=attr.converters.optional(tuple))
/home/deepak/anaconda3/lib/python3.6/site-packages/_pytest/fixtures.py:849: DeprecationWarning: The `convert` argument is deprecated in favor of `converter`.  It will be removed after 2019/01.
  ids = attr.ib(default=None, convert=_ensure_immutable_ids)

ed.plot_bar_Ke(eccoKe)
# ed.plot_bar_Ke(argoKe)
_images/116466a491697f37c968e6c35a573edb014e9d556e599d5a46bb371511b26780.png 
plt.figure()
argomean.dTdia.plot(x="rho")
eccomean.dTdia.plot(x="rho")
(trmean.dTdz).plot(x="rho")
plt.legend(["argo", "ecco", "transect"])

plt.figure()
argomean.dTiso.plot(x="rho")
eccomean.dTiso.plot(x="rho")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-30-9f1417f79a7b> in <module>()
      1 plt.figure()
----> 2 argomean.dTdia.plot(x='rho')
      3 eccomean.dTdia.plot(x='rho')
      4 (trmean.dTdz).plot(x='rho')
      5 plt.legend(['argo', 'ecco', 'transect'])

NameError: name 'argomean' is not defined

<Figure size 780x660 with 0 Axes>

(
    ecco.dTiso.sel(lat=slice(transect.lat.min(), transect.lat.max()))
    .sel(lon=transect.lon.mean(), method="nearest")
    .plot.contourf(robust=True, yincrease=False)
)

(
    ecco.Tmean.sel(lat=slice(transect.lat.min(), transect.lat.max()))
    .sel(lon=transect.lon.mean(), method="nearest")
    .plot.contour(robust=True, colors="w", yincrease=False)
)

plt.gca().set_ylim([200, 0])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-12-04509405bfa7> in <module>()
----> 1 (ecco.dTiso.sel(lat=slice(transect.lat.min(), transect.lat.max()))
      2  .sel(lon=transect.lon.mean(), method='nearest')
      3  .plot.contourf(robust=True, yincrease=False))
      4 
      5 (ecco.Tmean.sel(lat=slice(transect.lat.min(), transect.lat.max()))

NameError: name 'ecco' is not defined

Ferrari & Polzin (2005) approach#

They average all data on density surfaces, estimate mean depth of isopycnal \(z_n\) and calculate diapycnal gradient. The result is a vertical profile of eddy diffusivity. All spatial variations are averaged over.

They calculate an isopycnal gradient by doing a plane fit along isopycnals to all the data. This is required for the finestructure estimate

transKe

<xarray.Dataset>
Dimensions:  (dist: 26, rho: 14)
Coordinates:
  * dist     (dist) float64 0.0 37.02 74.47 111.3 148.0 185.1 222.2 259.1 ...
  * rho      (rho) float64 1.02e+03 1.02e+03 1.021e+03 1.021e+03 1.022e+03 ...
Data variables:
    Ke       (dist, rho) float64 nan nan nan -293.5 -517.0 -2.837e+03 ...
    KT       (dist, rho) float64 nan nan 0.0003355 6.889e-08 2.743e-08 ...
    KtTz     (dist, rho) float64 nan nan 1.664e-06 4.679e-09 2.365e-09 ...
    dTdz     (dist, rho) float64 nan nan 0.004763 0.06767 0.1226 0.1494 ...
    chi      (dist, rho) float64 nan nan 1.865e-08 6.313e-10 3.188e-10 ...
    dTmdz    (dist, rho) float64 nan nan nan 0.1624 0.3492 0.3736 0.1541 ...
    dTiso    (dist, rho) float64 nan nan -1.973e-06 -1.23e-06 1.135e-06 ...

# def transect_fp(trdens):

trdens = ed.to_density_space(transect)

dTiso = xr.DataArray(
    np.ones_like(trdens.rho) * np.nan, dims=["rho"], coords={"rho": trdens.rho}
)

# fit straight lines to get mean gradients
# FP do a plane fit in 2D
for idx, rr in enumerate(trdens.rho):
    Tvec = trdens["T"].sel(rho=rr)
    mask = ~np.isnan(Tvec)

    if mask.sum() < 3:
        continue

    dTiso[idx], _, _, _, _ = sp.stats.linregress(trdens.dist[mask] * 1000, Tvec[mask])

trmean = trdens.mean(dim="cast")
trmean

Tmean = ed.to_depth_space(
    trmean["T"].expand_dims(["cast"]), Pold=trmean["P"].expand_dims(["cast"])
).mean(dim="cast")
rhomean = ed.to_depth_space(
    trmean["rho"].expand_dims(["cast"]), Pold=trmean["P"].expand_dims(["cast"])
).mean(dim="cast")

_, Tsmooth = ed.fit_spline(Tmean.P, Tmean, k=2)

# Tmean.plot()
dTdz = ed.xgradient(-Tsmooth, dim="P")
dTdz.coords["rho"] = rhomean

dTdz = dTdz.expand_dims(["cast"])
dTdz["rho"] = dTdz.rho.expand_dims(["cast"])

dTmdz = ed.to_density_space(dTdz).mean(dim="cast")

fpKe = xr.Dataset()
fpKe["chi"] = trmean.chi
fpKe["dTmdz"] = trmean.dTdz
fpKe["dTiso"] = dTiso
fpKe["KtTz"] = trmean.KtTz
fpKe["KT"] = trmean.KT
fpKe["dTdz"] = trmean.dTdz

fpKe["Ke"] = np.abs(fpKe.chi / 2 - fpKe.KtTz * fpKe.dTmdz) / fpKe.dTiso**2

ed.plot_bar_Ke(fpKe.to_dataframe())

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-323-721f90866cb5> in <module>()
      1 # def transect_fp(trdens):
      2 
----> 3 trdens = ed.to_density_space(transect)
      4 
      5 dTiso = xr.DataArray(np.ones_like(trdens.rho) * np.nan,

~/work/eddydiff/eddydiff/eddydiff.py in to_density_space(da, rhonew)
    199     in_dens = xr.Dataset()
    200 
--> 201     for vv in da.variables:
    202         if vv in da.coords or vv == 'rho':
    203             continue

~/anaconda3/lib/python3.6/site-packages/pandas/core/generic.py in __getattr__(self, name)
   4370             if self._info_axis._can_hold_identifiers_and_holds_name(name):
   4371                 return self[name]
-> 4372             return object.__getattribute__(self, name)
   4373 
   4374     def __setattr__(self, name, value):

AttributeError: 'DataFrame' object has no attribute 'variables'

Without interpolating to convert to density space#

mask = ~np.isnan(trmean_rho["T"].isel(rho=5))
ed.fit_spline(
    trmean_rho["dist"][mask].values,
    trmean_rho["T"].isel(rho=5)[mask].values,
    debug=True,
)

(<scipy.interpolate.fitpack2.UnivariateSpline at 0x7f917a7ad6d8>,
 array([26.25453411, 26.27269832, 26.26237948, 26.22850841, 26.17573395,
        26.10729261, 26.02839219, 25.94359928, 25.85682278, 25.77302511,
        25.69583305, 25.63065117, 25.58101885, 25.55179105, 25.54613198,
        25.56238605, 25.59792679, 25.64976453, 25.7146301 , 25.79062798,
        25.79063393, 25.87376963, 25.95664316, 26.02779096, 26.08732754,
        26.10676851, 26.07836197]))
_images/a9b789acbd763cc6884b38569bc669c23d4474303f7e9e63629adc878905dd03.png 
# Run Section 1.3 first
# 2. Generate mean field T, P using cubic splines to smooth
#    Calculate Trms as RMS of deviations from smoothed T field

trmean_rho["Tmean"] = ed.smooth_cubic_spline(trmean_rho["T"])
trmean_rho.Tmean.attrs["name"] = "Smoothed temperature"

trmean_rho["Pmean"] = ed.smooth_cubic_spline(trmean_rho["P"])
trmean_rho.Pmean.attrs["name"] = "Smoothed pressure"


def check_spline_smoothing(T, Ts, P, Ps, Traw, Praw):
    f, ax = plt.subplots(3, 1)

    T.plot.contour(ax=ax[0], y="rho", levels=50, colors="r")
    Ts.plot.contour(ax=ax[0], y="rho", levels=50, colors="k", yincrease=False)

    P.plot.contour(ax=ax[1], y="rho", levels=50, colors="r")
    Ps.plot.contour(ax=ax[1], y="rho", levels=50, colors="k", yincrease=False)

    # cmat, rhomat = xr.broadcast(T.cast, T.rho)
    # ax[2].contour(cmat, P, T, levels=50, colors='r')

    Tz = ed.to_depth_space(T, P)
    Tzs = ed.to_depth_space(Ts, Ps)

    ax[2].contour(
        xr.broadcast(Traw.cast, Praw)[0], Praw.T, Traw.T, 50, cmap=mpl.cm.BuGn
    )
    Tz.plot.contour(ax=ax[2], x="cast", levels=30, colors="r")
    Tzs.plot.contour(ax=ax[2], x="cast", levels=30, colors="k", yincrease=False)

    f.set_size_inches((8, 10))
    ax[0].set_title("Smoothed temp")
    ax[1].set_title("Smoothed pressure")
    ax[2].set_title("Raw and smoothed temperature in depth space")
    ax[2].set_ylim([220, 0])
    plt.tight_layout()


check_spline_smoothing(
    trmean_rho["T"],
    trmean_rho["Tmean"],
    trmean_rho["P"],
    trmean_rho["Pmean"],
    transect["T"],
    transect["pres"],
)

trmean_rho["Trms"] = np.sqrt(
    ((trmean_rho["T"] - trmean_rho["Tmean"]) ** 2).mean(dim="cast")
)
trmean_rho.Trms.attrs["name"] = "RMS temp variations"

/home/deepak/anaconda3/lib/python3.6/site-packages/numpy/core/_methods.py:29: RuntimeWarning: invalid value encountered in reduce
  return umr_minimum(a, axis, None, out, keepdims)
/home/deepak/anaconda3/lib/python3.6/site-packages/numpy/core/_methods.py:26: RuntimeWarning: invalid value encountered in reduce
  return umr_maximum(a, axis, None, out, keepdims)
_images/100c2dc51a6702814f9d74aca4a6a344de7f3f4f7b62a7876be89cf4203205e4.png 
# 3. Calculate isopyncal and diapycnal gradients of the mean field
trmean_rho["dTiso"], trmean_rho["dTdia"] = ed.calc_iso_dia_gradients(
    trmean_rho["Tmean"], trmean_rho["Pmean"], debug=True
)
trmean_rho.dTiso.attrs["name"] = "Isopycnal ∇T"
trmean_rho.dTdia.attrs["name"] = "Diapycnal ∇T"
_images/e4ff9be4bd17f24e749a8e101a5273fd12f1e21c4daac640f69ab498edeaa899.png 
transKe = xr.Dataset()
transKe["KT"] = trmean_rho.KT
transKe["KtTz"] = trmean_rho.KtTz
transKe["dTdz"] = trmean_rho.dTdz
transKe["chi"] = trmean_rho.chi
transKe["dTmdz"] = trmean_rho.dTdia
transKe["dTiso"] = trmean_rho.dTiso
transKe["rho"].values = np.round(transKe.rho.values, decimals=2)
transKe["Tm"] = trmean_rho["Tmean"]
transKe.attrs["name"] = "BoB CTD χpod"


def calc_Ke(transKe, navg=None):

    transKe = transKe.copy()

    if navg is not None and np.isinf(navg):
        transKe = transKe.mean(dim="cast")

    elif navg is not None:
        cbins = np.arange(transKe.cast.min(), transKe.cast.max() + 1, navg)
        transKe = transKe.groupby_bins("cast", cbins, labels=cbins[:-1] + 1).mean(
            dim="cast"
        )

    transKe["Ke"] = (
        transKe.chi / 2 - transKe.KtTz * transKe.dTmdz
    ) / transKe.dTiso**2
    transKe.Ke.values[np.abs(transKe.dTiso.values) < 5e-7] = np.nan
    transKe.Ke["name"] = "$K_e$"

    transKe.attrs["navg"] = navg

    if "cast_bins" in transKe.coords:
        transKe = transKe.rename({"cast_bins": "cast"})

    return transKe


transKe = calc_Ke(transKe)

ed.plot_transect_Ke(transKe)

/home/deepak/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in less
_images/043fb2c34262b6890684a3fc357c04f519e392366dc19de0df173f59f3f68d17.png 
transKe2 = calc_Ke(transKe, navg=None)
ax, axback = ed.plot_transect_var(
    x="cast",
    y="rho",
    data=transKe2.Ke,
    fill=trmean_rho["dTiso"],
    contour=trmean_rho["P"],
    xlim=[1, 1e4],
    xticks=[1e3],
)
ax, axback = ed.plot_transect_var(
    x="cast",
    y="rho",
    data=transKe2.chi / 2,
    bar2=transKe2.KtTz * transKe2.dTmdz,
    fill=trmean_rho["T"],
    contour=trmean_rho["P"],
)

/home/deepak/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in less
/home/deepak/anaconda3/lib/python3.6/site-packages/matplotlib/contour.py:1173: UserWarning: No contour levels were found within the data range.
  warnings.warn("No contour levels were found"
_images/e6665114d38aa90446a37dfc65fbaaa2fd8100d1eb7e6217087b0a11a3224cb1.png _images/04d4effdeb3df90a4d77280e7cd6720a1ff50227da5260999d7386307095ef29.png 
transKe2.Ke.plot.line(y="rho", hue="cast", yincrease=False)
plt.gca().set_xscale("log")
plt.gca().set_xlim([1e2, 1e5])

cole = ed.read_cole()

region = ed.get_region_from_transect(transect)

cole_bay = (cole.sel(lat=region["lat"])).mean(dim="lon").mean(dim="lat")
cole_bay["density_mean_depth"] += 1000
cole_bay = cole_bay.set_coords("density_mean_depth")
cole_bay.diffusivity.plot(y="density_mean_depth", color="k", yincrease=False)

[<matplotlib.lines.Line2D at 0x7f91bf98f860>]
_images/d040a0aabfb4107262f697c824d23adf6782d43edd767d1318d4bfdd9b41555d.png 
transKe2.Ke.where(transKe2.Ke > 0).mean(dim="cast").plot.line(y="rho")
cole_bay = (cole.sel(lat=region["lat"])).mean(dim="lon").mean(dim="lat")
cole_bay["density_mean_depth"] += 1000
cole_bay = cole_bay.set_coords("density_mean_depth")
cole_bay.diffusivity.plot(y="density_mean_depth", color="k", yincrease=False)

plt.gca().set_xscale("log")
plt.gca().set_xlim([1e2, 1e5])
plt.figlegend(["Mean in along-transect direction", "Cole et al (2015)"])

<matplotlib.legend.Legend at 0x7f91bfb673c8>
_images/103d430040e296ff84913ec0309773341b60c68a2993f5a64d0ab9d5e4700dd7.png

Compare various mean fields / gradients#

argo = xr.open_dataset(
    "../datasets/argo_clim_iso_gradients.nc", decode_times=False, autoclose=True
).load()

ecco = xr.open_dataset("../datasets/ecco_annual_iso_gradient.nc", autoclose=True).load()

region = ed.get_region_from_transect(transect)

argo.dTiso.sel(**region).mean(dim="lon").plot(robust=True, yincrease=False)
plt.figure()
ecco.dTiso.sel(**region).mean(dim="lon").plot(robust=True, yincrease=False)

plt.figure()
np.abs(trmean_rho.dTiso).plot(robust=True, yincrease=False)

<matplotlib.collections.QuadMesh at 0x7f1c7a70b9e8>
_images/63d313ed6cc6dd9734cf6ceaa13d8a5e27fb66d4cfb29e88637cecfbaf546690.png _images/c96a601eb1bc1a23acac8b632057d48f41dd2efd7ea3de3cacefab4ed27fd359.png _images/787f0841a021a0902f4b462cbd433bf57ee714de4ae219e0abc5e5ba42e210b4.png

Argo estimate (Cole et al, 2015)#

region = ed.get_region_from_transect(transect)

cole_bay = cole.sel(lat=region["lat"], lon=region["lon"])

(
    cole_bay.diffusivity.mean(dim="lon").plot(
        y="depth", yincrease=False, norm=mpl.colors.LogNorm(), cmap=mpl.cm.Reds
    )
)
plt.gca().set_ylim([200, 0])

(200, 0)

/home/deepak/anaconda3/lib/python3.6/site-packages/matplotlib/colors.py:1031: RuntimeWarning: invalid value encountered in less_equal
  mask |= resdat <= 0
_images/0c032117bdb699245e4efd078ae5298caa2faae74189feac6747cdb26dc99380.png 
(
    cole_bay.diffusivity.mean(dim="lat")
    .mean(dim="lon")
    .plot.line(y="depth", yincrease=False)
)
plt.gca().set_xscale("log")
plt.gca().set_xlim([1e2, 3e4])
plt.gca().set_ylim([200, 0])
plt.gca().grid(True)
_images/5186383d545f09eefe7308194d53a3abc32cd67774ed0b9a558c76576ee86bbf.png 
(
    cole_bay.diffusivity.groupby_bins(
        cole_bay.density_mean_depth + 1000, ρbins, labels=(ρbins[:-1] + ρbins[1:]) / 2
    )
    .mean()
    .plot.line(y="density_mean_depth_bins", yincrease=False)
)

[<matplotlib.lines.Line2D at 0x7f91cf175dd8>]
_images/b8c8a2032606f147c58260b5931452b63e88bac160e65f5f30ecaeab07970392.png

All data combined#

Let’s ignore 3 for now, not sure why that overlaps.

transect1 = xr.open_dataset("../datasets/bob-ctd-chipod/transect_1.nc", autoclose=True)

transect2 = xr.open_dataset("../datasets/bob-ctd-chipod/transect_2.nc", autoclose=True)

transect3 = xr.open_dataset("../datasets/bob-ctd-chipod/transect_3.nc", autoclose=True)

plt.plot(transect1.lon, transect1.lat, ".", ms=20)
plt.plot(transect2.lon, transect2.lat, ".", ms=20)
plt.plot(transect3.lon, transect3.lat, "o")
plt.legend(("big(1)", "big(2)", "big(3)"))

<matplotlib.legend.Legend at 0x7fb11cb45390>
_images/c65a6c4c682a4f41b416040598c6bb380c02682b1a3a4959385eaffe09243b4e.png 
# merge transect 1 and transect2
transect = pd.concat(
    [tr.to_dataframe().reset_index() for tr in [transect1.sel(P=slice(0, 120))]]
)

transect["KtTz"] = transect["KT"] * transect["dTdz"]
transect["Jq"] = 1025 * 4200 * transect["KtTz"]

mask = np.logical_or(np.abs(transect["Jq"].values) > 2000, transect["KT"].values > 5e-3)

transect["KtTz"].values[mask] = np.nan
transect["Jq"].values[mask] = np.nan
transect["chi"].values[mask] = np.nan

eccoKe2 = ed.process_transect_1d(transect, eccograd, "ECCO")
argoKe2 = ed.process_transect_1d(transect, argograd, "ARGO")

ed.plot_bar_Ke(eccoKe2)
ed.plot_bar_Ke(argoKe2)

eccoKe2.Ke

/home/deepak/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in greater
  
/home/deepak/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:9: RuntimeWarning: invalid value encountered in greater
  if __name__ == '__main__':

rho
(1019.0, 1020.6]      3938.569746
(1020.6, 1021.2]     10729.992555
(1021.2, 1022.7]     46947.664696
(1022.7, 1023.9]     -2230.156607
(1023.9, 1024.7]     -5942.732091
(1024.7, 1025.4]     48277.446666
(1025.4, 1026.0]     89531.725066
(1026.0, 1026.5]    316205.497546
(1026.5, 1026.9]     47694.340619
(1026.9, 1027.5]    -50417.990562
Name: Ke, dtype: float64
_images/9eeadf8436cb3f3018c0999ef5cb2d02f53531306bbb224df5d268a8007de1c0.png _images/2bb05595846e4f209b3c77d46ff8f491a1feaac73248676cbd794b597733559e.png 
transect1["T"].sel(cast=slice(9, 40)).plot(cmap=mpl.cm.RdYlBu_r, yincrease=False)
transect1["rho"].sel(cast=slice(9, 40)).plot.contour(
    colors="k", levels=np.arange(1020, 1028), yincrease=False
)
plt.figure()

(
    eccograd["ρmean"]
    .sel(**ed.get_region_from_transect(transect1))
    .mean(dim="lon")
    .plot.contour(levels=np.arange(1020, 1028), colors="k", yincrease=False)
)
plt.gca().set_ylim([200, 0])

(200, 0)
_images/aa669f7420d6fc6bd538d91bd3e650beeada5fe16b7b30588902183d243c3a01.png _images/0f0b0f350e7e8a201d4e7e40fbb3417cf8d3f9496ad8990ddb2497b6ff25634b.png 
trmean, ρbins = ed.average_transect_1d(transect)

eccomean = ed.average_clim(eccograd, transect, ρbins)

eccomean
lat lon pres Tmean Smean RHOAnoma ρmean dTiso dTdia dSiso dSdia dTdz_local dSdz Pmean dTdz
ρmean
(1019.1, 1020.7] 16.507143 87.464286 10.333333 28.939016 32.836595 -8.500231 1020.499769 7.333585e-07 0.005335 1.824942e-06 0.013389 0.005335 -0.013389 10.333333 0.032508
(1020.7, 1021.2] 12.984742 87.032864 14.647887 28.798760 33.449447 -7.972470 1021.027530 6.038192e-07 0.011383 9.272963e-07 0.009727 0.011383 -0.009727 14.647887 0.032730
(1021.2, 1022.5] 12.422030 86.626856 39.381188 27.957727 33.956771 -7.209345 1021.790655 8.870507e-07 0.054069 5.606805e-07 0.020746 0.054069 -0.020746 39.381188 0.053711
(1022.5, 1023.4] 12.877078 86.702494 67.187101 25.848192 34.483206 -6.026892 1022.973108 9.335261e-07 0.107467 4.611139e-07 0.014878 0.107467 -0.014878 67.187101 0.100398
(1023.4, 1024.2] 13.092742 86.763441 84.984611 23.781918 34.670366 -5.187314 1023.812686 7.491000e-07 0.123242 3.328542e-07 0.008024 0.123242 -0.008024 84.984611 0.126546
(1024.2, 1024.7] 12.621930 86.674561 98.163965 22.012167 34.765387 -4.548118 1024.451882 6.128813e-07 0.117844 2.484168e-07 0.004501 0.117844 -0.004500 98.163965 0.120330
(1024.7, 1025.2] 12.814516 86.855572 113.221291 20.440335 34.814054 -4.011924 1024.988076 4.669423e-07 0.104717 1.765324e-07 0.003152 0.104717 -0.003152 113.221291 0.103577
(1025.2, 1025.7] 12.620301 86.759398 127.308123 18.991971 34.855904 -3.536118 1025.463882 3.916333e-07 0.093500 1.250950e-07 0.002424 0.093500 -0.002424 127.308123 0.111689
(1025.7, 1026.0] 11.579341 86.423653 136.347489 17.930917 34.884922 -3.205289 1025.794711 4.132742e-07 0.082132 9.679136e-08 0.001976 0.082132 -0.001976 136.347489 0.101579
(1026.0, 1026.5] 12.843915 86.653439 157.727486 16.558279 34.917639 -2.753373 1026.246627 3.859926e-07 0.072626 6.896777e-08 0.001861 0.072626 -0.001861 157.727486 0.064202 
sp, Tsm = ed.fit_spline(eccomean["Pmean"], eccomean["Tmean"], k=3, debug=False)
plt.plot(sp.derivative(1)(eccomean["Pmean"]))
plt.plot(-eccomean.dTdz.values)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-22-79acc47f9c17> in <module>()
----> 1 sp, Tsm = ed.fit_spline(eccomean['Pmean'], eccomean['Tmean'], k=3, debug=False)
      2 plt.plot(sp.derivative(1)(eccomean['Pmean']))
      3 plt.plot(-eccomean.dTdz.values)

NameError: name 'eccomean' is not defined

Lessons learned#

Don’t interpolate when converting to density space#

+The fields and gradients are different!+

FALSE ALARM: This is because transect.P and transect.pres are different!

# transform to depth space, calculate gradient and transform back
Tsmooth = ed.to_depth_space(Tdens_i, Pold=Pdens_i, Pnew=None)
dTmdz = -ed.xgradient(Tsmooth, "P")
ed.to_density_space(dTmdz).plot(yincrease=False)

# original calculation ΔT and ΔP are estimated in density space
plt.figure()
trmean_rho.dTdia.plot(yincrease=False, x="cast")

<matplotlib.collections.QuadMesh at 0x7fbee9919be0>
_images/39f06c179c4e9f9892016aea3ed5343b8853c72960de252f063a19946d471071.png _images/ca84ab97301d4ab6b8e093dec2600a912536d932c6f1c80337150691900f8185.png

Spline-smooth both temperature and pressure#

Below I calculate dT/dP in density space and compare that to dT/dP calculated by

  1. transforming spline smoothed T to pressure space

  2. differentiate to get dT/dP

  3. convert back to density space and plot.

trdens = ed.to_density_space(transect)

plt.close("all")

Tdens_i = ed.smooth_cubic_spline(trdens["T"], False)
Pdens_i = ed.smooth_cubic_spline(trdens["P"], False)

levels = np.linspace(trdens["T"].min(), trdens["T"].max(), 40)
trdens["T"].plot.contour(colors="r", levels=levels)
Tdens_i.plot.contour(colors="k", levels=levels, yincrease=False)
plt.title("In ρ space")

plt.figure()
Tsmooth = ed.to_depth_space(Tdens_i, Pold=Pdens_i, Pnew=None)
transect["T"].plot.contour(colors="r", levels=levels, yincrease=False)
Tsmooth.plot.contour(colors="k", levels=levels, yincrease=False)
plt.title("To depth space using spline smoothed pressure")
plt.ylim([120, 0])

plt.figure()
psmooth = dcpy.ts.xfilter(trdens["P"], dim="cast", flen=10)
Tsmooth = ed.to_depth_space(Tdens_i, Pold=psmooth, Pnew=None)
transect["T"].plot.contour(colors="r", levels=levels, yincrease=False)
Tsmooth.plot.contour(colors="k", levels=levels, yincrease=False)
plt.title("To depth space using hann smoothed pressure")
plt.ylim([120, 0])


plt.figure()
Tsmooth = ed.to_depth_space(Tdens_i, Pold=trdens["P"], Pnew=None)
transect["T"].plot.contour(colors="r", levels=levels, yincrease=False)
Tsmooth.plot.contour(colors="k", levels=levels, yincrease=False)
plt.title("To depth space using unsmoothed pressure")
plt.ylim([120, 0])

(120, 0)
_images/0af728e3887f6fa9e6b04d6f0dfa32ba10134b7be26ac9f6ce3b46c217ce1faa.png _images/c7ed410425d6623a2a1245f0ea72910b32edb7ffef1d1c26447468ac41e4636b.png _images/308bca41183ad14f1d6c74e18ddf74bd45cb796702d2c97c3d2db05565e436c3.png _images/b2a4e6d429f0d59f39903d98043ea7068e3e5e374b585c3e370e28e5bf5f0aec.png

Groupby on dataframe or xarray consistently#

xarray seems to use a different kind of index with groupby_bins.

Test gradients in isopycnal planes#

coords = {"cast": np.linspace(1, 50, 50), "P": np.linspace(0, 500, 200)}

T = xr.DataArray(
    np.ones((len(coords["P"]), len(coords["cast"]))) * np.nan,
    dims=["P", "cast"],
    coords=coords,
    name="T",
)

T = -0.6 * T.P + 0.05 * T.cast
T.name = "T"

rho = 1025 * (1 - 1.7e-4 * (T - 15))
rho.name = "$ρ$"

T.plot.contourf(levels=20)
rho.plot.contour(colors="k", yincrease=False)

dT = ed.gradient(T.rename({"P": "z", "cast": "x"}))
dT["dy"] = xr.zeros_like(dT["dx"])
dT.attrs["name"] = "dT"

drho = ed.gradient(rho.rename({"P": "z", "cast": "x"}))
drho["dy"] = xr.zeros_like(drho["dx"])
drho.attrs["name"] = "dρ"

dT

<xarray.Dataset>
Dimensions:  (x: 50, z: 200)
Coordinates:
  * z        (z) float64 0.0 2.513 5.025 7.538 10.05 12.56 15.08 17.59 20.1 ...
  * x        (x) float64 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 ...
Data variables:
    dz       (z, x) float64 -0.6 -0.6 -0.6 -0.6 -0.6 -0.6 -0.6 -0.6 -0.6 ...
    dx       (z, x) float64 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 ...
    mag      (z, x) float64 0.6021 0.6021 0.6021 0.6021 0.6021 0.6021 0.6021 ...
    dy       (z, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
Attributes:
    name:     dT
_images/df72893015485f7ce1608c86afc212a7a486c8d1d8ee2af3cdbf64bf5abd7797.png

Convert transect .mat to netcdf files#

ed.convert_mat_to_netcdf()

/home/deepak/anaconda3/lib/python3.6/_collections_abc.py:743: FutureWarning: iteration over an xarray.Dataset will change in xarray v0.11 to only include data variables, not coordinates. Iterate over the Dataset.variables property instead to preserve existing behavior in a forwards compatible manner.
  for key in self._mapping:
../eddydiff/eddydiff.py:571: RuntimeWarning: invalid value encountered in less
  mask2d = np.logical_or(transect['T'].values < 1,
../eddydiff/eddydiff.py:572: RuntimeWarning: invalid value encountered in less
  transect['S'].values < 1)
/home/deepak/anaconda3/lib/python3.6/_collections_abc.py:743: FutureWarning: iteration over an xarray.Dataset will change in xarray v0.11 to only include data variables, not coordinates. Iterate over the Dataset.variables property instead to preserve existing behavior in a forwards compatible manner.
  for key in self._mapping:
../eddydiff/eddydiff.py:571: RuntimeWarning: invalid value encountered in less
  mask2d = np.logical_or(transect['T'].values < 1,
../eddydiff/eddydiff.py:572: RuntimeWarning: invalid value encountered in less
  transect['S'].values < 1)
/home/deepak/anaconda3/lib/python3.6/_collections_abc.py:743: FutureWarning: iteration over an xarray.Dataset will change in xarray v0.11 to only include data variables, not coordinates. Iterate over the Dataset.variables property instead to preserve existing behavior in a forwards compatible manner.
  for key in self._mapping:
../eddydiff/eddydiff.py:571: RuntimeWarning: invalid value encountered in less
  mask2d = np.logical_or(transect['T'].values < 1,
../eddydiff/eddydiff.py:572: RuntimeWarning: invalid value encountered in less
  transect['S'].values < 1)