P06 CTD χpod analysis

P06 CTD χpod analysis#

import dcpy
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import seawater as sw

import eddydiff as ed
import xarray as xr

Read CTD χpod data and convert to netCDF#

ed.sections.to_netcdf(
    infile="/home/deepak/work/eddydiff/datasets/P06/P06.mat",
    outfile="/home/deepak/work/eddydiff/datasets/P06/p06-try-again.nc",
    transect_name="P06",
)

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-102-a748ec0f2bb7> in <module>
      2     infile="/home/deepak/work/eddydiff/datasets/P06/P06.mat",
      3     outfile="/home/deepak/work/eddydiff/datasets/P06/p06-try-again.nc",
----> 4     transect_name="P06",
      5 )

~/work/eddydiff/eddydiff/sections.py in to_netcdf(infile, outfile, transect_name)
     24     import seawater as sw
     25 
---> 26     mat = loadmat(infile, squeeze_me=True)["grd"]
     27 
     28     ctd = xr.Dataset()

KeyError: 'grd'

matold = loadmat("/home/deepak/work/eddydiff/datasets/P06/P06-means.mat")
matold.keys()

dict_keys(['__header__', '__version__', '__globals__', 'lat', 'lon', 'dnum', 'TPvar', 'eps', 'chi', 'KT', 'dTdz', 'N2', 't', 's', 'P', 'Jq', 'dist'])

ed.sections.to_netcdf(
    infile="/home/deepak/work/eddydiff/datasets/P06/P06-new.mat",
    outfile="/home/deepak/work/eddydiff/datasets/P06/p06-new.nc",
    transect_name="P06",
)

Read transect + ancillary datasets#

import eddydiff as ed

p06 = xr.open_dataset("/home/deepak/work/eddydiff/datasets/P06/p06.nc")
bathy = p06.bathy.copy()
p06 = p06.where(p06["KT"] < 3e-3)
p06["bathy"] = bathy
p06["KtTz"] = p06.KT * p06.dTdz
p06

eccograd, argograd, cole, aber = ed.read_all_datasets(kind="monthly", transect=p06)

/home/deepak/miniconda3/envs/dcpy/lib/python3.7/site-packages/ipykernel_launcher.py:3: FutureWarning: The autoclose argument is no longer used by xarray.open_dataset() and is now ignored; it will be removed in a future version of xarray. If necessary, you can control the maximum number of simultaneous open files with xarray.set_options(file_cache_maxsize=...).
  This is separate from the ipykernel package so we can avoid doing imports until
/home/deepak/miniconda3/envs/dcpy/lib/python3.7/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide
  x = np.divide(x1, x2, out)

argograd.Smean.plot(x="cast", yincrease=False)

<matplotlib.collections.QuadMesh at 0x7f173e54d470>

/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)

_images/37af54c526b08c2d3a8313867d2ce0ee895373ec685547b1363bbc28aafab6dc.png

plots#

np.log10(p06.eps).plot(x="cast")
p06.rho.plot.contour(
    colors="k",
    x="cast",
    yincrease=False,
    levels=pd.qcut(p06.rho.values.ravel(), 20, retbins=True)[1],
)
plt.gca().fill_between(p06.cast, 6100, -p06.bathy, color="k", zorder=10)
plt.gcf().set_size_inches(10, 10 / 1.6)

_images/a8993c8af6b6c4a3fb9921463fd766c86d288c2ed1ccfeb696d50a2cdf667af7.png

p06["T"].plot(x="cast", cmap=mpl.cm.RdYlBu)
p06.rho.plot.contour(
    colors="k",
    x="cast",
    yincrease=False,
    levels=pd.qcut(p06.rho.values.ravel(), 20, retbins=True)[1],
)
plt.gca().fill_between(p06.cast, 6100, -p06.bathy, color="k", zorder=10)
plt.gcf().set_size_inches(10, 10 / 1.6)

