When generalizing a structured dataset with UxDataset.from_structured, return information about faces

[zarzycki]

Proposed new feature or change:

I have a script where I am trying to either load unstructured CAM-SE data (unstructured) or AORC data (structured). The data all represents hourly precipitation.

I want to be able to use uxarray subsetting to calculate regional statistics. The most logical way I can think of is to have an interface layer at the top of the script that reads either the CAM-SE or the AORC data, standardizes them to a ux DataArray, and then goes from there. That way everything below the interface layer just purely works with "unstructured" data (the structured mesh is just a special case of unstructured data)

Currently, when I use the AORC data:

data_filenames = "./aorc/PrecipitationAORC_2004*.nc"
xrds = xr.open_mfdataset(data_filenames)
uxgrid = ux.Grid.from_structured(lon=xrds.longitude, lat=xrds.latitude)
print(uxgrid.face_areas)
uxds = ux.UxDataset.from_xarray(xrds, uxgrid=uxgrid)
print(uxds)

The uxgrid from_structured data does indeed contain a column of face_areas. However, when creating a new uxds using UxDataset.from_xarray the notion of the data being face centered doesn't populate.

<uxarray.Grid>
Original Grid Type: structured
Grid Dimensions:
  * n_node: 35305204
  * n_face: 35292601
  * n_max_face_nodes: 4
Grid Coordinates (Spherical):
  * node_lon: (35305204,)
  * node_lat: (35305204,)
Grid Coordinates (Cartesian):
Grid Connectivity Variables:
  * face_node_connectivity: (35292601, 4)
Grid Descriptor Variables:

<xarray.UxDataset> Size: 68GB
Dimensions:       (time: 240, latitude: 4201, longitude: 8401)
Coordinates:
  * time          (time) datetime64[ns] 2kB 2004-09-04 ... 2004-09-13T23:00:00
  * latitude      (latitude) float64 34kB 20.0 20.01 20.02 ... 54.98 54.99 55.0
  * longitude     (longitude) float64 67kB -130.0 -130.0 -130.0 ... -60.01 -60.0
Data variables:
    APCP_surface  (time, latitude, longitude) float64 68GB dask.array<chunksize=(1, 400, 400), meta=np.ndarray>

This leads to errors when a function needs to query the form of the data (i.e., face-centered, edge-centered, etc.) such as:

radial_subset = hourly_max.subset.bounding_circle(center_coord, r)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[55], line 4
      1 # Subsetting
      2 center_coord = [landfall_lon, landfall_lat]
----> 4 radial_subset = hourly_max.subset.bounding_circle(center_coord, r)
      6 print(radial_subset)

