Regridding sensitivity to alignment of axis/bounds orientation

Hi there! I'm a core developer of the xCDAT package, which provides a horizontal regridding API that extends the capabilities of the xESMF regridding interface.

Overview

While investigating an issue in the E3SM Diags repository, we discovered a key detail about axis orientation during regridding. Specifically, in this comment, we observed that if the orientation of the coordinate axis and its bounds don’t match (e.g., one is ascending while the other is descending), it can lead to incorrect regridding results.

In contrast, as noted in this follow-up comment, the CDAT regrid2 package handles this more robustly. It derives 1D arrays for each direction (south, north, west, east) directly from the bounds. This approach makes the axis orientation irrelevant, producing consistent results regardless of the order.

Proposed Solution

Is it reasonable to add a function that can handle extract bounds similar to regrid2? We think the regridder should be robust to handle different orientations across grids. Thanks!

Log Output

Notice how the first two rows with misaligned orientations have results that are off, compared to aligned orientations represented by the last four rows. For the last two rows that have no lat_bnds, we drop existing bounds and generate new bounds using xCDAT .

Statistical Comparison:
                              Min        Max      Mean           Sum       Std
asc lat, desc lat_bnds   0.165452   9.341431  1.423206  20494.160156  1.052469
desc lat, asc lat_bnds   0.165457   9.341392  1.374420  19791.644531  1.021888
asc lat, asc lat_bnds    0.164364  10.046874  1.398573  20139.449219  1.039437
desc lat, desc lat_bnds  0.164364  10.046874  1.398573  20139.449219  1.039437
desc lat, no lat_bnds    0.164364  10.046874  1.398573  20139.449219  1.039437
asc lat, no lat_bnds     0.164364  10.046874  1.398573  20139.449219  1.039437

Minimum Script

The script is reproducible as the data is publicly accessible.

import cdms2
import numpy as np
import pandas as pd
import xarray as xr
import xcdat as xc
import xesmf as xe

from regrid2 import Regridder


def print_stats(*arrays, labels=None):
    """Prints statistical comparison of multiple arrays using a pandas DataFrame."""
    if labels is None:
        labels = [f"Array {i + 1}" for i in range(len(arrays))]
    elif len(labels) != len(arrays):
        raise ValueError("Number of labels must match the number of arrays.")

    stats = {
        "Min": [np.min(arr) for arr in arrays],
        "Max": [np.max(arr) for arr in arrays],
        "Mean": [np.mean(arr) for arr in arrays],
        "Sum": [np.sum(arr) for arr in arrays],
        "Std": [np.std(arr) for arr in arrays],
    }

    # Create a DataFrame from the stats dictionary
    df = pd.DataFrame(stats, index=labels)

    # Print the DataFrame
    print("\nStatistical Comparison:")
    print(df)

def regrid_with_xesmf(ds, output_grid, method="conservative_normed"):
    """Regrid a dataset using xESMF."""
    print("Regridding with xESMF...")
    print("--" * 20)
    # Print lat and lat_bnds ascending or descending with first 5 values
    lat_ascending = np.all(np.diff(ds["lat"].values) > 0)
    lat_bnds_ascending = np.all(np.diff(ds["lat_bnds"].values, axis=-1) > 0)

    print("Latitude is ascending:", lat_ascending)
    print("First 5 lat values:", ds["lat"].values[:5])
    print("Latitude bounds are ascending:", lat_bnds_ascending)
    print("First 5 lat_bnds values:", ds["lat_bnds"].values[:5])

    regridder = xe.Regridder(ds, output_grid, method=method)
    return regridder(ds["PRECT"])


# %%
# Load datasets
ds_a = xr.open_dataset(
    "/lcrc/group/e3sm/public_html/cdat-migration-fy24/25-02-14-branch-930-polar-after/polar/GPCP_v3.2/test.nc"
)
ds_b = xr.open_dataset(
    "/lcrc/group/e3sm/public_html/cdat-migration-fy24/25-02-14-branch-930-polar-after/polar/GPCP_v3.2/ref.nc"
)

# %%
# 1. xCDAT + xESMF (ascending lat, descending lat_bnds -- default)
# ----------------------------------------------------
output_grid_xesmf = ds_a.regridder.grid
ds_b1 = ds_b.copy(deep=True)
unaligned1 = regrid_with_xesmf(ds_b1, output_grid_xesmf)