_images/a7e3ca6181de8436ac8dc6accc4f754e0399e86bd8b7294813bc262fadefbbc6.png

import sciviscolor as svc

svc.cm.

p06

%matplotlib inline

### p06['KT'].attrs['long_name'] = '$K_T$'
p06["KT"].attrs["units"] = "m$^2$/s"
p06["P"].attrs["long_name"] = "Pressure"
p06["P"].attrs["units"] = "dbar"
p06["lon"].attrs["long_name"] = "Longitude"
p06["lon"].attrs["units"] = "deg. East"

logKT = np.log10(p06["KT"])
logKT.attrs = p06.KT.attrs
logKT.attrs["long_name"] = "log$_{10}$ $K_T$"
logKT.plot(x="lon", y="P", cmap=mpl.cm.RdYlBu_r, robust=True)
(p06.sigma_4 - 1000).plot.contour(
    colors="k",
    x="lon",
    y="P",
    yincrease=False,
    labels=False,
    levels=pd.qcut(p06.sigma_4.values.ravel() - 1000, 15, retbins=True)[1],
)
plt.gca().fill_between(p06.lon, 6100, -p06.bathy, color="k", zorder=10)
plt.gcf().set_size_inches(10, 10 / 2.5)

plt.gca().set_title(
    "Potential density contours and log $K_T$ for the P06 GOSHIP section"
)

plt.savefig("../images/p06-kt.png", bbox_inches="tight")

/home/deepak/anaconda3/lib/python3.6/site-packages/matplotlib/contour.py:1000: UserWarning: The following kwargs were not used by contour: 'labels'
  s)

_images/cfbb118ae666d7a9a7d4ca1e130931b9b4dd4ca745e67e27034b0e941917bde6.png

bins = [
6,
6,
8,
0,
25,
5,
75,
0,
1,
2,
3,
4,
5,
6,
7,
8,
9,
]
dcpy.oceans.TSplot(p06.S, p06["T"], p06.pres, rho_levels=bins)

_images/df6485de19454766fcb56402a501543b26d208c9dfa03e97793969a312fa78e2.png

Calculate#

nbins = 18

p06["rho"] = p06["sigma_0"]
eccograd["ρmean"] = eccograd["sigma_0"]


eccoKe = ed.process_transect_1d(p06, eccograd, "ECCO", nbins=nbins)
ed.plot_bar_Ke(eccoKe, Ke_log=True)

argoKe = ed.process_transect_1d(p06, argograd, "ARGO", nbins=nbins)
ed.plot_bar_Ke(argoKe, Ke_log=True)

_images/e668c445150083e18d71abe9df542a16945eea9a411fb945c4373e8e7e7a9cd0.png

_images/03cc59fa59b6e097264ba5509b28ef930443f814c2f2eced40767f13b12b1bdd.png

Calculate 2#

import eddydiff as ed

p06 = xr.open_dataset("/home/deepak/work/eddydiff/datasets/P06/p06-new.nc")

##### get lat, lon, dist from older version
p06old = xr.open_dataset("/home/deepak/work/eddydiff/datasets/P06/p06.nc")
p06old = p06old.where(p06old.time.notnull(), drop=True)
p06 = p06.where(p06.time.notnull(), drop=True)
p06.update(
    p06old[["lon", "lat", "bathy"]]
    .swap_dims({"cast": "time"})
    .reindex(time=p06.time.values, method="nearest")
    .drop("cast")
    .swap_dims({"time": "cast"})
)
####

bathy = p06.bathy.copy()
p06 = p06.where(p06["KT"] < 1e-3)
p06["bathy"] = bathy
p06["KtTz"] = p06.KT * p06.dTdz
p06

if "eccograd" not in locals():
    eccograd, argograd, cole, aber = ed.read_all_datasets(kind="monthly", transect=p06)

# use sigma_0
p06["rho"] = p06["sigma_0"]
eccograd["ρmean"] = eccograd["sigma_0"]