File [~/miniconda3/envs/uxarray/lib/python3.11/site-packages/uxarray/subset/dataarray_accessor.py:93](http://localhost:8888/lab/workspaces/auto-h/tree/~/miniconda3/envs/uxarray/lib/python3.11/site-packages/uxarray/subset/dataarray_accessor.py#line=92), in DataArraySubsetAccessor.bounding_circle(self, center_coord, r, element, inverse_indices, **kwargs)
     73 """Subsets an unstructured grid by returning all elements within some
     74 radius (in degrees) from a center coord.
     75 
   (...)     88     - False: No index storage (default)
     89 """
     90 grid = self.uxda.uxgrid.subset.bounding_circle(
     91     center_coord, r, element, inverse_indices=inverse_indices, **kwargs
     92 )
---> 93 return self.uxda._slice_from_grid(grid)

File [~/miniconda3/envs/uxarray/lib/python3.11/site-packages/uxarray/core/dataarray.py:1706](http://localhost:8888/lab/workspaces/auto-h/tree/~/miniconda3/envs/uxarray/lib/python3.11/site-packages/uxarray/core/dataarray.py#line=1705), in UxDataArray._slice_from_grid(self, sliced_grid)
   1701     da_sliced = self.isel(
   1702         n_node=sliced_grid._ds["subgrid_node_indices"], ignore_grid=True
   1703     )
   1705 else:
-> 1706     raise ValueError(
   1707         "Data variable must be either node, edge, or face centered."
   1708     )
   1710 return UxDataArray(da_sliced, uxgrid=sliced_grid)

ValueError: Data variable must be either node, edge, or face centered.

Maybe I'm missing something obvious!

[philipc2]

@zarzycki

Thanks for sharing!

I think a quick fix would be to flatten APCP_surface to a (time, latitude x longitude) array with a final dimension of n_face. UXarray requires all data variables to be mapped to one of the Grid dimensions, otherwise we can't recognize where the data is mapped to.

Let me know if this works.

[zarzycki]

Thanks, @philipc2! Glad to see you still come around!

I tried this:

data_filenames = "./aorc/PrecipitationAORC_2004*.nc"
xrds = xr.open_mfdataset(data_filenames, combine='by_coords')
xrds = xrds.sel(time=slice(t0, t1))
uxgrid = ux.Grid.from_structured(lon=xrds.longitude, lat=xrds.latitude)
xrds = xrds.stack(n_face=('latitude', 'longitude'))
xrds = xrds.reset_index('n_face', drop=True)
uxds = ux.UxDataset.from_xarray(xrds, uxgrid=uxgrid)
print(uxgrid)
print("--")
print(xrds)
print("--")
print(uxds)

which returns:

<uxarray.Grid>
Original Grid Type: structured
Grid Dimensions:
  * n_node: 35305204
  * n_face: 35292601
  * n_max_face_nodes: 4
Grid Coordinates (Spherical):
  * node_lon: (35305204,)
  * node_lat: (35305204,)
Grid Coordinates (Cartesian):
Grid Connectivity Variables:
  * face_node_connectivity: (35292601, 4)
Grid Descriptor Variables:

--
<xarray.Dataset> Size: 34GB
Dimensions:       (time: 120, n_face: 35292601)
Coordinates:
  * time          (time) datetime64[ns] 960B 2004-09-04T12:00:00 ... 2004-09-...
Dimensions without coordinates: n_face
Data variables:
    APCP_surface  (time, n_face) float64 34GB dask.array<chunksize=(1, 151218), meta=np.ndarray>
Attributes:
    history:  Fri Feb 28 06:22:42 2025: ncks -4 -L 1 --cnk_dmn time,1 --cnk_d...
    NCO:      netCDF Operators version 5.2.6 (Homepage = http://nco.sf.net/, C...
--
<xarray.UxDataset> Size: 34GB
Dimensions:       (time: 120, n_face: 35292601)
Coordinates:
  * time          (time) datetime64[ns] 960B 2004-09-04T12:00:00 ... 2004-09-...
Dimensions without coordinates: n_face
Data variables:
    APCP_surface  (time, n_face) float64 34GB dask.array<chunksize=(1, 151218), meta=np.ndarray>

So it still seems from_xarray cannot quite link the uxgrid and the data together.

[rajeeja]

@zarzycki do you have files to reproduce this? I can take a look, not sure if it is fixed already..

[zarzycki]

Hi @rajeeja.

Seems like this issue still exists. Here is a repro code that should work on Derecho/Casper. I am currently just using the vanilla npl stack.

import uxarray as ux
import xarray as xr

print(f"UXarray version: {ux.__version__}")
print(f"Xarray version: {xr.__version__}")

# Open structured data
data_filenames = "/glade/work/zarzycki/uxarray/appendfaces/PrecipitationAORC_20050828*.nc"
xrds = xr.open_mfdataset(data_filenames)

# Use lat/lon to generate ux Grid object
uxgrid = ux.Grid.from_structured(lon=xrds.longitude, lat=xrds.latitude)

# Has face areas!
print(uxgrid.face_areas)

# Now convert the xr dataset to ux dataset using the uxgrid to translate
uxds = ux.UxDataset.from_xarray(xrds, uxgrid=uxgrid)
print(uxds)

# However, trying to apply subset.bounding_circle doesn't work
# because face_areas didn't carry over to uxds
center_coord = [-87.6298, 41.8781]
r = 2
radial_subset = uxds["APCP_surface"][0].subset.bounding_circle(center_coord, r)

returns:

(npl) zarzycki@crlogin3:/glade/work/zarzycki/uxarray/appendfaces> python repro_appendfaces.py 
UXarray version: 2025.6.0
Xarray version: 2025.7.1
<xarray.DataArray 'face_areas' (n_face: 35292601)> Size: 282MB
array([1.98766557e-08, 1.98766557e-08, 1.98766557e-08, ...,
       1.21328814e-08, 1.21328814e-08, 1.21328814e-08])
Dimensions without coordinates: n_face
Attributes:
    cf_role:    face_areas
    long_name:  Area of each face.
<xarray.UxDataset> Size: 7GB
Dimensions:       (time: 24, latitude: 4201, longitude: 8401)
Coordinates:
  * latitude      (latitude) float64 34kB 20.0 20.01 20.02 ... 54.98 54.99 55.0
  * longitude     (longitude) float64 67kB -130.0 -130.0 -130.0 ... -60.01 -60.0
  * time          (time) datetime64[ns] 192B 2005-08-28 ... 2005-08-28T23:00:00
Data variables:
    APCP_surface  (time, latitude, longitude) float64 7GB dask.array<chunksize=(1, 400, 400), meta=np.ndarray>
Traceback (most recent call last):
  File "/glade/work/zarzycki/uxarray/appendfaces/repro_appendfaces.py", line 25, in <module>
    radial_subset = uxds["APCP_surface"][0].subset.bounding_circle(center_coord, r)
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/apps/opt/conda/envs/npl/lib/python3.12/site-packages/uxarray/subset/dataarray_accessor.py", line 96, in bounding_circle
    return self.uxda._slice_from_grid(grid)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/apps/opt/conda/envs/npl/lib/python3.12/site-packages/uxarray/core/dataarray.py", line 1378, in _slice_from_grid
    raise ValueError(
ValueError: Data variable must be either node, edge, or face centered.

[rajeeja]

@zarzycki Reproduced and fixed — PR #1569.

Root cause was in the from_structured(lon=, lat=) path: it recorded no source-dim mapping and used a lowercase "structured" spec, so UxDataset.from_xarray never flattened the (lon, lat) data variables onto n_face, leaving them unmapped and breaking subset. Your from_structured(lon=xrds.longitude, lat=xrds.latitude) → from_xarray(xrds, uxgrid=uxgrid) workflow now maps air/APCP_surface onto n_face and subsets as expected.

The PR adds an end-to-end regression test (build grid from lon/lat, attach the Dataset, then bounding_circle) and a new section in the structured-grids user guide demonstrating exactly this interface-layer pattern.