Compare Hycom3d with a GDP drifter¶

In this notebook, we show :

how to decode dataset so that it is easy to access generic coordinates and variables,
how to compute depths from layer thicknesses,
how to interpolate currents from U and V positions to T position on an arakawa C grid,
how to perform a 4D interpolation with a variable depth coordinate to random positions,
how to make an horizontal slice of a 4D variable with a variable depth coordinate.

Initialisations¶

[1]:

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xoa
from xoa.grid import dz2depth, shift
from xoa.regrid import grid2loc, regrid1d
import xoa.cf as xcf
import xoa.geo as xgeo

Register the xoa xarray accessor :

[2]:

xoa.register_accessors()

Register the Hycom naming specifications:

[3]:

xcf.register_cf_specs(xoa.get_data_sample("hycom.cfg"))

Here is what these CF specifications contain:

[4]:

with open(xoa.get_data_sample("hycom.cfg")) as f:
    print(''.join([line for line in f if len(line) > 1]))

[register]
name=hycom
    [[attrs]]
    nsigma = "[0-9]"
[vertical]
positive=down
type=dz
[data_vars]
    [[bathy]]
    name=bathymetry
    [[dz]]
    name=h
[coords]
    [[x]]
    name=X
    [[y]]
    name=Y
    [[z]]
    name=lev

Note that most of these specifications are not needed to decode the dataset, i.e to find known variables and coordinates, since the default specifications help for that.

Read and decode the datasets¶

Compute and interpolate velocities¶

We compute the drifter velocity components

[15]:

drifter["u"] = drifter.lon.differentiate("time", datetime_unit="s") * xgeo.EARTH_RADIUS*np.pi/180
drifter["u"] *= np.cos(np.radians(drifter.lat.values))
drifter["v"] = drifter.lat.differentiate("time", datetime_unit="s") * xgeo.EARTH_RADIUS*np.pi/180

We interpolate the Hycom velocity to the drifter positions with xoa.regrid.grid2loc. note that this is a 4D linear interpolation.

[16]:

uloc = grid2loc(hycom["u"], drifter)
vloc = grid2loc(hycom["v"], drifter)

Instead of just showing the velocities along the drifter positions, we can plot the mean model velocities over the same period as a background. So we interpolate them at 15 m with xoa.regrid.regrid1d and compute a time average.

[17]:

d15 = xr.DataArray([15.], dims="depth", name="depth")
uh15 = regrid1d(hycom["u"], d15).squeeze(drop=True).mean(dim="time")
vh15 = regrid1d(hycom["v"], d15).squeeze(drop=True).mean(dim="time")

Plot¶

The geographic extent is easily computed with the xoa.geo.get_extent function.

[18]:

pmerc = ccrs.Mercator()
pcarr = ccrs.PlateCarree()
fig, ax = plt.subplots(figsize=(7, 7), subplot_kw={"facecolor": "teal", "projection": pmerc})
ax.gridlines(draw_labels=True, dms=True)
ax.set_extent(xgeo.get_extent(uh15))
kwqv = dict(scale_units="dots", scale=0.1/20, units="dots", transform=pcarr)
qv = ax.quiver(uh15.lon.values, uh15.lat.values, uh15.values, vh15.values,
          width=3, color="w", alpha=0.3, **kwqv)
ax.quiver(uloc.lon.values, uloc.lat.values, uloc.values, vloc.values,
          color="w", width=2, label="Model", **kwqv)
ax.plot(drifter.lon.values, drifter.lat.values, '-', color="C1", transform=pcarr, lw=.5)
ax.quiver(drifter.lon.values, drifter.lat.values, drifter.u.values, drifter.v.values,
          color="C1", label="Drifter", width=2, **kwqv)
plt.quiverkey(qv, 0.1, 1.06, 0.1, r"0.1 $m\,s^{-1}$", color="k", alpha=1, labelpos="E")
plt.legend();

[18]:

<matplotlib.legend.Legend at 0x7f1bf4926b20>

Compare Hycom3d with a GDP drifter¶

Initialisations¶

Read and decode the datasets¶

Hycom velocities¶

GDP drifter¶

Compute and interpolate velocities¶

Plot¶