Saving subset ugrid to disk

[Chrismarsh]

What is the approach to save a subset ugrid to disk? For example, suppose I do

# massive dataset
ds = ux.open_mfdataset('output.2017.partial/output.2017.nc',
                       'output.2017.partial/output.2017.nc', 
                       parallel=False, lock=False)

# subset part of it
elev = ds.Elevation.subset.bounding_box(lon, lat)

# save subset
elev.to_xarray().to_netcdf('test2.nc', engine="netcdf4")

The test2.nc file does not contain valid ugrid topology:

$ ncdump -h test2.nc
netcdf test2 {
dimensions:
        n_face = 142546 ;
variables:
        double Mesh2_face_x(n_face) ;
                Mesh2_face_x:_FillValue = NaN ;
                Mesh2_face_x:standard_name = "latitude" ;
                Mesh2_face_x:long_name = "Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate)." ;
                Mesh2_face_x:units = "degrees_north" ;
        double Mesh2_face_y(n_face) ;
                Mesh2_face_y:_FillValue = NaN ;
                Mesh2_face_y:standard_name = "latitude" ;
                Mesh2_face_y:long_name = "Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate)." ;
                Mesh2_face_y:units = "degrees_north" ;
        double face_lon(n_face) ;
                face_lon:_FillValue = NaN ;
                face_lon:standard_name = "latitude" ;
                face_lon:long_name = "Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate)." ;
                face_lon:units = "degrees_north" ;
        double face_lat(n_face) ;
                face_lat:_FillValue = NaN ;
                face_lat:standard_name = "latitude" ;
                face_lat:long_name = "Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate)." ;
                face_lat:units = "degrees_north" ;
        double Elevation(n_face) ;
                Elevation:_FillValue = NaN ;
                Elevation:_QuantizeBitGroomNumberOfSignificantDigits = 4 ;
                Elevation:mesh = "Mesh2" ;
                Elevation:location = "face" ;
                Elevation:coordinates = "Mesh2_face_x Mesh2_face_y face_lat face_lon" ;
}

[philipc2]

Hi @Chrismarsh

UXarray currently operates with the assumption that Grid and Data files are separated. This means that the UGRID topology is typically stored in one file, and any other data variables are stored separately.

uxds = ...

subset = uxds['var'].subset.bounding_box(lon, lat)

# 1) Save the data
subset.to_xarray().to_netcdf("data.nc")

# 2) Save the grid
subset.uxgrid.to_xarray().to_netcdf("grid.nc")

I'm going to take a look into this and we can include a usage notebook better detailing this.

[Chrismarsh]

I've let this elev.uxgrid.to_xarray().to_netcdf('grid.nc', engine="netcdf4") run for 30 min and eventually killed it off. The other to_netcdf took <1s to run.

[Chrismarsh]

Tried again to day and the subset.uxgrid call fails deep in xarray with

File src/netCDF4/_netCDF4.pyx:5079, in netCDF4._netCDF4.Variable.__getitem__()

File src/netCDF4/_netCDF4.pyx:6051, in netCDF4._netCDF4.Variable._get()

File src/netCDF4/_netCDF4.pyx:2164, in netCDF4._netCDF4._ensure_nc_success()

RuntimeError: NetCDF: HDF error

which I can't say I've seen before!

[rajeeja]

I

What is the approach to save a subset ugrid to disk? For example, suppose I do

# massive dataset
ds = ux.open_mfdataset('output.2017.partial/output.2017.nc',
                       'output.2017.partial/output.2017.nc', 
                       parallel=False, lock=False)

# subset part of it
elev = ds.Elevation.subset.bounding_box(lon, lat)

# save subset
elev.to_xarray().to_netcdf('test2.nc', engine="netcdf4")

With the current implmentation, I thinkto_xarray call gets rid of uxgrid object, return type is a pure xarray object, so what you are seeing is expected. Try elev.to_netcdf(). I have WIP PR, where with to_xarray we can pass in a grid format to save the grid as an expected file format.

[philipc2]

+1

[Chrismarsh]

I have also tried elev.to_netcdf('test3.nc', engine="netcdf4") and it still has the same missing topology

 ncdump -h test3.nc                                                                                                                                                                                                                                