bins = [
    1024.6,
    1025.6,
    1025.8,
    1026.0,
    1026.5,
    1026.9,
    1027.0,
    1027.1,
    1027.2,
    1027.3,
    1027.4,
    1027.5,
    1027.6,
    1027.7,
    1027.8,
    1027.9,
]

trdens, bins = ed.bin_to_density_space(p06, bins=bins)

eccograd["dist"] = p06.dist
eccodens, _ = ed.bin_to_density_space(eccograd.rename({"ρmean": "rho"}), bins)

newKe = xr.Dataset()
newKe["KT"] = trdens.KT.mean(dim="cast")
newKe["chi"] = trdens.chi.mean(dim="cast")
newKe["KtTz"] = trdens.KtTz.mean(dim="cast")
newKe["dTiso"] = np.abs(eccodens.dTiso.mean(dim="cast"))
newKe["dTmdz"] = np.abs(eccodens.dTdz.mean(dim="cast"))
newKe["dTdz"] = trdens.dTdz.mean(dim="cast")

newKe["Ke"] = (newKe["chi"] / 2 - np.abs(newKe["KtTz"] * newKe["dTmdz"])) / (
    newKe["dTiso"] ** 2
)
df = newKe.to_dataframe()
df.index = np.round(newKe.to_dataframe().index, 3)
ed.plot_bar_Ke(df)

_images/8094f0e4c74aba564c5c6fbf6ca99709f1f7574d666d00b982a327e378a416e9.png

proposal version#

bins = [
6,
6,
8,
0,
15,
3,
45,
6,
8,
0,
1,
2,
3,
4,
5,
6,
7,
8,
9,
]

ed.plot_bar_Ke(newKe.to_dataframe())

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

_images/895392d2b03c71c60c7335178101b51486787ccc4a9ef6a58b7ce620e8215c76.png

The above really awesome plot results after a few tweaks.

Doing things in σ_0.
Handcrafted density bins using the T-S diagram as guide. This is really important. Had to fiddle around slightly to not get negative values. Biggest improvement is that I can now resolve thermocline and coarsen the abyss. This is necessary because the big signal is in the thermocline and the abyssal values need to account for topography. Automatically choosing bins by pandas.cut/pandas.qcut would put too many bins down deep and fewer in the thermocline.
Throwing out all KT > 1e-3. This is really important. Will need to QC the values coming out of the CTD_chipod analysis code.
Redid (vectorized) the bin-averaging by cast. This is a better way to do it I think. Have to use pandas so I can groupby using multiple variables.
I tried fitting straight lines to the bin-averaged ECCO field to get dTiso but this seems to under-estimate values. Currently I bin-average dTiso and use that. This might not be crazy because it mirrors what I do with χ and we know that χ has to go along with the appropriate local gradients. i.e. if high χ coincides with high dTiso locally, we want the averaged dTiso to be biased high since averaged χ will be biased high.

f, ax = plt.subplots(2, 2, constrained_layout=True)

Trms = trdens["T"].std(dim="cast")
Trms.name = "RMS potential temp."
Trms.attrs["units"] = "°C"
Trms.plot(ax=ax[0, 0], y="rho", yincrease=False)
newKe.Ke.where(newKe.Ke > 0).plot(marker=".", ax=ax[0, 1], y="rho", yincrease=False)
ax[0, 0].set_ylabel("$σ_0$ [kg/m$^3$]")
ax[0, 1].set_ylabel("$σ_0$ [kg/m$^3$]")
ax[0, 1].set_xlabel("$K_e [m^2/s^2]$ ")

ax[1, 0].plot(newKe.Ke.where(newKe.Ke > 0), trdens.P.mean(dim="cast"), ".-")
ax[1, 0].plot(cole.diffusivity.mean(dim="cast"), cole.P)
ax[1, 0].legend(["χpod estimate", "Cole et al. (2015)"])
ax[1, 0].invert_yaxis()
ax[1, 0].set_ylabel("Pressure [dbar]")
ax[1, 0].set_xlabel("$K_e [m^2/s]$ ")

