Regridding

Regridding#

HYCOM output lives on two non-standard grids:

  • Horizontal: curvilinear, with 2-D lon/lat and velocities on a staggered, rotated C-grid.

  • Vertical: hybrid isopycnal, with layer thicknesses stored as pressure (thknss).

xhycom.regrid maps both onto a regular lon/lat/depth grid, conservatively by default (area-conservative horizontally, depth-integral-conserving vertically), for comparison with reanalyses like GLORYS.

import numpy as np
import xarray as xr
import xhycom
import matplotlib.pyplot as plt

grid = "/cluster/home/nlo043/NERSC-HYCOM-CICE/TP2a0.10/topo/regional.grid"
DATA_PATH = "/nird/datalake/NS9481K/shuang/TP2_output/expt_02.8/"

Open HYCOM output#

Load the 2020 archive year with its grid. postprocess=True converts thicknesses and SSH to metres, and chunks={"time": 1} loads the data lazily, one time step at a time. Learn more about choosing chunk sizes and working with lazy-loaded data in Lazy loading & chunking.

ds = xhycom.open_mfdataset(DATA_PATH + "archm.2020*", grid=grid, chunks={"time": 1}, postprocess=True)

ds

<xarray.Dataset> Size: 1TB
Dimensions:     (y: 380, x: 400, time: 366, k: 50, ki: 51)
Coordinates:
    lon         (y, x) float64 1MB dask.array<chunksize=(380, 400), meta=np.ndarray>
    lat         (y, x) float64 1MB dask.array<chunksize=(380, 400), meta=np.ndarray>
    lon_u       (y, x) float64 1MB dask.array<chunksize=(380, 400), meta=np.ndarray>
    lat_u       (y, x) float64 1MB dask.array<chunksize=(380, 400), meta=np.ndarray>
    lon_v       (y, x) float64 1MB dask.array<chunksize=(380, 400), meta=np.ndarray>
    lat_v       (y, x) float64 1MB dask.array<chunksize=(380, 400), meta=np.ndarray>
  * time        (time) object 3kB 2020-01-01 00:00:00 ... 2020-12-31 00:00:00
  * k           (k) int64 400B 1 2 3 4 5 6 7 8 9 ... 42 43 44 45 46 47 48 49 50
    dens        (k) float64 400B dask.array<chunksize=(50,), meta=np.ndarray>
  * ki          (ki) int64 408B 0 1 2 3 4 5 6 7 8 ... 42 43 44 45 46 47 48 49 50
