Access WOC Total Surface Current at 15 m

* Original Data Source: https://www.worldoceancirculation.org/Products#/metadata/b72f4ac1-ffc6-4ba9-a729-670b6e879d52
* Reference:  https://doi.org/10.25423/cmcc/4dmedsea_biophys_rep_3d
* OSC entry: https://opensciencedata.esa.int/products/total-surface-current-15/collection
* License: CC-BY-4.0

import xarray as xr 

zarr_href = 'https://s3.waw4-1.cloudferro.com/EarthCODE/OSCAssets/ocean_datasets/woc.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

# Choose timestamp
timestamp = "2016-06-1T15:00:00"

# Select nearest available time and remove singleton depth dimension
ds_t = ds.sel(time=timestamp, method="nearest").squeeze("depth")

# Vector components
u = ds_t["ut_woc"]
v = ds_t["vt_woc"]

# Current speed
speed = np.hypot(u, v)


fig, ax = plt.subplots(
    figsize=(12, 7),
    subplot_kw={"projection": ccrs.PlateCarree()}
)

# Plot speed as background
pcm = ax.pcolormesh(
    ds_t.lon,
    ds_t.lat,
    speed,
    cmap="viridis",
    shading="auto",
    transform=ccrs.PlateCarree()
)

# Subsample vectors so arrows do not overlap
step = 5

q = ax.quiver(
    ds_t.lon.values[::step],
    ds_t.lat.values[::step],
    u.values[::step, ::step],
    v.values[::step, ::step],
    transform=ccrs.PlateCarree(),
    scale=10,
    width=0.002,
    color="black"
)

# Reference arrow
ax.quiverkey(
    q,
    X=0.9,
    Y=1.03,
    U=0.5,
    label="0.5 m s$^{-1}$",
    labelpos="E"
)

ax.coastlines(resolution="10m")
ax.add_feature(cfeature.LAND, facecolor="lightgray")
ax.add_feature(cfeature.BORDERS, linewidth=0.5)


# Colorbar
cbar = fig.colorbar(
    pcm,
    ax=ax,
    orientation="vertical",
    pad=0.02
)
cbar.set_label("Current speed [m s$^{-1}$]")

actual_time = np.datetime_as_string(ds_t.time.values, unit="h")

ax.set_title(
    f"Total ocean current (CMEMS)\n{actual_time}"
)

plt.show()