# 2. xCDAT + xESMF (descending lat, ascending lat_bnds)
# ----------------------------------------------------
ds_b2 = ds_b.copy(deep=True).sortby("lat", ascending=False)
ds_b2["lat_bnds"].values = np.sort(ds_b2["lat_bnds"], axis=-1)
unaligned2 = regrid_with_xesmf(ds_b2, output_grid_xesmf)

# 3. xCDAT + xESMF (ascending lat, ascending lat_bnds)
# ----------------------------------------------------
ds_b3 = ds_b.copy(deep=True)
ds_b3["lat_bnds"].values = np.sort(ds_b3["lat_bnds"], axis=-1)
aligned1 = regrid_with_xesmf(ds_b3, output_grid_xesmf)

# 4. xCDAT + xESMF (descending lat, descending lat_bnds)
# ----------------------------------------------------
ds_b4 = ds_b.copy(deep=True)
ds_b4 = ds_b4.sortby("lat", ascending=False)
ds_b4["lat_bnds"].values = np.sort(ds_b4["lat_bnds"], axis=-1)[:, ::-1]
aligned2 = regrid_with_xesmf(ds_b4, output_grid_xesmf)

# 5. Desc lat, no lat_bnds
# ----------------------------------------------------
ds_b5 = ds_b.copy(deep=True).sortby("lat", ascending=False)
ds_b5 = ds_b5.drop_vars("lat_bnds")
ds_b5 = ds_b5.bounds.add_bounds(axis="Y")
nobounds1 = regrid_with_xesmf(ds_b5, output_grid_xesmf)

# 6. Asc lat, no lat_bnds
# ----------------------------------------------------
ds_b6 = ds_b.copy(deep=True).sortby("lat", ascending=True)
ds_b6 = ds_b6.drop_vars("lat_bnds")
ds_b6 = ds_b6.bounds.add_bounds(axis="Y")
nobounds2 = regrid_with_xesmf(ds_b6, output_grid_xesmf)

# %%
# Compare statistics
# ----------------------------------------------------
print_stats(
    unaligned1.values,
    unaligned2.values,
    aligned1.values,
    aligned2.values,
    nobounds1.values,
    nobounds2.values,
    labels=[
        "asc lat, desc lat_bnds",
        "desc lat, asc lat_bnds",
        "asc lat, asc lat_bnds",
        "desc lat, desc lat_bnds",
        "desc lat, no lat_bnds",
        "asc lat, no lat_bnds",

    ],
)

@aulemahal

Thanks for the detailed report.

One suspect for this behavior could be https://github.com/xarray-contrib/cf-xarray/blob/4be96ff2e40e7a7fc7ff6d57744c4a42b6cbaa94/cf_xarray/helpers.py#L190
which xESMF relies on to go from 1D bounds like yours to vertices expected by ESMF. The helper function assumes the bounds follow CF-Conventions, which mention that the ordering for bounds should be consistent with the coordinates. See https://cfconventions.org/Data/cf-conventions/cf-conventions-1.12/cf-conventions.html#bounds-one-d

@aulemahal Is this something we should raise with cf-xarray, or add a check for within xESMF ?

import xarray as xr
import cf_xarray as cfxr

bounds = xr.DataArray([[50.5, 50. ], [51.,  50.5], [51.5, 51. ], [52.,  51.5], [52.5, 52.]], dims=('lat', 'bounds'))
cfxr.bounds_to_vertices(bounds, "bounds")

yields

<xarray.DataArray (lat_vertices: 6)> Size: 48B
array([50.5, 51. , 51.5, 52. , 52.5, 52. ])
Dimensions without coordinates: lat_vertices

instead of starting at 50 and ending at 52.5.

This is certainly something that should be raised to cf-xarray. I contributed to that part of cf-xarray, I'll see how we can fix this.

@aulemahal

Hi @aulemahal, let me know if I should open this issue on the cf-xarray repo and if you need help with a fix. I might have time to contribute via a PR.

Yes please open an issue on cf-xarray. I do not have as much time as I expected to fix this, sorry. I expect the fix to be a simple transpose, but I haven't tested anything...

@huard

@huard @aulemahal I opened xarray-contrib/cf-xarray#579 to fix cf_xarray.bounds_to_vertices() so it can handle descending bounds correctly.

Example

import xarray as xr
import cf_xarray as cfxr