Dimensions without coordinates: y, x
Data variables: (12/83)
    montg1      (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    srfhgt      (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    oneta       (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    surflx      (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    wtrflx      (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    salflx      (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    ...          ...
    CO2_wind    (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    ECO_bots    (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    surface__1  (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    surface__2  (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    si_u        (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    si_v        (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
Attributes:
    iversn:        23
    iexpt:         28
    yrflag:        3
    archive_type:  mean

Surface temperature on the native curvilinear grid:

ds["temp"].isel(time=0, k=0).plot(cmap="inferno", center=False)

<matplotlib.collections.QuadMesh at 0x15047bb930e0>
_images/f56345cf173eb9fdb5680fb7869b79556040677217e1b7eca76e71b07d25d198.png

Regrid onto a GLORYS grid#

GLORYS ships a regular lon/lat/depth grid with a land mask. Pass it straight to regrid as target=, lon/lat/depth and the mask come from it.

GRIDS = "/nird/datapeak/NS9481K/MERCATOR_DATA/REGULAR_GRID_COORD"
glorys = xr.open_dataset(f"{GRIDS}/GLO-MFC_001_030_mask_bathy.nc")
glorys

<xarray.Dataset> Size: 511MB
Dimensions:     (depth: 50, latitude: 2041, longitude: 4320)
Coordinates:
  * depth       (depth) float32 200B 0.494 1.541 2.646 ... 5.275e+03 5.728e+03
  * latitude    (latitude) float32 8kB -80.0 -79.92 -79.83 ... 89.83 89.92 90.0
  * longitude   (longitude) float32 17kB -180.0 -179.9 -179.8 ... 179.8 179.9
Data variables:
    mask        (depth, latitude, longitude) int8 441MB ...
    deptho      (latitude, longitude) float32 35MB ...
    deptho_lev  (latitude, longitude) float32 35MB ...
Attributes:
    title:                Bathymetry and mask for product GLOBAL_MULTIYEAR_PH...
    references:           http://marine.copernicus.eu
    credit:               E.U. Copernicus Marine Service Information (CMEMS)
    licence:              http://marine.copernicus.eu/services-portfolio/serv...
    contact:              servicedesk.cmems@mercator-ocean.eu
    producer:             CMEMS - Global Monitoring and Forecasting Centre
    institution:          Mercator Ocean
    Conventions:          CF-1.6
    area:                 GLOBAL
    source:               GLORYS12V1
    product_user_manual:  http://marine.copernicus.eu/documents/PUM/CMEMS-GLO...
# Conservative horizontal + vertical, in one call, onto the GLORYS grid.
%time ds_g = xhycom.regrid(ds, target=glorys, grid=grid)

CPU times: user 54.7 s, sys: 1.91 s, total: 56.6 s
Wall time: 2min 1s

This operation was slow: it computes the xESMF weight matrix that maps the HYCOM grid onto the GLORYS grid. The regridding itself is lazy, regrid returns dask-backed fields and the actual remap is deferred until you load them, so the time above is essentially all weight generation. Those weights depend only on the source and target grids, not on the field values or the time step, so the cost is the same whether you regrid a single slice or a whole year. That also means they can be built once and reused, which is done in the Caching the regrid weights section below.

ds_g # repr shows dask chunks, nothing computed yet

<xarray.Dataset> Size: 52TB
Dimensions:     (time: 366, lat: 624, lon: 4320, depth: 50, ki: 51)
Coordinates:
  * time        (time) object 3kB 2020-01-01 00:00:00 ... 2020-12-31 00:00:00
  * lat         (lat) float32 2kB 38.08 38.17 38.25 38.33 ... 89.83 89.92 90.0
  * lon         (lon) float32 17kB -180.0 -179.9 -179.8 ... 179.8 179.8 179.9
  * depth       (depth) float64 400B 0.494 1.541 2.646 ... 5.275e+03 5.728e+03
  * ki          (ki) int64 408B 0 1 2 3 4 5 6 7 8 ... 42 43 44 45 46 47 48 49 50
Data variables: (12/82)
    montg1      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    srfhgt      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    oneta       (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    surflx      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    wtrflx      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    salflx      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    ...          ...
    light_sw    (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    attenuat    (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    total_ch    (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    total_ca    (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    u-vel.      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    v-vel.      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
Attributes:
    iversn:         23
    iexpt:          28
    yrflag:         3
    archive_type:   mean
    regrid_method:  conservative
%time ds_g.temp.isel(time=0, depth=0).plot(figsize=(10, 4), cmap="inferno", center=False)

CPU times: user 8.36 s, sys: 1.24 s, total: 9.6 s
Wall time: 9.72 s

<matplotlib.collections.QuadMesh at 0x15046a685e50>
_images/bd7d8a07e45b13e1895111ff10d8cd3352169c39dc181edfb9bede3a893ce45b.png 
%time ds_g.temp.isel(time=0).sel(depth=200, method="nearest").plot(figsize=(10, 4), cmap="inferno", center=False)

CPU times: user 8.37 s, sys: 1.09 s, total: 9.46 s
Wall time: 9.53 s

<matplotlib.collections.QuadMesh at 0x15046953d1d0>
_images/3a8861b24689584c3aba87264a2c41ad718202e457c76f119664f4bad34a5a97.png

Comparing to GLORYS#

With HYCOM on the GLORYS grid, a like-for-like comparison against the reanalysis is a straightforward subtraction, open a GLORYS field for the matching date and difference it.

# GLORYS reanalysis fields for the same date (thetao = potential temperature).
# Adjust the path / variable names to your GLORYS product if they differ.
glorys_data = xr.open_dataset(
    "/nird/datapeak/NS9481K/MERCATOR_DATA/PHY/2020/MERCATOR-PHY-24-2020-01-01-12.nc"
)

hycom_temp = ds_g["temp"].isel(time=0).sel(depth=250, method="nearest").compute()
glorys_temp = glorys_data["thetao"].isel(time=0).sel(depth=250, method="nearest").rename(longitude="lon", latitude="lat").compute()

vmin = hycom_temp.min().values
vmax = hycom_temp.max().values

fig, axs = plt.subplots(3, 1, figsize=(9, 11))

hycom_temp.plot(cmap="inferno", vmin=vmin, vmax=vmax, ax=axs[0])
glorys_temp.plot(cmap="inferno", vmin=vmin, vmax=vmax, ax=axs[1])
(hycom_temp - glorys_temp).plot()

for ax in axs:
    ax.set_ylim(50, None)      # latitude > 50

axs[0].set_title("HYCOM")
axs[1].set_title("GLORYS")
axs[2].set_title("HYCOM - GLORYS")
fig.suptitle(f"Temperature at {glorys_temp.depth.values.round()}m depth on {str(glorys_temp.time.values)[:10]}", y=1.01)
plt.tight_layout()
_images/5b33b2f2f24d1995ad27c3723f9830c358a3b8796d80f444b4e35e5d1308ff7f.png

Custom target grid (arrays)#

Instead of a target= grid you can pass lon/lat/depth arrays directly:

lon = np.arange(-40, 40, 0.1)
lat = np.arange(55, 82, 0.1)
depth = np.array([0, 10, 20, 50, 100, 200, 500, 1000.0])
ds_box = xhycom.regrid(ds, lon=lon, lat=lat, depth=depth, grid=grid)

ds_box

<xarray.Dataset> Size: 305GB
Dimensions:     (time: 366, lat: 270, lon: 800, ki: 51, depth: 8)
Coordinates:
  * time        (time) object 3kB 2020-01-01 00:00:00 ... 2020-12-31 00:00:00
  * lat         (lat) float64 2kB 55.0 55.1 55.2 55.3 ... 81.6 81.7 81.8 81.9
  * lon         (lon) float64 6kB -40.0 -39.9 -39.8 -39.7 ... 39.7 39.8 39.9
  * ki          (ki) int64 408B 0 1 2 3 4 5 6 7 8 ... 42 43 44 45 46 47 48 49 50
  * depth       (depth) float64 64B 0.0 10.0 20.0 50.0 100.0 200.0 500.0 1e+03
Data variables: (12/82)
    montg1      (time, lat, lon) float64 632MB dask.array<chunksize=(1, 270, 800), meta=np.ndarray>
    srfhgt      (time, lat, lon) float64 632MB dask.array<chunksize=(1, 270, 800), meta=np.ndarray>
    oneta       (time, lat, lon) float64 632MB dask.array<chunksize=(1, 270, 800), meta=np.ndarray>
    surflx      (time, lat, lon) float64 632MB dask.array<chunksize=(1, 270, 800), meta=np.ndarray>
    wtrflx      (time, lat, lon) float64 632MB dask.array<chunksize=(1, 270, 800), meta=np.ndarray>
    salflx      (time, lat, lon) float64 632MB dask.array<chunksize=(1, 270, 800), meta=np.ndarray>
    ...          ...
    light_sw    (time, lat, lon, depth) float64 5GB dask.array<chunksize=(1, 270, 800, 8), meta=np.ndarray>
    attenuat    (time, lat, lon, depth) float64 5GB dask.array<chunksize=(1, 270, 800, 8), meta=np.ndarray>
    total_ch    (time, lat, lon, depth) float64 5GB dask.array<chunksize=(1, 270, 800, 8), meta=np.ndarray>
    total_ca    (time, lat, lon, depth) float64 5GB dask.array<chunksize=(1, 270, 800, 8), meta=np.ndarray>
    u-vel.      (time, lat, lon, depth) float64 5GB dask.array<chunksize=(1, 270, 800, 8), meta=np.ndarray>
    v-vel.      (time, lat, lon, depth) float64 5GB dask.array<chunksize=(1, 270, 800, 8), meta=np.ndarray>
Attributes:
    iversn:         23
    iexpt:          28
    yrflag:         3
    archive_type:   mean
    regrid_method:  conservative
ds_box.salin.isel(time=10).sel(depth=20).plot()

<matplotlib.collections.QuadMesh at 0x1503a0025bd0>
_images/89ac65cf0b56a4005703065a5d2449eb4476934b32cb8c0c9302a960421866d2.png

Run the two regridding steps separately#

regrid is just regrid_horizontal followed by regrid_vertical, and each is useful on its own.

To change only the vertical: put HYCOM’s hybrid layers onto fixed depth levels while keeping the native curvilinear horizontal grid, call regrid_vertical by itself. It needs no grid=, only thknss:

# Vertical only: hybrid layers -> z-levels; native curvilinear lon/lat (y, x) kept.
ds_v = xhycom.regrid_vertical(ds, depth=glorys["depth"].values)
ds_v

<xarray.Dataset> Size: 1TB
Dimensions:     (y: 380, x: 400, time: 366, ki: 51, depth: 50)
Coordinates:
    lon         (y, x) float64 1MB dask.array<chunksize=(380, 400), meta=np.ndarray>
    lat         (y, x) float64 1MB dask.array<chunksize=(380, 400), meta=np.ndarray>
    lon_u       (y, x) float64 1MB dask.array<chunksize=(380, 400), meta=np.ndarray>
    lat_u       (y, x) float64 1MB dask.array<chunksize=(380, 400), meta=np.ndarray>
    lon_v       (y, x) float64 1MB dask.array<chunksize=(380, 400), meta=np.ndarray>
    lat_v       (y, x) float64 1MB dask.array<chunksize=(380, 400), meta=np.ndarray>
  * time        (time) object 3kB 2020-01-01 00:00:00 ... 2020-12-31 00:00:00
  * ki          (ki) int64 408B 0 1 2 3 4 5 6 7 8 ... 42 43 44 45 46 47 48 49 50
  * depth       (depth) float64 400B 0.494 1.541 2.646 ... 5.275e+03 5.728e+03
Dimensions without coordinates: y, x
Data variables: (12/82)
    montg1      (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    srfhgt      (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    oneta       (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    surflx      (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    wtrflx      (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    salflx      (time, y, x) float64 445MB dask.array<chunksize=(1, 380, 400), meta=np.ndarray>
    ...          ...
    ECO_c2ch_2  (time, y, x, depth) float64 22GB dask.array<chunksize=(1, 380, 400, 50), meta=np.ndarray>
    ECO_c2ch_3  (time, y, x, depth) float64 22GB dask.array<chunksize=(1, 380, 400, 50), meta=np.ndarray>
    light_sw    (time, y, x, depth) float64 22GB dask.array<chunksize=(1, 380, 400, 50), meta=np.ndarray>
    attenuat    (time, y, x, depth) float64 22GB dask.array<chunksize=(1, 380, 400, 50), meta=np.ndarray>
    total_ch    (time, y, x, depth) float64 22GB dask.array<chunksize=(1, 380, 400, 50), meta=np.ndarray>
    total_ca    (time, y, x, depth) float64 22GB dask.array<chunksize=(1, 380, 400, 50), meta=np.ndarray>
Attributes:
    iversn:        23
    iexpt:         28
    yrflag:        3
    archive_type:  mean
ds_v.salin.isel(time=0).sel(depth=50, method="nearest").plot()

<matplotlib.collections.QuadMesh at 0x1503d739b390>
_images/2e918e04f1037d2298e2786fada20bb43f881648cf78cae6103c77dd2af955b1.png

Or change only the horizontal: map the curvilinear grid to regular lon/lat while keeping the hybrid layers (k). This step needs grid= (it supplies pang to rotate velocities and qlon/qlat for conservative cell corners):

ds_h = xhycom.regrid_horizontal(ds, target=glorys, grid=grid)   # curvilinear -> regular lon/lat
ds_h

<xarray.Dataset> Size: 21TB
Dimensions:     (time: 366, lat: 624, lon: 4320, k: 50, ki: 51, y_b: 381,
                 x_b: 401)
Coordinates:
  * time        (time) object 3kB 2020-01-01 00:00:00 ... 2020-12-31 00:00:00
  * lat         (lat) float32 2kB 38.08 38.17 38.25 38.33 ... 89.83 89.92 90.0
  * lon         (lon) float32 17kB -180.0 -179.9 -179.8 ... 179.8 179.8 179.9
  * k           (k) int64 400B 1 2 3 4 5 6 7 8 9 ... 42 43 44 45 46 47 48 49 50
    dens        (k) float64 400B dask.array<chunksize=(50,), meta=np.ndarray>
  * ki          (ki) int64 408B 0 1 2 3 4 5 6 7 8 ... 42 43 44 45 46 47 48 49 50
    lon_b       (y_b, x_b) float64 1MB -94.76 -94.6 -94.43 ... 99.11 98.97 98.84
    lat_b       (y_b, x_b) float64 1MB 38.94 39.04 39.14 ... 56.3 56.21 56.12
Dimensions without coordinates: y_b, x_b
Data variables: (12/83)
    montg1      (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    srfhgt      (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    oneta       (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    surflx      (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    wtrflx      (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    salflx      (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    ...          ...
    u_btrop     (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    v_btrop     (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    umix        (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    vmix        (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    si_u        (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    si_v        (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
Attributes:
    iversn:         23
    iexpt:          28
    yrflag:         3
    archive_type:   mean
    regrid_method:  conservative
ds_h.temp.isel(time=0, k=0).plot(figsize=(10, 4), cmap="inferno", center=False)

<matplotlib.collections.QuadMesh at 0x1503d5a0e350>
_images/8646763ab6b6f01102f6426685257a7e30fa5e08beacb8970bd10279b2fe3de9.png

Chaining the two standalone pieces, horizontal then vertical, is exactly what regrid does by default. (Whether to run in this order or the reverse is a real choice, compared in the next section.)

# Chaining the standalone pieces == xhycom.regrid(ds, target=glorys, grid=grid):
ds_hv = xhycom.regrid_vertical(ds_h, depth=glorys["depth"].values)
ds_hv

<xarray.Dataset> Size: 20TB
Dimensions:     (time: 366, lat: 624, lon: 4320, ki: 51, depth: 50)
Coordinates:
  * time        (time) object 3kB 2020-01-01 00:00:00 ... 2020-12-31 00:00:00
  * lat         (lat) float32 2kB 38.08 38.17 38.25 38.33 ... 89.83 89.92 90.0
  * lon         (lon) float32 17kB -180.0 -179.9 -179.8 ... 179.8 179.8 179.9
  * ki          (ki) int64 408B 0 1 2 3 4 5 6 7 8 ... 42 43 44 45 46 47 48 49 50
  * depth       (depth) float64 400B 0.494 1.541 2.646 ... 5.275e+03 5.728e+03
Data variables: (12/82)
    montg1      (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    srfhgt      (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    oneta       (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    surflx      (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    wtrflx      (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    salflx      (time, lat, lon) float64 8GB dask.array<chunksize=(1, 624, 4320), meta=np.ndarray>
    ...          ...
    light_sw    (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    attenuat    (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    total_ch    (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    total_ca    (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    u-vel.      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    v-vel.      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
Attributes:
    iversn:         23
    iexpt:          28
    yrflag:         3
    archive_type:   mean
    regrid_method:  conservative

Which regridding order? Mix along isopycnals#

regrid runs the lateral and vertical steps in sequence, and the order is not neutral. HYCOM’s hybrid layers are isopycnals in the stratified interior, and temperature and salinity mix along isopycnals, so where the lateral averaging happens matters:

  • order="horizontal_first" (default) blends along the native layers → mixes same-density water, preserving water masses and the T–S relationship.

  • order="vertical_first" blends at constant depth → where isopycnals tilt (fronts, slopes) it mixes across density surfaces and smears water masses, though it resolves more cells on shelves and slopes.

Both conserve the global integral; they differ locally. Run both and compare:

g_hf = xhycom.regrid(ds, target=glorys, grid=grid, order="horizontal_first")
g_vf = xhycom.regrid(ds, target=glorys, grid=grid, order="vertical_first")

# extra ocean cells the depth-first order resolves near topography:
int(np.isfinite(g_vf["temp"]).isel(time=0).sum() - np.isfinite(g_hf["temp"]).isel(time=0).sum())

128308

The two agree across the smooth interior but diverge along fronts and topography. A meridional section of the difference shows constant-depth averaging mixing across tilted isopycnals, the dipoles sit on sloping density surfaces:

temp_vf = g_vf["temp"].isel(time=0).sel(lon=0, method="nearest").compute()
temp_hf = g_hf["temp"].isel(time=0).sel(lon=0, method="nearest").compute()

fig, axs = plt.subplots(3, 1, figsize=(8, 12))

temp_vf.plot(cmap="inferno", x="lat", y="depth", center=False, yincrease=False, ax=axs[0])
temp_hf.plot(cmap="inferno", x="lat", y="depth", center=False, yincrease=False, ax=axs[1])
(temp_vf - temp_hf).plot(x="lat", y="depth", vmin=-0.3, vmax=0.3, cmap="RdBu_r", center=0, yincrease=False, ax=axs[2])

for ax in axs:
    ax.set_xlim(50, None)      # latitude > 50
    ax.set_ylim(4200, None)    # depth < 4200 (see note below)

axs[0].set_title("Horizontal first")
axs[1].set_title("Vertical first")
axs[2].set_title("Difference")

plt.tight_layout()
_images/80a91555c675bce1b384624e4eec936861a53ed227c4aa059fa0b5c73eeaac15.png

Recommendation: horizontal_first#

For temperature, salinity, and any water-mass analysis, use horizontal_first (the default): it averages the way tracers actually mix. Reach for vertical_first only when geometric fidelity to a z-level product near topography outweighs water-mass integrity.

horizontal_first (default)

vertical_first

Lateral averaging

along native isopycnal layers

at constant depth

Water masses / T–S

preserved

smeared across tilted isopycnals

Shelf / slope coverage

good

best (per-column bathymetry)

Global integral

conserved

conserved

Use when

T/S & water-mass budgets (most cases)

matching z-level geometry near topography

Method options#

Argument

Default

Alternatives

method (horizontal)

"conservative"

"bilinear", "patch"

z_method (vertical)

"conservative"

"linear"

  • Conservative (default) conserves area / depth integrals, use for budgets and fluxes. Horizontal needs source cell corners, so pass grid=.

  • Bilinear / linear interpolate point values, smoother, no qlon/qlat needed.

  • Target mask: when target= carries a mask, its land points become NaN; disable with apply_target_mask=False.

The regridding is conservative#

Vertical: the transform thickness-weights each field, so a column’s depth integral (e.g. heat content \(\int T\,dz\)) is preserved.

zlev = np.arange(2.0, 6500, 4.0)
ds_z = xhycom.regrid_vertical(ds, depth=zlev)                 # conservative (default)

# Column integral ∫T dz on native hybrid layers vs on z-levels (one snapshot):
native  = (ds["temp"] * ds["thknss"]).isel(time=0).sum("k").compute()   # hybrid layers
on_zlev = (ds_z["temp"] * 4.0).isel(time=0).sum("depth").compute()     # z-levels (Δz = 4 m)

fig, axs = plt.subplots(3, 1, figsize=(7, 15))

native.plot(ax=axs[0])
on_zlev.plot(ax=axs[1])
(native - on_zlev).plot(ax=axs[2])

axs[0].set_title("On native hybrid layers")
axs[1].set_title("On z-levels")
axs[2].set_title("Difference")

plt.tight_layout()
_images/f588232c23d70abb76aff6092c073ce5c8b24025c29580eef270c8e4dda83c45.png

The absolute difference between the integral computed on native hybrid layers and on regridded z-levels is orders of magnitude smaller than the signal (compare first and second amplitudes versus the third).

Horizontal: conservative remapping reproduces a uniform field exactly (a defining property), so it conserves the area integral. Regrid a field of ones and map the deviation from 1, it sits at machine precision (\(\sim 10^{-12}\)):

ones = ds[["temp"]].isel(time=0).copy()
ones["temp"] = xr.ones_like(ones["temp"])
out = xhycom.regrid_horizontal(ones, target=glorys, grid=grid, apply_target_mask=False)        # conservative (default)

# A uniform field must remap back to 1 everywhere: deviation is at machine precision:
(out["temp"].isel(k=0) - 1).plot(robust=True)                 # ~1e-12

<matplotlib.collections.QuadMesh at 0x14d4f47aed50>
_images/c85c86fcaffffd6ff2a90120d0802746b40fc1f970cd99a670343dd5c6ae3a09.png

Velocities are rotated to east/north#

HYCOM stores velocities on a staggered Arakawa C-grid, with components along the model grid axes: u-vel. runs along the grid’s x-direction and v-vel. along y. On a curvilinear grid those axes are not east/north, they rotate across the domain (sharply near the grid’s poles). regrid / regrid_horizontal de-stagger the components to the T-points and rotate them onto the true geographic axes using the grid angle pang, so the output u-vel. / v-vel. are genuinely eastward / northward.

First, the native components along the model axes:

# Native components, plotted on their own staggered lon/lat: along the model axes.
fig, axs = plt.subplots(1, 2, figsize=(14, 4))
ds["u-vel."].isel(time=0, k=0).plot(ax=axs[0], cmap="RdBu_r", center=0)
ds["v-vel."].isel(time=0, k=0).plot(ax=axs[1], cmap="RdBu_r", center=0)
axs[0].set_title("native u-vel., along model x-axis")
axs[1].set_title("native v-vel., along model y-axis")
fig.tight_layout()
_images/67b6aea7248f927c49f266a0a0e87ddcf610041049c2ebe6bb8c8ca66424341d.png

After regridding, u-vel. / v-vel. carry standard_name = eastward/northward_sea_water_velocity (and a comment recording the de-stagger and rotation). The real test that the rotation is physical, not just a relabel, is that on the shared GLORYS grid they line up with GLORYS’s own eastward (uo) / northward (vo) velocities:

# Surface velocities on the shared GLORYS grid. 
hu = ds_g["u-vel."].isel(time=0, depth=0).compute()
hv = ds_g["v-vel."].isel(time=0, depth=0).compute()
gu = glorys_data["uo"].isel(time=0, depth=0)
gv = glorys_data["vo"].isel(time=0, depth=0)

fig, axs = plt.subplots(2, 2, figsize=(22, 8))
kw = dict(robust=True, cmap="RdBu_r", center=0)

hu.plot(ax=axs[0, 0], **kw)
gu.plot(ax=axs[0, 1], **kw)
hv.plot(ax=axs[1, 0], **kw)
gv.plot(ax=axs[1, 1], **kw)

titles = [["HYCOM eastward (u-vel.)", "GLORYS eastward (uo)"],
          ["HYCOM northward (v-vel.)", "GLORYS northward (vo)"]]

for r in range(2):
    for c in range(2):
        axs[r][c].set_title(titles[r][c])
fig.tight_layout()


for ax in axs.flatten():
    ax.set_ylim(50, None)      # latitude > 50
_images/2bbc6a8de76dd0101519d4aad48a43b43765820d7fb965c2f48df98031241bbb.png

The GLORYS reanalysis shows substantially more mesoscale eddy activity than the HYCOM TOPAZ2 simulation due to its higher spatial resolution. Nevertheless, the large-scale velocity patterns and flow directions remain consistent between the two (compare left and right panels), confirming that xhycom performs the velocity rotation via .regrid and .regrid_horizontal correctly.

Caching the regrid weights#

Since the weights depend only on the grids, you can build them once and reuse them. Pass weights=True to cache the weight matrix to disk; later regrids of other files (or other runs in the same session) read it back instead of regenerating it.

The cache is keyed by grid geometry, the source idm/jdm and lon/lat extent (straight from regional.grid) plus the GLORYS shape/extent, so TP0/TP2/TP5 × target × method each get their own file under $XHYCOM_CACHE_DIR (default ~/.cache/xhycom/regrid_weights). Point XHYCOM_CACHE_DIR at project or scratch space to share weights across jobs.

import os

# Where to cache the xESMF weights. Set this *before* the first weights=True call.
os.environ["XHYCOM_CACHE_DIR"] = f"/cluster/projects/nn2993k/{os.environ['USER']}/.xhycom-cache-dir"

# First call: build the weights and write them to the cache.
%time ds_g = xhycom.regrid(ds, target=glorys, grid=grid, weights=True)

CPU times: user 49 s, sys: 7.13 s, total: 56.1 s
Wall time: 56.5 s

Re-running with the same grids reads the weights back from disk instead of regenerating them. Constructing the regrid is then cheap, and the field values stay lazy: weight generation is the only eager step; the actual remap streams per time chunk when you load.

call

weight matrix

field application

first (weights=True)

built once, written to cache

lazy, per time chunk

later files / sessions

read from cache

lazy, per time chunk

 
# Second call reuses the cached weights; the result is dask-backed (lazy).
%time ds_g = xhycom.regrid(ds, target=glorys, grid=grid, weights=True)
ds_g # repr shows dask chunks, nothing computed yet

CPU times: user 25.4 s, sys: 7.11 s, total: 32.5 s
Wall time: 32.7 s

<xarray.Dataset> Size: 52TB
Dimensions:     (time: 366, lat: 624, lon: 4320, depth: 50, ki: 51)
Coordinates:
  * time        (time) object 3kB 2020-01-01 00:00:00 ... 2020-12-31 00:00:00
  * lat         (lat) float32 2kB 38.08 38.17 38.25 38.33 ... 89.83 89.92 90.0
  * lon         (lon) float32 17kB -180.0 -179.9 -179.8 ... 179.8 179.8 179.9
  * depth       (depth) float64 400B 0.494 1.541 2.646 ... 5.275e+03 5.728e+03
  * ki          (ki) int64 408B 0 1 2 3 4 5 6 7 8 ... 42 43 44 45 46 47 48 49 50
Data variables: (12/82)
    montg1      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    srfhgt      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    oneta       (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    surflx      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    wtrflx      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    salflx      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    ...          ...
    light_sw    (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    attenuat    (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    total_ch    (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    total_ca    (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    u-vel.      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    v-vel.      (time, lat, lon, depth) float64 395GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
Attributes:
    iversn:         23
    iexpt:          28
    yrflag:         3
    archive_type:   mean
    regrid_method:  conservative

The other direction: regrid GLORYS onto the HYCOM grid#

The comparison above brought HYCOM onto the GLORYS grid. regrid_to_hycom does the inverse: it interpolates a regular lon/lat product like GLORYS onto HYCOM’s native curvilinear (y, x) grid. This is the natural direction when the model is coarser than the product, regridding HYCOM up onto the finer GLORYS grid mostly interpolates and adds no information, whereas bringing GLORYS down onto the model grid keeps the comparison in the model’s own space.

Only the horizontal grid changes; fields keep their own depth coordinate. Pass like= a HYCOM field to mask the result to the model’s land/sea mask, and weights=True to cache the remap weights (see Caching the regrid weights).

# GLORYS (regular lon/lat) -> HYCOM curvilinear (y, x). bilinear because the
# model grid is coarser; like=ds masks to HYCOM land; weights=True caches.
glorys_on_hycom = xhycom.regrid_to_hycom(
    glorys_data, grid, method="bilinear", periodic=True, like=ds, weights=True
)
glorys_on_hycom

<xarray.Dataset> Size: 247MB
Dimensions:  (time: 1, depth: 50, y: 380, x: 400)
Coordinates:
  * time     (time) datetime64[ns] 8B 2020-01-01
  * depth    (depth) float32 200B 0.494 1.541 2.646 ... 5.275e+03 5.728e+03
    lon      (y, x) float64 1MB -94.75 -94.58 -94.42 ... 99.09 98.95 98.82
    lat      (y, x) float64 1MB 39.06 39.16 39.26 39.36 ... 56.39 56.29 56.2
Dimensions without coordinates: y, x
Data variables:
    so       (time, depth, y, x) float64 61MB nan nan nan nan ... nan nan nan
    thetao   (time, depth, y, x) float64 61MB nan nan nan nan ... nan nan nan
    vo       (time, depth, y, x) float64 61MB nan nan nan nan ... nan nan nan
    uo       (time, depth, y, x) float64 61MB nan nan nan nan ... nan nan nan
    zos      (time, y, x) float64 1MB nan nan nan nan nan ... nan nan nan nan
Attributes:
    source:                    MERCATOR GLORYS12V1
    title:                     daily mean fields from Global Ocean Physics An...
    Conventions:               CF-1.4
    comment:                   CMEMS product
    institution:               MERCATOR OCEAN
    history:                   2023/06/01 16:20:05 MERCATOR OCEAN Netcdf crea...
    references:                http://www.mercator-ocean.fr
    copernicusmarine_version:  2.0.1
    regrid_method:             bilinear
# Scalars (T/S/SSH) are exact. Both fields now live on HYCOM's (y, x) grid,
# so they difference directly: no separate alignment step.
hycom_sst = ds["temp"].isel(time=0, k=0).compute()
glorys_sst = glorys_on_hycom["thetao"].isel(time=0, depth=0).compute()

fig, axs = plt.subplots(3, 1, figsize=(7, 15))

hycom_sst.plot(ax=axs[0], cmap="inferno", center=False)
glorys_sst.plot(ax=axs[1], cmap="inferno", center=False)
(hycom_sst - glorys_sst).plot(ax=axs[2], cmap="RdBu_r", center=0, robust=True)

axs[0].set_title("HYCOM SST (native grid)")
axs[1].set_title("GLORYS SST on HYCOM grid")
axs[2].set_title("HYCOM − GLORYS")
fig.tight_layout()
_images/5adc7bf4d7be5043bd1397f68c35ed1acb25f22954f246a8256a9e0944eaf868.png

Velocities need rotating to a common frame#

regrid_to_hycom interpolates uo/vo as scalars, so GLORYS velocities land on the HYCOM T-points but stay as geographic eastward/northward components, they are not HYCOM’s native model-axis components on the staggered C-grid. HYCOM’s own u-vel./v-vel. run along the (rotated) model axes, so differencing them against GLORYS directly is wrong wherever the grid turns, i.e. across the whole Arctic cap.

velocities_east_north closes the gap: it de-staggers HYCOM’s velocities to the T-points and rotates them to true east/north using the grid angle pang, keeping the native (y, x) grid. Both fields are then eastward/northward on the same grid and difference directly.

ds["u-vel."].attrs

{'long_name': 'sea water x velocity',
 'units': 'm s-1',
 'hycom_velocity': 'total',
 'comment': 'total current (baroclinic + barotropic); the barotropic part was summed in when the archm time mean was formed'}

# Rotate HYCOM's model-axis velocities to true east/north on the native grid.
hyc = xhycom.velocities_east_north(ds, grid)

hyc["u-vel."].attrs

{'long_name': 'eastward component of sea water x velocity',
 'units': 'm s-1',
 'hycom_velocity': 'total',
 'comment': 'de-staggered to T-points and rotated to geographic axes',
 'standard_name': 'eastward_sea_water_velocity'}

# Both are now eastward velocity on HYCOM's (y, x) grid -> difference directly.
hu = hyc["u-vel."].isel(time=0, k=0).compute()
gu = glorys_on_hycom["uo"].isel(time=0, depth=0).compute()

fig, axs = plt.subplots(3, 1, figsize=(7, 15))
hu.plot(ax=axs[0], cmap="RdBu_r", center=0, robust=True)
gu.plot(ax=axs[1], cmap="RdBu_r", center=0, robust=True)
(hu - gu).plot(ax=axs[2], cmap="RdBu_r", center=0, robust=True)
axs[0].set_title("HYCOM eastward (rotated)")
axs[1].set_title("GLORYS eastward (uo) on HYCOM grid")
axs[2].set_title("HYCOM − GLORYS")
fig.tight_layout()
_images/0d1dc6f196a4d75e43b7320486023935b8144b6f3688588ce19f2888e26b6826.png