UXARRAY / uxarray

Xarray-styled package for reading and directly operating on unstructured grid datasets following UGRID conventions
https://uxarray.readthedocs.io/
Apache License 2.0
142 stars 30 forks source link

subset.bounding_circle error during plotting #848

Open mathomp4 opened 2 weeks ago

mathomp4 commented 2 weeks ago

Version

v2024.06.0

How did you install UXarray?

Pip

What happened?

With the support for GEOS CS files via #802 in v2024.06.0, I decided to try it out and duplicate the plots by @philipc2 as seen here.

I was able to do the global plot, but the regional plot failed with:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[6], line 1
----> 1 uxds['RH'][0][0].subset.bounding_circle((-90, 25), 10).plot.rasterize(method='polygon', exclude_antimeridian=True, pixel_ratio=4.0) * features

File ~/installed/Core/GEOSpyD/24.4.0-0_py3.11/2024-06-11/lib/python3.11/site-packages/uxarray/subset/dataarray_accessor.py:85, in DataArraySubsetAccessor.bounding_circle(self, center_coord, r, element, **kwargs)
     70 """Subsets an unstructured grid by returning all elements within some
     71 radius (in degrees) from a center coord.
     72 
   (...)
     80     Element for use with `coords` comparison, one of `nodes`, `face centers`, or `edge centers`
     81 """
     82 grid = self.uxda.uxgrid.subset.bounding_circle(
     83     center_coord, r, element, **kwargs
     84 )
---> 85 return self.uxda._slice_from_grid(grid)

File ~/installed/Core/GEOSpyD/24.4.0-0_py3.11/2024-06-11/lib/python3.11/site-packages/uxarray/core/dataarray.py:1092, in UxDataArray._slice_from_grid(self, sliced_grid)
   1087 else:
   1088     raise ValueError(
   1089         "Data variable must be either node, edge, or face centered."
   1090     )
-> 1092 return UxDataArray(
   1093     uxgrid=sliced_grid,
   1094     data=d_var,
   1095     name=self.name,
   1096     coords=self.coords,
   1097     dims=self.dims,
   1098     attrs=self.attrs,
   1099 )

File ~/installed/Core/GEOSpyD/24.4.0-0_py3.11/2024-06-11/lib/python3.11/site-packages/uxarray/core/dataarray.py:71, in UxDataArray.__init__(self, uxgrid, *args, **kwargs)
     68 else:
     69     self.uxgrid = uxgrid
---> 71 super().__init__(*args, **kwargs)

File ~/installed/Core/GEOSpyD/24.4.0-0_py3.11/2024-06-11/lib/python3.11/site-packages/xarray/core/dataarray.py:455, in DataArray.__init__(self, data, coords, dims, name, attrs, indexes, fastpath)
    453 data = _check_data_shape(data, coords, dims)
    454 data = as_compatible_data(data)
--> 455 coords, dims = _infer_coords_and_dims(data.shape, coords, dims)
    456 variable = Variable(dims, data, attrs, fastpath=True)
    458 if not isinstance(coords, Coordinates):

File ~/installed/Core/GEOSpyD/24.4.0-0_py3.11/2024-06-11/lib/python3.11/site-packages/xarray/core/dataarray.py:194, in _infer_coords_and_dims(shape, coords, dims)
    191             var.dims = (dim,)
    192             new_coords[dim] = var.to_index_variable()
--> 194 _check_coords_dims(shape, new_coords, dims_tuple)
    196 return new_coords, dims_tuple

File ~/installed/Core/GEOSpyD/24.4.0-0_py3.11/2024-06-11/lib/python3.11/site-packages/xarray/core/dataarray.py:128, in _check_coords_dims(shape, coords, dim)
    126 for d, s in v.sizes.items():
    127     if s != sizes[d]:
--> 128         raise ValueError(
    129             f"conflicting sizes for dimension {d!r}: "
    130             f"length {sizes[d]} on the data but length {s} on "
    131             f"coordinate {k!r}"
    132         )

ValueError: conflicting sizes for dimension 'n_face': length 18392 on the data but length 437400 on coordinate 'lons'

As I'm using the same code and the same file, my guess is maybe I'm missing some package? (I'm still a newbie at uxarray).

What did you expect to happen?

I was expecting to duplicate the image: image

Can you provide a MCVE to repoduce the bug?

import uxarray as ux
import geoviews.feature as gf
import cartopy.crs as ccrs

features = gf.coastline(
    projection=ccrs.PlateCarree(), line_width=1, scale="50m"
) * gf.states(projection=ccrs.PlateCarree(), line_width=1, scale="50m")

file_path = r"test-fixapp-c270.geosgcm_prog_nat.20200415_0000z.nc4"
uxds = ux.open_dataset(file_path, file_path)
uxds['RH'][0][0].subset.bounding_circle((-90, 25), 10).plot.rasterize(method='polygon', exclude_antimeridian=True, pixel_ratio=4.0) * features
philipc2 commented 2 weeks ago

Do you mind sharing what uxds['RH'] looks like?

EDIT:

Forgot I had this dataset, I'll look into it this afternoon and see if I can replicate the issue