[Bug]: lon_orient=(0, 360) tacking on an extra longitude

[mzelinka]

What happened?

I am reading in a dataset with 360 longitudes that span -179.5 to 179.5 by 1 degree. Instead I want the data to run from 0.5 to 359.5 by 1 degree. Using ds = xc.open_dataset(filename, use_cftime = True, lon_orient=(0, 360)), it successfully reorients the longitudes to be positive definite (in degrees East), but it adds an additional longitude such that I have both 359.5 and 360.

What did you expect to happen? Are there are possible answers you came across?

I expected this to keep the same longitudes, but orient the data from 0.5 to 359.5 without adding that 360 degree point.

Minimal Complete Verifiable Example (MVCE)

On NERSC:

filename='/global/cfs/cdirs/m4581/obs4MIPs/obs4MIPs_input/MOHC/HadISST1-1/v20231114/HadISST.1.1.sst.nc'
ds0 = xc.open_mfdataset(filename, use_cftime = True)
ds1 = xc.open_mfdataset(filename, use_cftime = True, lon_orient=(0, 360))
ds0.longitude,ds1.longitude

Relevant log output

(<xarray.DataArray 'longitude' (longitude: 360)> Size: 1kB
 array([-179.5, -178.5, -177.5, ...,  177.5,  178.5,  179.5], dtype=float32)
 Coordinates:
   * longitude  (longitude) float32 1kB -179.5 -178.5 -177.5 ... 178.5 179.5
 Attributes:
     units:          degrees_east
     long_name:      Longitude
     standard_name:  longitude
     bounds:         longitude_bnds,
 <xarray.DataArray 'longitude' (longitude: 361)> Size: 1kB

 array([  0.5,   1.5,   2.5, ..., 358.5, 359.5, 360. ], dtype=float32)
 Coordinates:
   * longitude  (longitude) float32 1kB 0.5 1.5 2.5 3.5 ... 358.5 359.5 360.0
 Attributes:
     units:          degrees_east
     long_name:      Longitude
     standard_name:  longitude
     bounds:         longitude_bnds)

Anything else we need to know?

No response

Environment

xCDAT version 0.8.0

INSTALLED VERSIONS

commit: None
python: 3.13.2 | packaged by conda-forge | (main, Feb 17 2025, 14:10:22) [GCC 13.3.0]
python-bits: 64
OS: Linux
OS-release: 5.14.21-150500.55.97_13.0.78-cray_shasta_c
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: ('en_US', 'UTF-8')
libhdf5: 1.14.4
libnetcdf: 4.9.2

xarray: 2025.3.0
pandas: 2.2.3
numpy: 2.1.3
scipy: 1.15.2
netCDF4: 1.7.2
pydap: None
h5netcdf: None
h5py: None
zarr: None
cftime: 1.6.4
nc_time_axis: 1.4.1
iris: None
bottleneck: 1.4.2
dask: 2025.3.0
distributed: 2025.3.0
matplotlib: 3.10.1
cartopy: 0.24.1
seaborn: None
numbagg: None
fsspec: 2025.3.2
cupy: None
pint: 0.24.4
sparse: 0.15.4
flox: None
numpy_groupies: None
setuptools: 78.1.1
pip: 25.1
conda: None
pytest: None
mypy: None
IPython: 9.0.2
sphinx: None

[pochedls]

I haven't fully tracked this down, but the extra longitude value emerges here when xcdat converts longitude_bnds to (0, 360). A temporary work around is:

# imports
import xcdat as xc
# parameters
filename='/global/cfs/cdirs/m4581/obs4MIPs/obs4MIPs_input/MOHC/HadISST1-1/v20231114/HadISST.1.1.sst.nc'
# open dataset
# I think use_cftime = True is the default and doesn't need to be specified
# also note you don't need to use mfdataset since this is a single file
ds = xc.open_dataset(filename)  
# drop the bounds to avoid the problem
ds = ds.drop_vars('longitude_bnds')
# do conversion
ds = xc.swap_lon_axis(ds, to=(0, 360))
# add the bounds back
ds = ds.bounds.add_missing_bounds()

[mzelinka]

Thanks for the temporary workaround!

[tomvothecoder]

Thanks for reporting this issue @mzelinka, and thank you @pochedls for a temporary workaround.

I think it might be related to the code that handles the prime meridian, which extends longitude and bounds by one cell:

xcdat/xcdat/axis.py

Lines 375 to 384 in d9d0e76

	# Handle cases where a prime meridian cell exists, which can occur
	# after swapping longitude bounds to (0, 360). This involves extending
	# the longitude and bounds by one cell to take into account the prime
	# meridian. It also results in extending the data variables by one
	# value.
	if to == (0, 360):
	p_meridian_index = _get_prime_meridian_index(ds[key])
	if p_meridian_index is not None:
	ds = _align_lon_to_360(ds, ds[key], p_meridian_index)