# Descending bounds
bounds = xr.DataArray([[50.5, 50. ], [51.,  50.5], [51.5, 51. ], [52.,  51.5], [52.5, 52.]], dims=('lat', 'bounds'))
cfxr.bounds_to_vertices(bounds, "bounds") # array([52.5, 52. , 51.5, 51. , 50.5, 50. ])

# Ascending bounds
bounds = xr.DataArray(
    [[50.0, 50.5], [50.5, 51.0], [51.0, 51.5], [51.5, 52.0], [52.0, 52.5]],
    dims=("lat", "bounds"),
)
cfxr.bounds_to_vertices(bounds, "bounds") # array([50. , 50.5, 51. , 51.5, 52. , 52.5])

According to the tests in cf-xarray, bounds_to_vertices extracts the vertices with an order that is preserved from the original bounds. I think this is the correct behavior.

@huard Does xESMF expects the output vertices to always be ascending instead?

@tomvothecoder

Thanks @tomvothecoder for the PR! What xESMF expects is that the bounds will we in the same order as the coordinate. In your "descending" example, this means we expect the coord to be : 52.25, 51.75, 51.25, 50.75, 50.25. And the opposite for the "ascending" example. Maybe, we need a check for that ? This would be in xESMF indeed.

Thanks @tomvothecoder for the PR! What xESMF expects is that the bounds will we in the same order as the coordinate. In your "descending" example, this means we expect the coord to be : 52.25, 51.75, 51.25, 50.75, 50.25. And the opposite for the "ascending" example. Maybe, we need a check for that ? This would be in xESMF indeed.

Of course! My PR should handle this correctly, as the tests still pass.

I re-ran the minimum example with xarray-contrib/cf-xarray#579 and confirmed it fixes this issue.

Once a new version of cf-xarray is released, we can close this issue.

Statistical Comparison:
                              Min        Max      Mean           Sum       Std
asc lat, desc lat_bnds   0.164364  10.046874  1.398573  20139.449219  1.039437
desc lat, asc lat_bnds   0.164364  10.046874  1.398573  20139.449219  1.039437
asc lat, asc lat_bnds    0.164364  10.046874  1.398573  20139.449219  1.039437
desc lat, desc lat_bnds  0.164364  10.046874  1.398573  20139.449219  1.039437
desc lat, no lat_bnds    0.164364  10.046874  1.398573  20139.449219  1.039437
asc lat, no lat_bnds     0.164364  10.046874  1.398573  20139.449219  1.039437

Min example

import numpy as np
import pandas as pd
import xarray as xr
import xcdat as xc
import xesmf as xe


def print_stats(*arrays, labels=None):
    """Prints statistical comparison of multiple arrays using a pandas DataFrame."""
    if labels is None:
        labels = [f"Array {i + 1}" for i in range(len(arrays))]
    elif len(labels) != len(arrays):
        raise ValueError("Number of labels must match the number of arrays.")

    stats = {
        "Min": [np.min(arr) for arr in arrays],
        "Max": [np.max(arr) for arr in arrays],
        "Mean": [np.mean(arr) for arr in arrays],
        "Sum": [np.sum(arr) for arr in arrays],
        "Std": [np.std(arr) for arr in arrays],
    }

    # Create a DataFrame from the stats dictionary
    df = pd.DataFrame(stats, index=labels)

    # Print the DataFrame
    print("\nStatistical Comparison:")
    print(df)


def regrid_with_xesmf(ds, output_grid, method="conservative_normed"):
    """Regrid a dataset using xESMF."""
    print("\n")
    print("--" * 20)
    # Print lat and lat_bnds ascending or descending with first 5 values
    lat_ascending = np.all(np.diff(ds["lat"].values) > 0)
    lat_bnds_ascending = np.all(np.diff(ds["lat_bnds"].values, axis=-1) > 0)

    print("Latitude is ascending:", lat_ascending)
    print("First 5 lat values:", ds["lat"].values[:5])
    print("Latitude bounds are ascending:", lat_bnds_ascending)
    print("First 5 lat_bnds values:", ds["lat_bnds"].values[:5])

    regridder = xe.Regridder(ds, output_grid, method=method)
    return regridder(ds["PRECT"])


# %%
# Load datasets
ds_a = xr.open_dataset(
    "/lcrc/group/e3sm/public_html/cdat-migration-fy24/25-02-14-branch-930-polar-after/polar/GPCP_v3.2/test.nc"
)
ds_b = xr.open_dataset(
    "/lcrc/group/e3sm/public_html/cdat-migration-fy24/25-02-14-branch-930-polar-after/polar/GPCP_v3.2/ref.nc"
)