netcdf test3 {
dimensions:
        n_face = 142546 ;
variables:
        double Mesh2_face_x(n_face) ;
                Mesh2_face_x:_FillValue = NaN ;
                Mesh2_face_x:standard_name = "latitude" ;
                Mesh2_face_x:long_name = "Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate)." ;
                Mesh2_face_x:units = "degrees_north" ;
        double Mesh2_face_y(n_face) ;
                Mesh2_face_y:_FillValue = NaN ;
                Mesh2_face_y:standard_name = "latitude" ;
                Mesh2_face_y:long_name = "Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate)." ;
                Mesh2_face_y:units = "degrees_north" ;
        double face_lon(n_face) ;
                face_lon:_FillValue = NaN ;
                face_lon:standard_name = "latitude" ;
                face_lon:long_name = "Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate)." ;
                face_lon:units = "degrees_north" ;
        double face_lat(n_face) ;
                face_lat:_FillValue = NaN ;
                face_lat:standard_name = "latitude" ;
                face_lat:long_name = "Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate)." ;
                face_lat:units = "degrees_north" ;
        double Elevation(n_face) ;
                Elevation:_FillValue = NaN ;
                Elevation:_QuantizeBitGroomNumberOfSignificantDigits = 4 ;
                Elevation:mesh = "Mesh2" ;
                Elevation:location = "face" ;
                Elevation:coordinates = "Mesh2_face_x Mesh2_face_y face_lat face_lon" ;
}

where elev looks like

[rajeeja]

I have also tried elev.to_netcdf('test3.nc', engine="netcdf4") and it still has the same missing topology

can you try the latest main? also share something to reproduce this bug on our end.

[Chrismarsh]

@rajeeja @philipc2 replying here so I don't have to duplicate replies in two different threads.

• @rajeeja I built the latest main ux.__version__ = '2025.6.1.dev77+g743f0673' against xr.__version__='2025.7.1' and dask.__version__ = '2025.7.0'
• I am just trying to load data via ux, and then write it out to disk. Skipping the subset to keep it simple.

If I do

uxds = ux.open_mfdataset('output_2017/output_2017.nc',
                       'output_2017/output_2017.nc')
# 1) Save the data
uxds['t'].to_xarray().to_netcdf("t.nc")

I can then load this data against my original netcdf that contains the ugrid topology and it works correctly. ✅

ds= ux.open_mfdataset('output_2017/output_2017.nc',
                       't.nc')

However, if I do:

uxds.uxgrid.to_xarray().to_netcdf("grid.nc")

to write the ugrid topology, it just spins inside to_netcdf to acquire a write lock. I know I have had issues with netcdf before, e.g., pydata/xarray#3961, but given the individual data variable write works, this seems weird.

Poking around trying a few things, I thought I would compute uxds.uxgrid.to_xarray().compute() and write that to disk. In doing so , I realized that the uxgrid contains the time variable! I understand that grid topology can be time variant, but my original file does not have a time dimension to the ugrid topology so I don't understand where this is coming from.

uxarray.Grid

Dimensions:
n_face: 28495325n_max_face_nodes: 3n_node: 14638167time: 8737

however I can't select out a single time dim

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[16], [line 1](vscode-notebook-cell:?execution_count=16&line=1)
----> [1](vscode-notebook-cell:?execution_count=16&line=1) ds.uxgrid.isel(time=0)

File spack/environments/analysis/.spack-env/view/lib/python3.12/site-packages/uxarray/grid/grid.py:2268, in Grid.isel(self, inverse_indices, **dim_kwargs)
   2263     return _slice_face_indices(
   2264         self, dim_kwargs["n_face"], inverse_indices=inverse_indices
   2265     )
   2267 else:
-> [2268](spack/environments/analysis/.spack-env/view/lib/python3.12/site-packages/uxarray/grid/grid.py:2268)     raise ValueError(
   2269         "Indexing must be along a grid dimension: ('n_node', 'n_edge', 'n_face')"
   2270     )

ValueError: Indexing must be along a grid dimension: ('n_node', 'n_edge', 'n_face')

But time is totally there because .compute() tried to allocated 1.8TiB! ds.isel(time=0).uxgrid produces a null output.

[philipc2]

I'm happy to hear it worked. We'll try to improve our Grid reader. Thank you for all the feedback! Converting this to an issue.

[philipc2]

I think another option would be to drop any variables that contain the time dimension too, which might be a bit cleaner .

[philipc2]

ds = xr.open_mfdataset('output_2017/output_2017.nc').drop_dims("time")

[Chrismarsh]

Yes that's much better thanks.

[Chrismarsh]

It looks like grid write is missing the following
Mesh2:face_coordinates = "Mesh2_face_x Mesh2_face_y" ;
which is required for dual mesh, e.g., ESMF regridding. My original mesh had this attribute which seems to have been lost somewhere along the line