[Bug]: "More than one grid cell spans prime meridian" message when doing spatial.average

[jlee1046]

What happened?

I am trying to compute an area weighted global mean of precipitation using original monthly mean IMERG GPM precipitation data. The data is rectilinear (1800x3600). While spatial.average function is executed, I get the error message as shown below:
ValueError: More than one grid cell spans prime meridian.

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

I expected to get 12 (monthly) values of area weighted precipitation mean stored in 'case_data'

Minimal Complete Verifiable Example (MVCE)

case_dir = "/global/cfs/cdirs/m3312/jlee1046/GPM/GPM_Cess_2019Sep-2020Aug/Monthly"
file2search = f'{case_dir}/3B-MO.MS.MRG.3IMERG.*.V07B.HDF5.nc4'
case_flist = glob.glob(file2search)
case_flist.sort()

case_ds = xcdat.open_mfdataset(case_flist)
case_ds.lat.attrs["axis"] = "Y"
case_ds.lon.attrs["axis"] = "X"
case_data = case_ds.spatial.average('precipitation')["precipitation"]

Relevant log output

Traceback (most recent call last):
  File "/global/cfs/cdirs/m3312/jlee1046/SCREAM_Cess/plotting_scripts/compute_area_weighted_mean.py", line 26, in <module>
    case_data = case_ds.spatial.average('precipitation')["precipitation"]
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/global/homes/j/jlee1046/.conda/envs/xcdat_env/lib/python3.11/site-packages/xcdat/spatial.py", line 191, in average
    self._weights = self.get_weights(axis, lat_bounds, lon_bounds, data_var)
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/global/homes/j/jlee1046/.conda/envs/xcdat_env/lib/python3.11/site-packages/xcdat/spatial.py", line 286, in get_weights
    weights = axis_bounds[key]["weights_method"](d_bounds, r_bounds)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/global/homes/j/jlee1046/.conda/envs/xcdat_env/lib/python3.11/site-packages/xcdat/spatial.py", line 429, in _get_longitude_weights
    p_meridian_index = _get_prime_meridian_index(d_bounds)
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/global/homes/j/jlee1046/.conda/envs/xcdat_env/lib/python3.11/site-packages/xcdat/axis.py", line 429, in _get_prime_meridian_index
    raise ValueError("More than one grid cell spans prime meridian.")
ValueError: More than one grid cell spans prime meridian.

Anything else we need to know?

I changed the file permission. You should be able to access the IMERG monthly mean files.

Environment

xCDAT version = 0.6.1

INSTALLED VERSIONS

commit: None
python: 3.11.7 | packaged by conda-forge | (main, Dec 23 2023, 14:43:09) [GCC 12.3.0]
python-bits: 64
OS: Linux
OS-release: 5.14.21-150400.24.81_12.0.87-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.3
libnetcdf: 4.9.2

xarray: 2023.12.0
pandas: 2.1.4
numpy: 1.26.3
scipy: 1.11.4
netCDF4: 1.6.5
pydap: None
h5netcdf: None
h5py: None
Nio: None
zarr: None
cftime: 1.6.3
nc_time_axis: None
iris: None
bottleneck: None
dask: 2024.1.0
distributed: 2024.1.0
matplotlib: None
cartopy: None
seaborn: None
numbagg: None
fsspec: 2023.12.2
cupy: None
pint: None
sparse: 0.15.1
flox: None
numpy_groupies: None
setuptools: 69.0.3
pip: 23.3.2
conda: None
pytest: None
mypy: None
IPython: 8.20.0
sphinx: None

[pochedls]

@jlee1046 – could you inspect the longitude bounds to see if more than one of the bounds spans 0 degrees longitude? Or to see if they are strange in some other way?

If so, you might be able to drop the longitude bounds and re-generate them with add_missing_bounds.

[jlee1046]

@pochedls longitude bound of first point is ( -180, -179.9) and the bound for the last point is (179.9, 180). in the middle, the bounds are (-0.09999999, 3.72529e-09) and (-3.72529e-09, 0.09999999). Does this mean that I have more than on of the bounds spans 0 degree?

[pochedls]

Yes – the grid cells that have the left bound west of 0^o longitude and the right bound east of 0^o longitude both encompass the prime meridian (and have overlapping bounds). You could manually fix the bounds or drop the bounds and have xcdat regenerate them (and then take the spatial average).

[lee1043]

@jlee1046 Following @pochedls's suggestion, I have drafted a potential code you might be able to use. Can you try if this would work for you? variable name "lat_bnds" and "lon_bnds" may or may not need to be fixed depending how your variable name defined in the file.

case_ds = xcdat.open_mfdataset(case_flist)
case_ds.lat.attrs["axis"] = "Y"
case_ds.lon.attrs["axis"] = "X"
# Drop bounds
case_ds_bnds_dropped = case_ds.drop(["lat_bnds", "lon_bnds"])
# Regenerate bounds
case_ds = case_ds_bnds_dropped.bounds.add_missing_bounds()
# Spatial average
case_data = case_ds.spatial.average('precipitation')["precipitation"]