# %%
# 1. xCDAT + xESMF (ascending lat, descending lat_bnds -- default)
# ----------------------------------------------------
output_grid_xesmf = ds_a.regridder.grid
ds_b1 = ds_b.copy(deep=True)
unaligned1 = regrid_with_xesmf(ds_b1, output_grid_xesmf)

# 2. xCDAT + xESMF (descending lat, ascending lat_bnds)
# ----------------------------------------------------
ds_b2 = ds_b.copy(deep=True).sortby("lat", ascending=False)
ds_b2["lat_bnds"].values = np.sort(ds_b2["lat_bnds"], axis=-1)
unaligned2 = regrid_with_xesmf(ds_b2, output_grid_xesmf)

# 3. xCDAT + xESMF (ascending lat, ascending lat_bnds)
# ----------------------------------------------------
ds_b3 = ds_b.copy(deep=True)
ds_b3["lat_bnds"].values = np.sort(ds_b3["lat_bnds"], axis=-1)
aligned1 = regrid_with_xesmf(ds_b3, output_grid_xesmf)

# 4. xCDAT + xESMF (descending lat, descending lat_bnds)
# ----------------------------------------------------
ds_b4 = ds_b.copy(deep=True)
ds_b4 = ds_b4.sortby("lat", ascending=False)
ds_b4["lat_bnds"].values = np.sort(ds_b4["lat_bnds"], axis=-1)[:, ::-1]
aligned2 = regrid_with_xesmf(ds_b4, output_grid_xesmf)

# 5. Desc lat, no lat_bnds
# ----------------------------------------------------
ds_b5 = ds_b.copy(deep=True).sortby("lat", ascending=False)
ds_b5 = ds_b5.drop_vars("lat_bnds")
ds_b5 = ds_b5.bounds.add_bounds(axis="Y")
nobounds1 = regrid_with_xesmf(ds_b5, output_grid_xesmf)

# 6. Asc lat, no lat_bnds
# ----------------------------------------------------
ds_b6 = ds_b.copy(deep=True).sortby("lat", ascending=True)
ds_b6 = ds_b6.drop_vars("lat_bnds")
ds_b6 = ds_b6.bounds.add_bounds(axis="Y")
nobounds2 = regrid_with_xesmf(ds_b6, output_grid_xesmf)

# %%
# Compare statistics
# ----------------------------------------------------
print_stats(
    unaligned1.values,
    unaligned2.values,
    aligned1.values,
    aligned2.values,
    nobounds1.values,
    nobounds2.values,
    labels=[
        "asc lat, desc lat_bnds",
        "desc lat, asc lat_bnds",
        "asc lat, asc lat_bnds",
        "desc lat, desc lat_bnds",
        "desc lat, no lat_bnds",
        "asc lat, no lat_bnds",
    ],
)

I've confirmed this issue is now resolved with cf-xarray v0.10.7.

Statistical Comparison:
                              Min        Max      Mean           Sum       Std
asc lat, desc lat_bnds   0.164364  10.046874  1.398573  20139.449219  1.039437
desc lat, asc lat_bnds   0.164364  10.046874  1.398573  20139.449219  1.039437
asc lat, asc lat_bnds    0.164364  10.046874  1.398573  20139.449219  1.039437
desc lat, desc lat_bnds  0.164364  10.046874  1.398573  20139.449219  1.039437
desc lat, no lat_bnds    0.164364  10.046874  1.398573  20139.449219  1.039437
asc lat, no lat_bnds     0.164364  10.046874  1.398573  20139.449219  1.039437

Thank you everyone for the help! Closing this issue now.

Navigation Menu

Repository navigation

Regridding sensitivity to alignment of axis/bounds orientation #433

Overview

Proposed Solution

Log Output

Minimum Script

Example

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Regridding sensitivity to alignment of axis/bounds orientation #433

Description

Overview

Proposed Solution

Log Output

Minimum Script

Activity

huard commented on May 22, 2025

huard commented on May 22, 2025

aulemahal commented on May 22, 2025

tomvothecoder commented on Jun 26, 2025

aulemahal commented on Jun 27, 2025

tomvothecoder commented on Jul 9, 2025

Example

aulemahal commented on Jul 9, 2025

tomvothecoder commented on Jul 9, 2025

tomvothecoder commented on Aug 12, 2025

tomvothecoder commented on Aug 27, 2025

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions