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
from xoa.plot import plot_flow

xr.set_options(display_style="text");

[1]:

<xarray.core.options.set_options at 0x7f378908b520>

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(f.read())

[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¶

Hycom velocities¶

U and V are stored in separate files at U and V locations on the Hycom Arakawa C-Grid. We choose to rename the dataset contents so that we can use generic names to access it. We call for that the decode method the xoa accessor.

[5]:

# U velocity component
u_hycom = xoa.open_data_sample("hycom.gdp.u.nc").xoa.decode().u.squeeze(drop=True)
# V velocity component
v_hycom = xoa.open_data_sample("hycom.gdp.v.nc").xoa.decode().v.squeeze(drop=True)
# Layer thicknesses
hycom = xoa.open_data_sample("hycom.gdp.h.nc").xoa.decode().squeeze(drop=True)

We now compute the depths from the layer thicknesses at the T location thanks to the xoa.grid.dz2depth function: we make integrate from a null SSH and interpolate the depth from W locations to T locations.

[6]:

hycom.coords["depth"] = dz2depth(hycom.dz, centered=True)

Now we interpolate the velocity components to the T location with the xoa.grid.shift function.

[7]:

ut_hycom = shift(u_hycom, {"x": "left", "y": "left"})
vt_hycom = shift(v_hycom, {"x": "left", "y": "left"})
hycom["u"] = ut_hycom.assign_coords(**hycom.coords)
hycom["v"] = vt_hycom.assign_coords(**hycom.coords)

So, finally we obtain:

[8]:

hycom

[8]:

<xarray.Dataset>
Dimensions:  (time: 50, x: 13, y: 11, z: 6)
Coordinates:
    lon      (y, x) float32 -3.2 -3.125 -3.05 -2.975 ... -2.45 -2.375 -2.3
    lat      (y, x) float32 44.69 44.69 44.69 44.69 ... 45.22 45.22 45.22 45.22
  * z        (z) int32 1 2 3 4 5 6
  * time     (time) datetime64[ns] 2021-03-01 ... 2021-03-25T12:00:00
    depth    (time, z, y, x) float32 1.511 1.511 1.51 ... 20.11 20.09 20.06
Dimensions without coordinates: x, y
Data variables:
    bathy    (y, x) float32 ...
    dz       (time, z, y, x) float32 3.021 3.021 3.02 3.019 ... 5.22 5.216 5.209
    u        (time, z, y, x) float32 0.03078 0.02984 0.03803 ... 0.1292 0.1364
    v        (time, z, y, x) float32 0.0497 0.07527 0.1111 ... 0.0498 0.05054
Attributes: (12/26)
    history:                   Fri Apr  9 12:15:22 2021: ncrcat -d X,472,508,...
    Conventions:               COARDS
    nhybrd:                    20
    nsigma:                    19
    kapref:                    2
    hybflg:                    0
    ...                        ...
    dp00ft:                    [1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 ...
    dp00xt:                    [21.874 21.414 21.414 13.197  7.85   3.85   4....
    ds00t:                     [3.555 3.026 3.026 1.547 0.828 0.828 0.968 0.8...
    ds00ft:                    [1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 ...
    ds00xt:                    [21.874 21.414 21.414  7.197  3.85   3.85   4....
    nco_openmp_thread_number:  1

GDP drifter¶

The drifter comes as a csv file and we read it as pandas.DataFrame instance.

[9]:

drifter = xoa.open_data_sample("gdp-6203641.csv", header=0, skiprows=[1], parse_dates=[2], index_col=2)

Since the sampling is not that nice, we resample it to 3-hour intervals.

[10]:

drifter = drifter.resample("3H").mean()

We convert it to and xarray.Dataset, fix time and decode it.

[11]:

drifter = drifter.to_xarray().assign_coords(time=drifter.index.values).xoa.decode()

We drop missing values.

[12]:

drifter = drifter.where(~drifter.lon.isnull() & ~drifter.lat.isnull(), drop=True)

We add a constant depth of 15 m.

[13]:

drifter.coords["depth"] = drifter.lon*0 + 15

Here is what we obtain.

[14]:

drifter

[14]:

<xarray.Dataset>
Dimensions:        (time: 159)
Coordinates:
  * time           (time) datetime64[ns] 2021-03-03T18:00:00 ... 2021-03-24T0...
    lat            (time) float64 44.68 44.72 44.72 44.72 ... 44.72 44.71 44.71
    lon            (time) float64 -2.544 -2.583 -2.592 ... -2.523 -2.514 -2.505
    depth          (time) float64 15.0 15.0 15.0 15.0 ... 15.0 15.0 15.0 15.0
Data variables:
    platform_code  (time) float64 6.204e+06 6.204e+06 ... 6.204e+06 6.204e+06
    sst            (time) float64 13.09 12.79 12.77 12.78 ... 12.56 12.53 12.61
    slp            (time) float64 1.027e+03 1.024e+03 ... 1.025e+03 1.026e+03
    lon360         (time) float64 357.5 357.4 357.4 357.4 ... 357.5 357.5 357.5

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. The background currents are plotted with the xoa.plot.plot_flow 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)
plot_flow(uh15, vh15, axes=ax, transform=pcarr, color="w", alpha=(.3, 1), linewidth=.6)
qv = 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 0x7f377ee912e0>

The discrepancies between the lagrangian and mean eulerian currents highlight the variability due to the mesocale activity. The differences between the observed and modeled currents are partly due to the lack of data assimilation in the model.