Access 4DMED-SEA 3D Physical Fields

import xarray as xr 

zarr_href = 'https://s3.waw4-1.cloudferro.com/EarthCODE/OSCAssets/ocean_datasets/tsgeo150.zarr/'

ds = xr.open_zarr(zarr_href)
ds

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# -----------------------------
# Settings
# -----------------------------
date = "2020-06-15"

depths = [10, 50, 100]  

salinity = ds["so"].sel(
    time=date,
    method="nearest"
)

# Select nearest available depths
salinity_depths = [
    salinity.sel(depth=d, method="nearest")
    for d in depths
]

# Load only these three 2-D slices into memory
salinity_depths = [x.compute() for x in salinity_depths]

# Actual date selected
actual_date = np.datetime_as_string(
    salinity_depths[0].time.values,
    unit="D"
)

# Using the same scale makes the depths directly comparable
vmin = min(float(x.min(skipna=True)) for x in salinity_depths)
vmax = max(float(x.max(skipna=True)) for x in salinity_depths)

fig, axes = plt.subplots(
    1, 3,
    figsize=(18, 6),
    subplot_kw={"projection": ccrs.PlateCarree()},
    constrained_layout=True
)

for ax, data, requested_depth in zip(
    axes,
    salinity_depths,
    depths
):
    actual_depth = float(data.depth.values)

    pcm = ax.pcolormesh(
        ds["longitude"],
        ds["latitude"],
        data,
        transform=ccrs.PlateCarree(),
        shading="auto",
        cmap="viridis",
        vmin=vmin,
        vmax=vmax
    )

    # Basemap
    ax.add_feature(
        cfeature.LAND,
        facecolor="lightgray",
        zorder=10
    )

    ax.coastlines(
        resolution="10m",
        linewidth=0.8,
        zorder=11
    )

    ax.add_feature(
        cfeature.BORDERS,
        linewidth=0.4,
        zorder=11
    )

    # Mediterranean extent from dataset
    ax.set_extent(
        [
            float(ds.longitude.min()),
            float(ds.longitude.max()),
            float(ds.latitude.min()),
            float(ds.latitude.max())
        ],
        crs=ccrs.PlateCarree()
    )
    
    ax.set_title(
        f"Depth = {abs(actual_depth):.0f} m"
    )

cbar = fig.colorbar(
    pcm,
    ax=axes,
    orientation="horizontal",
    fraction=0.06,
    pad=0.08
)

cbar.set_label(
    r"sea water salinity (practical salinity unit)"
)

fig.suptitle(
    f"Sea water salinity - {actual_date}",
    fontsize=16
)

plt.show()