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

jlee1046 commented 4 months ago

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 commented 4 months ago

@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 commented 4 months ago

@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 commented 4 months ago

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 commented 4 months ago

@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"]

tomvothecoder commented 2 weeks ago

Error is correct and workaround provided above. Closing this issue.

xCDAT / xcdat