dcpy.oceans.TSplot(
    p06.S,
    p06["T"],
    rho_levels=bins,
    hexbin=False,
    ax=ax[1, 1],
    size=1,
    plot_distrib=False,
)

f.suptitle(p06.transect_name)
f.set_size_inches(8, 10)

plt.savefig("../images/p06.png", bbox_inches="tight")

/home/deepak/work/python/xarray/xarray/core/nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)

_images/6fd180c0d954ad2fd576f9803bc79a789e2a7966a78418a4f28e0065322a6864.png

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

zvec = trdens.P.mean(dim="cast")
rhovec = newKe.rho

ax[0].plot(newKe.Ke, zvec, ".-")
ax[0].plot(cole.diffusivity.mean(dim="cast"), cole.P)
ax[0].invert_yaxis()
ax[0].set_ylabel("Pressure [dbar]")
ax[0].set_xlabel("Transect-wide mean $K_e [m^2/s]$ ")

ax[0].spines["right"].set_visible(True)
ax[0].spines["top"].set_visible(True)

# https://matplotlib.org/gallery/subplots_axes_and_figures/fahrenheit_celsius_scales.html
ax_z = ax[0]
ax_rho = ax_z.twinx()
ax_rho.set_yticks(ax_z.get_yticks())
ax_rho.invert_yaxis()
ax_rho.set_yticklabels(
    [
        str(np.round(rr - 1000, 2))
        for rr in np.interp(ax_z.get_yticks(), zvec, rhovec.values)
    ]
)
ax_rho.set_ylabel("$σ_0$ [kg/m$^3$]")
ax_rho.set_ylim(ax_z.get_ylim())

ax[0].plot(
    [
        aber.mean(dim="cast").K_OC_LAT,
        aber.mean(dim="cast").K_OC_PSI,
        aber.mean(dim="cast").K_OC_SST,
    ],
    [0, 0, 0],
    "ro",
)
ax[0].legend(
    ["χpod estimate", "Cole et al. (2015)", "Abernathey & Marshall (2013)"],
    loc="lower right",
)

dcpy.oceans.TSplot(
    p06.S, p06.T, 0, rho_levels=bins, ax=ax[1], size=1, hexbin=False, plot_distrib=False
)

f.suptitle(p06.transect_name)
f.set_size_inches(8, 5)

# plt.savefig("../images/p06-Ke.png", bbox_inches="tight")

/home/deepak/work/python/xarray/xarray/core/nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)

_images/d5d127f75f138be54da8ee5d8e04b0c8de62229fd877d87a896f18aff6ef238c.png

newKe.Ke.where(newKe.Ke > 0).isel(
    cast=slice(240)
)  # .coarsen(cast=5).mean().plot(x="cast", norm=mpl.colors.LogNorm(0.1, 5e3), robust=True)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-259-9f6bb7b7d953> in <module>
----> 1 newKe.Ke.where(newKe.Ke > 0).isel(cast=slice(240)) #.coarsen(cast=5).mean().plot(x="cast", norm=mpl.colors.LogNorm(0.1, 5e3), robust=True)

~/work/python/xarray/xarray/core/dataarray.py in isel(self, indexers, drop, missing_dims, **indexers_kwargs)
   1051         # lists, or zero or one-dimensional np.ndarray's
   1052 
-> 1053         variable = self._variable.isel(indexers, missing_dims=missing_dims)
   1054 
   1055         coords = {}

~/work/python/xarray/xarray/core/variable.py in isel(self, indexers, missing_dims, **indexers_kwargs)
   1062 
   1063         key = tuple(indexers.get(dim, slice(None)) for dim in self.dims)
-> 1064         return self[key]
   1065 
   1066     def squeeze(self, dim=None):

~/work/python/xarray/xarray/core/variable.py in __getitem__(self, key)
    702         array `x.values` directly.
    703         """
--> 704         dims, indexer, new_order = self._broadcast_indexes(key)
    705         data = as_indexable(self._data)[indexer]
    706         if new_order:

~/work/python/xarray/xarray/core/variable.py in _broadcast_indexes(self, key)
    541 
    542         if all(isinstance(k, BASIC_INDEXING_TYPES) for k in key):
--> 543             return self._broadcast_indexes_basic(key)
    544 
    545         self._validate_indexers(key)

~/work/python/xarray/xarray/core/variable.py in _broadcast_indexes_basic(self, key)
    569             dim for k, dim in zip(key, self.dims) if not isinstance(k, integer_types)
    570         )
--> 571         return dims, BasicIndexer(key), None
    572 
    573     def _validate_indexers(self, key):

~/work/python/xarray/xarray/core/indexing.py in __init__(self, key)
    376             new_key.append(k)
    377 
--> 378         super().__init__(new_key)
    379 
    380 

~/work/python/xarray/xarray/core/indexing.py in __init__(self, key)
    329         if type(self) is ExplicitIndexer:
    330             raise TypeError("cannot instantiate base ExplicitIndexer objects")
--> 331         self._key = tuple(key)
    332 
    333     @property

TypeError: descriptor '_key' for 'ExplicitIndexer' objects doesn't apply to 'BasicIndexer' object

Linear fits to ECCO T along ρ surfaces in the along-transect direction doesn’t work. The field has too much curvature

# recalculate gradients from mean fields in density space
dTiso = xr.ones_like(eccodens.rho) * np.nan
dTiso.name = "dTiso"
dTdz = xr.ones_like(eccodens.rho) * np.nan
dTdz.name = "dTdz"

for idx, rr in enumerate(eccodens.rho):
    Tvec = eccodens.Tmean.sel(rho=rr)
    mask = np.isnan(Tvec)
    if len(Tvec[~mask]) > 0:
        slope, intercept, r_value, p_value, std_err = sp.stats.linregress(
            eccodens.dist[~mask] * 1000, Tvec[~mask]
        )
    dTiso[idx] = slope

dTdz.values = -np.gradient(eccodens.Tmean.mean(dim="cast"), eccodens.P.mean(dim="cast"))

f, ax = plt.subplots(1, 2)
np.abs(dTiso).plot.line(ax=ax[0], y="rho")
np.abs(eccodens.dTiso).mean(dim="cast").plot.line(ax=ax[0], y="rho", yincrease=False)
ax[0].set_xlim([0, 5e-6])
f.legend(
    [
        "recalculated from density-averaged field",
        "along-transect mean of local gradient",
    ]
)

dTdz.plot(ax=ax[1], y="rho", yincrease=False)
eccodens.dTdz.mean(dim="cast").plot.line(ax=ax[1], y="rho", yincrease=False)

[<matplotlib.lines.Line2D at 0x7f1727a84630>]

_images/87f5e401045cbd34804e9aa75110ba244b19f1cdac03d5c4bced1ecf6deaa5e2.png

nn = 6

Tvec = eccodens.Tmean.isel(rho=nn)
mask = np.isnan(Tvec)
plt.plot(eccodens.dist[~mask] * 1000, Tvec[~mask], ".")

plt.plot(
    eccodens.dist * 1000,
    Tvec.mean() + dTiso.isel(rho=nn) * (eccodens.dist - eccodens.dist.mean()) * 1000,
)

[<matplotlib.lines.Line2D at 0x7f1725c2a6d8>]

_images/43c16f2060f0e8a1bb12cac14164f5ca84eb20b4c4a2c2f94ecc6866f9f79e3e.png

argograd.dTiso.mean("time").plot(x="cast", robust=True, yincrease=False)

<matplotlib.collections.QuadMesh at 0x7fd75806b350>

_images/a8e496c7a3c48944155b6509d8d565325cf6543fc94d210c389a68eff6e441e5.png

eccograd.dTiso.plot(x="cast", yincrease=False, robust=True)

<matplotlib.collections.QuadMesh at 0x7fd75a8a8510>

_images/74ef3628c907bcb0c0b9e7290ce11c7f8d8bf3e034544acf42eb8e248897d599.png

clim = eccograd.copy()

clim = clim.to_dataframe().reset_index()

trmean, ρbins = ed.average_transect_1d(p06, nbins=15)

climrho = clim.groupby(pd.cut(clim.ρmean, ρbins, precision=1))

cole.diffusivity.plot(norm=mpl.colors.LogNorm(), x="cast", yincrease=False)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-55-05ef08ea103d> in <module>
----> 1 cole.diffusivity.plot(norm=mpl.colors.LogNorm(), x='cast', yincrease=False)

NameError: name 'mpl' is not defined

cole.diffusivity.mean(dim="cast").plot()

[<matplotlib.lines.Line2D at 0x7f1730d097b8>]

_images/de07399ec4854796827a1d9ee355905a5fe018bdc4b849697d1aaa061f9a38c2.png

clim = eccograd.sel(time=9)
transect = p06.sel(P=slice(0, 2000))

trmean, ρbins = ed.average_transect_1d(transect, nbins=20)
gradmean = ed.average_clim(clim, transect, ρbins)

Ke = ed.estimate_Ke(trmean, gradmean)
Ke

	chi	KtTz	dTdz	dTmdz	dTiso	KT	Ke
rho
(1024.8000000000002, 1026.7]	1.340903e-07	1.843108e-06	0.021105	0.021170	4.694987e-06	0.000118	1271.439129
(1026.7, 1027.5]	8.350073e-08	1.431120e-06	0.026548	0.023706	3.095114e-06	0.000058	816.732447
(1027.5, 1028.2]	6.809042e-08	1.229691e-06	0.024866	0.026552	2.792055e-06	0.000053	178.802085
(1028.2, 1028.9]	3.888472e-08	9.499299e-07	0.018603	0.019622	1.978210e-06	0.000053	205.196687
(1028.9, 1029.4]	1.696346e-08	6.028373e-07	0.011949	0.011681	1.093078e-06	0.000053	1205.030222
(1029.4, 1029.9]	1.025537e-08	4.112075e-07	0.008578	0.009867	4.934073e-07	0.000047	4396.180518
(1029.9, 1030.4]	6.098227e-09	3.541930e-07	0.006906	0.006893	3.572366e-07	0.000053	4761.760099
(1030.4, 1030.9]	5.922741e-09	3.396340e-07	0.006503	0.005814	3.952476e-07	0.000051	6315.357062
(1030.9, 1031.4]	4.479737e-09	3.283833e-07	0.006106	0.005641	3.961862e-07	0.000056	2468.021433
(1031.4, 1031.9]	4.335862e-09	3.358884e-07	0.005741	0.005445	3.654930e-07	0.000059	2536.931307
(1031.9, 1032.4]	3.778694e-09	3.191589e-07	0.005229	0.006691	3.008868e-07	0.000064	-2719.446279
(1032.4, 1032.9]	2.663402e-09	2.663522e-07	0.004435	0.005343	2.868961e-07	0.000064	-1111.907562
(1032.9, 1033.4]	2.050237e-09	2.458359e-07	0.003670	0.002568	1.781828e-07	0.000071	12403.643288
(1033.4, 1033.9]	1.563165e-09	2.117540e-07	0.003097	0.003163	1.434013e-07	0.000068	5439.135721
(1033.9, 1034.4]	1.084065e-09	1.862677e-07	0.002532	0.002644	1.283914e-07	0.000076	3006.651780
(1034.4, 1034.9]	8.051773e-10	1.646219e-07	0.002085	0.002116	1.193445e-07	0.000082	3804.070508
(1034.9, 1035.4]	5.643653e-10	1.386612e-07	0.001779	0.001680	1.091446e-07	0.000080	4130.684457
(1035.4, 1035.9]	3.849335e-10	1.173186e-07	0.001456	NaN	1.000824e-07	0.000083	NaN
(1035.9, 1036.4]	2.681242e-10	1.045344e-07	0.001217	NaN	NaN	0.000090	NaN
(1036.4, 1036.9]	2.172361e-10	9.460649e-08	0.001071	NaN	8.956461e-08	0.000094	NaN

def compare_transect_clim(transect, clim):

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

    for aa, (trvar, climvar) in enumerate(
        zip(["T", "S", "rho"], ["Tmean", "Smean", "ρmean"])
    ):
        _, levels = pd.qcut(transect[trvar].values.ravel(), 20, retbins=True)

        (
            clim[climvar]
            .sel(**ed.get_region_from_transect(p06))
            .isel(lat=1)
            .plot.contour(
                ax=ax[aa],
                levels=levels,
                x="lon",
                add_colorbar=True,
                yincrease=False,
                cmap=mpl.cm.RdYlBu,
            )
        )
        (
            transect[trvar].plot.contour(
                ax=ax[aa], levels=levels, colors="k", x="lon", y="P", yincrease=False
            )
        )

    trname = (
        transect.attrs["transect_name"]
        if "transect_name" in transect.attrs
        else "transect"
    )
    f.suptitle("Compare " + trname + " vs " + clim.dataset, y=0.9)
    plt.gcf().set_size_inches((10, 12))


compare_transect_clim(p06, eccograd.sel(time=8))
plt.gca().set_ylim((2400, 0))
#  compare_transect_clim(p06, argograd)
# plt.gca().set_ylim((2400, 0))

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-19-6a0f5cf774f5> in <module>()
     21 
     22 
---> 23 compare_transect_clim(p06, eccograd.sel(time=8))
     24 plt.gca().set_ylim((2400, 0))
     25 #  compare_transect_clim(p06, argograd)

~/work/python/xarray/xarray/core/dataset.py in sel(self, indexers, method, tolerance, drop, **indexers_kwargs)
   1507         indexers = either_dict_or_kwargs(indexers, indexers_kwargs, 'sel')
   1508         pos_indexers, new_indexes = remap_label_indexers(
-> 1509             self, indexers=indexers, method=method, tolerance=tolerance)
   1510         result = self.isel(indexers=pos_indexers, drop=drop)
   1511         return result._replace_indexes(new_indexes)

~/work/python/xarray/xarray/core/coordinates.py in remap_label_indexers(obj, indexers, method, tolerance, **indexers_kwargs)
    353 
    354     pos_indexers, new_indexes = indexing.remap_label_indexers(
--> 355         obj, v_indexers, method=method, tolerance=tolerance
    356     )
    357     # attach indexer's coordinate to pos_indexers

~/work/python/xarray/xarray/core/indexing.py in remap_label_indexers(data_obj, indexers, method, tolerance)
    235     new_indexes = {}
    236 
--> 237     dim_indexers = get_dim_indexers(data_obj, indexers)
    238     for dim, label in iteritems(dim_indexers):
    239         try:

~/work/python/xarray/xarray/core/indexing.py in get_dim_indexers(data_obj, indexers)
    203     if invalid:
    204         raise ValueError("dimensions or multi-index levels %r do not exist"
--> 205                          % invalid)
    206 
    207     level_indexers = defaultdict(dict)

ValueError: dimensions or multi-index levels ['time'] do not exist

P06 CTD χpod analysis

Contents

P06 CTD χpod analysis#

Read CTD χpod data and convert to netCDF#

Read transect + ancillary datasets#

plots#

Calculate#

Calculate 2#

proposal version#

Calculate 3#