[pochedls]

The issue has roots in the spatial averager. Some of that functionality was applied for this case (swap_lon_axis). The spatial averager gets the longitudinal component of the weights by computing: lon_bounds[:, 1] - lon_bounds[:, 0]. But this is problematic if a bound spans the prime meridian (e.g., [359, 1] because the weight for that grid cell would be -358 instead of 2). The spatial averager handles this by breaking this into two pieces: [0, 1] and [359, 360]. This results in the weights dataarray of having an extra longitude index, but the extra index is added to the original prime meridian cell and then the extra index is dropped in the weights dataarray.

In the _swap_lon_axis (within swap_lon_axis) implementation this final step isn't done (dropping the extra index) for the longitude bounds array.

Adding the following here mostly fixes the problem:

                dim = get_dim_keys(ds[key], "X")
                ds = ds.isel({dim: slice(0, -1)})

But the bounds are not correct for the prime meridian index.

But I'm not entirely sure any of this is needed. If you drop this block, I think it yields the correct bounds. So we might be able to drop this block of code. One residual issue is that if you do this, the last bounds is [359, 0] and it should be [359, 360]. Maybe this could be fixed with a separate function?

[tomvothecoder]

@pochedls I followed your suggestions above to:

1. remove the prime meridian code here
2. add a separate function to fix [359, 0] with [359, 360].

These changes fix the issue in this GitHub description. However, it creates another axis conversion issue for this unit test case (which has a cell wrapping the prime meridian):

                 [-180.5, -179.5],
                 [-1.5, -0.5],
                 [-0.5, 0.5],
                 [0.5, 1.5],
                 [180.5, 179.5],

So, for the bounds [-0.5, 0.5]:

• -0.5 becomes 359.5 (-0.5 + 360)
• 0.5 stays 0.5

Therefore, the interval [-0.5, 0.5] in [-180, 180] becomes [359.5, 0.5] in [0, 360]. This interval now wraps around the 0-degree meridian.

Maybe we still need the prime meridian code and to follow your suggestion here to drop the last value with additional code to fix the bounds:

                dim = get_dim_keys(ds[key], "X")
                ds = ds.isel({dim: slice(0, -1)})

[tomvothecoder]

After consulting with ChatGPT:

TL;DR: For longitude orientation swaps, don’t handle prime-meridian cells. It should be a pure modulo + sort.

• Keep the cell count unchanged; centers in a half-open range ([0,360) or (−180,180]).
• Map centers with % 360 (or to target interval) and sortby.
• Map bounds element-wise with the same modulo; it’s fine if the PM cell has lo > hi in [0,360).
• Don’t split cells or append a 360.0 center—splitting is only needed in area-aware ops (spatial averaging/regridding).
• The extra 360.0 came from PM-split logic leaking into orientation; fix by dropping that last appended cell (short term) and removing PM-split from the orientation path (long term), gating it in the averager instead.

@pochedls Does this sound right?

[tomvothecoder]

One residual issue is that if you do this, the last bounds is [359, 0] and it should be [359, 360].

I'm reading that it is fine for the last longitude cell to have a bounds like [359, 0] (for cell 359.5) after swapping from (-180, 180) -> (0, 360).

There are two equivalent ways to describe that last cell, both CF-compatible:

• [359, 0] → emphasizes the wrap across the prime meridian (lower > upper).
• [359, 360] → emphasizes a closed upper edge at 360, which is the same meridian as 0.

Both are valid because longitude is modulo 360 — 0° and 360° refer to the same line. CF only requires that:

• the longitude coordinates themselves are strictly monotonic, and
• the bounds follow the same sense as the coordinates.

So [359, 0] for the cell centered at 359.5° is fine and interpretable by CF-compliant software. Some tools normalize it to [359, 360], others leave it as [359, 0], but they mean the same thing.

👉 If you want to be maximally “conventional,” you can normalize to [359, 360]. But if you leave it [359, 0], it’s still correct and doesn’t break CF compliance.

@pochedls Is this right too? And does the spatial average correctly handle these bounds ([359, 0]) using the prime meridian code?

[pochedls]

I haven't fully digested this thread, but to give quick feedback: I don't think ChatGPT is right on all of this (see this thread).

I suspect that we should follow the final part of my comment here: drop the prime meridian index code and then write a separate function to fix bounds that span the prime meridian (e.g., [359, 360] stays as-is; [359, 1] goes to [-1, 1] or [359, 361]).

[tomvothecoder]

I haven't fully digested this thread, but to give quick feedback: I don't think ChatGPT is right on all of this (see this thread).

I suspect that we should follow the final part of my comment here: drop the prime meridian index code and then write a separate function to fix bounds that span the prime meridian (e.g., [359, 360] stays as-is; [359, 1] goes to [-1, 1] or [359, 361]).

@pochedls Thanks. That makes sense to me. I'm going to push code to #791 with these changes.