Filtering when reading a Grid file, avoid time leaking into the Grid

[dylannelson]

Closes#1444 (and others (?) planned to update)

Overview

1. time was shown to exist in a grid from a user reading in a grid and data from a single file
2. Later when the user tried to use .subset there was a crash related to time being different on the grid vs the data
3. time shouldn't realistically be on the grid to start with, so this was a concern, and was causing downstream issues not addressed when reading in the data
4. There was 2 main ideas of where to place filtering to avoid this, Grid.__init__ or _read_ugrid, both were tested extensively
   • Grid.__init__ is run in multiple scenarios, including when a file is read from disk, but in other scenarios like _slice_from_grid which can occur when subsetting
   • the path to _read_ugrid occurs primarily when a file is read from disk, like:
     • ux.open_grid(..) -> Grid.from_dataset(..) -> _read_ugrid(ds)
     • _read_ugrid(ds) main function is: "Parses an unstructured grid dataset and encodes it in the UGRID" so this seems like the right place for cleaning
5. This has been tested in multiple ways and is shown to:
   • Scalar time gone from grid coords and grid variables
   • node_lon, node_lat, face_node_connectivity preserved
   • Grid non-empty (n_node, n_face > 0)
   • bounding_box runs, no error
   • Data time preserved
   • Subset returned faces
   • Sliced grid keeps subgrid_face_indices
   • sel(time=...) still works
   • Dimensional time axis + coord stripped
   • Leaked data variable stripped

Expected Usage

Changes are all located in uxarray\uxarray\io\_ugrid.py

Test these changes on an old version and the new version. This test likely fails in old uxarray and succeeds in this branch

importpandasaspd, xarrayasxr, uxarrayasux# This file exists in uxarray repo, update with your local path added if neededpath_to_nc="uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc"# Opens a grid, adds scalar time, saves to local
(xr.open_dataset(path_to_nc)
.assign_coords(time=pd.Timestamp("2016-10-01"))
.to_netcdf("grid_with_time.nc"))
uxgrid=ux.open_grid("grid_with_time.nc")
print("grid coords:", list(uxgrid._ds.coords))
print("time removed?", "time"notinuxgrid._ds.coords) # True = fix working

PR Checklist

General

• An issue is created and linked
• Added appropriate labels (if your uxarray repo permissions allow it)
• Filled out Overview and Expected Usage (if applicable) sections

Testing & Benchmarking

• Adequate tests are created if there is new functionality
• Tests are not too basic (such as simply calling a function and nothing else)
• Tests cover all major paths in your new functions
• If this PR could affect performance, ran ASV benchmarks and confirmed they show expected behavior (add a new benchmark if necessary)

Documentation

• Docstrings have been added to all new functions
• Docstrings have been updated with any function changes
• User (public) functions have been added to docs/api.rst
• Internal (private) function names start with an underscore (_)

Examples

• All notebook examples cleared the output of all cells before committing
• New notebook examples added to appropriate folder (gallery: docs/examples/; guide: docs/user-guide/; quickstart: docs/getting-started/)
• New notebook examples referenced in appropriate .rst file (gallery: docs/gallery.rst; guide: docs/userguide.rst; quickstart: docs/quickstart.rst)
• New notebook gallery examples added entry in docs/gallery.yml with appropriate thumbnail photo in docs/_static/thumbnails/

AI Disclosure

AI Usage:

• I take responsibility for all AI-generated content in my PR.
• I have tested all AI-generated content in my PR.

[erogluorhan]

Maybe too early for a review, but I'd like to share some thoughts:

1. time was shown to exist in a grid from a user reading in a grid and data from a single file

   A little more clarity could be helpful here. IIRC, time is not shown directly in the Grid object but in the Grid._ds.

2. That said, I'd like to know what actual purpose we had with that object attribute, i.e. _ds. If it was to keep complete track of the original file content through xarray.Dataset, I'd say this dropping at the level of _ds might break that.

3. Maybe we'd want to fix it at the Grid itself rather than altering _ds, but would it be possible?

4. Finally, ugrid (i.e. _ugrid.py) is only one of the several formats we support. What about other formats in our I/O, e.g. _mpas.py etc. Their grid files can come in a similar way. Maybe we will need to do this fix at the Grid level rather than particular I/O modules?

[dylannelson]

#1444 — Solution history

This is a bit about each version of the idea, why it's located where it is, and what the final code looks like.

---

Version 1 — .values guard in _slice_from_grid

# uxarray/core/dataarray.py — _slice_from_grid()n_face=sliced_grid._ds["subgrid_face_indices"].values# was: [...["subgrid_face_indices"]]

• Location: within _slice_from_grid in core/dataarray.py.
• Idea: Pass the face indexer as .values (a plain array) so xarray does positional indexing with no coordinate alignment. It fixes the crash for every subset method, but only guards the symptom — the stray time stays on the grid, so it never addresses why the grid carries a time.

---

Version 2 — Name-based clean build in _read_ugrid

# uxarray/io/_ugrid.py — end of _read_ugrid()ds=_keep_only_grid_vars(ds) # keep only recognized grid vars/coords, drop the rest

• Location: at the end of _read_ugrid in io/_ugrid.py.
• Idea: Filter the parsed dataset down to recognized grid variables, dropping extras where the grid is built. This removes the stray time at its source (the leaky reader), so nothing downstream can collide. It's essentially the final answer — a simpler-looking alternative was proposed next, which prompted the comparison below.

---

Version 3 — drop_dims in Grid.__init__ (team's idea)

# uxarray/grid/grid.py — Grid.__init__extra_dims=set(grid_ds.dims) -set(DIM_NAMES)
grid_ds=grid_ds.drop_dims(extra_dims)

• Location: at the top of Grid.__init__ in grid/grid.py.
• Idea: Drop any dimension not in DIM_NAMES, centrally for every reader, to remove stray non-grid data by axis. It captured a dimensionaltime, but the leaked time here is a scalar (0-d) coord with no dimension — so drop_dims finds nothing to remove and the bug survives.

---

Version 4 — "Best of both" (drop_dims + 0-d/name drop)

grid_ds=grid_ds.drop_dims(set(grid_ds.dims) -set(DIM_NAMES)) # dimensional extrasdrop= [nforn, vingrid_ds.variables.items() ifnnotinKEEPandv.ndim==0]
grid_ds=grid_ds.drop_vars(drop, errors="ignore") # scalar extras

• Location: at the top of Grid.__init__ in grid/grid.py.
• Idea: Keep V3's drop_dims and add a second pass dropping 0-d non-grid vars, catching the scalar time it missed. It works, but uses two mechanisms for what one can do — V5's name-based keep-list covers scalar, dimensional, and data-var leaks in a single pass.

---

Version 5 — Single name-based keep-list in _read_ugrid (final)

# uxarray/io/_ugrid.py — end of _read_ugrid()def_keep_only_grid_vars(ds):
keep= {"grid_topology", *SPHERICAL_COORD_NAMES, *CARTESIAN_COORD_NAMES,
*CONNECTIVITY_NAMES, *DESCRIPTOR_NAMES}
returnds.drop_vars([nforninds.variablesifnnotinkeep], errors="ignore")

• Location: at the end of _read_ugrid in io/_ugrid.py — different location (see next section).
• Idea: One name-based keep-list that drops every kind of non-grid extra (scalar, dimensional, data vars) in a single pass, collapsing V2–V4 into the simplest complete fix. Could be paired with the .values option, but didn't seem necessary.

---

Why V5 lives in _read_ugrid, not Grid.__init__

Two code paths build a Grid. They share Grid.__init__ but not _read_ugrid, so why is it in _read_ugrid?

Path A — opening a grid from a file (where the time leaks in):

ux.open_grid(...) core/api.py
└─ _open_dataset_with_fallback() # full xr.Dataset — stray `time` included
└─ Grid.from_dataset(...) grid.py
└─ no source_grid_spec
└─ _parse_grid_type()
└─ _read_ugrid(...) io/_ugrid.py ← renames in place, returns EVERYTHING (leak enters here)
└─ Grid.__init__(...) grid.py # self._ds = grid_ds

Path B — building a sliced grid during a subset (trying not to disturb this path):

uxda.subset.bounding_box(...)
└─ Grid.isel(...) grid.py
└─ _slice_face_indices(...) grid/slice.py
ds = grid._ds.isel(...) # copy of cleaned grid
ds["subgrid_face_indices"] = ...
└─ Grid.from_dataset(...) slice.py
└─ grid_ds = ds (_read_ugrid is skipped)
└─ Grid.__init__(...) grid.py # self._ds holds subgrid_*_indices

	Through _read_ugrid?	Through Grid.__init__?	_ds holds subgrid_*_indices?
Path A (file open)	Yes	Yes	No
Path B (subset)	No	Yes	Yes (needed by _slice_from_grid)

• Filter in _read_ugrid (V5): Path A only. time is dropped as the grid is read; Path B never enters _read_ugrid, so its subgrid_*_indices survive, and its sliced grid (copied from the already-clean grid._ds) comes out clean for free.
• Filter in Grid.__init__: both paths. It would strip subgrid_face_indices from Path B's sliced grid → _slice_from_grid hits a KeyError. Fixes the leak but breaks subsetting.

---

Final code preview

Basically came down to two edits, both in uxarray/io/_ugrid.py.

Edit 1 — add one import at the top of the file (ugrid is already imported; only DESCRIPTOR_NAMES is new):

importnumpyasnpimportxarrayasxrimportuxarray.conventions.ugridasugridfromuxarray.constantsimportINT_DTYPE, INT_FILL_VALUEfromuxarray.grid.connectivityimport_replace_fill_valuesfromuxarray.conventions.descriptorsimportDESCRIPTOR_NAMES# <---- this line

Edit 2 — filter at the end of _read_ugrid, then define the helper directly after it:

def_read_ugrid(ds):
"""Parses an unstructured grid dataset and encodes it in the UGRID conventions."""# ... (topology parse, coord/connectivity renames, dim swaps — unchanged) ...dim_dict[ds["face_node_connectivity"].dims[1]] =ugrid.N_MAX_FACE_NODES_DIMds=ds.swap_dims(dim_dict)
# ===== NEW =====# Strip non-grid extras (e.g. a stray scalar `time` coordinate, or unrelated data variables)ds=_keep_only_grid_vars(ds)
# ===== end NEW =====returnds, dim_dict# ===== NEW helper — place directly after _read_ugrid =====def_keep_only_grid_vars(ds):
"""Return ``ds`` with only recognized UGRID grid variables/coordinates. Anything else on the dataset (a stray scalar ``time`` coordinate, unrelated data variables, etc.) is dropped so it cannot leak onto ``grid._ds``. Runs on the file-read path only. """# uxarray's own canonical grid-variable names (the same lists Grid filters against)keep= {"grid_topology"}
keep.update(ugrid.SPHERICAL_COORD_NAMES) # node/edge/face lon-latkeep.update(ugrid.CARTESIAN_COORD_NAMES) # node/edge/face x-y-zkeep.update(ugrid.CONNECTIVITY_NAMES) # face_node_connectivity, edge_node_connectivity, ...keep.update(DESCRIPTOR_NAMES) # n_nodes_per_face, face_areas, boundary_*_indices, ...# drop_vars removes variables/coords by name (never bare dimensions), so grid# dims survive with their variables; errors="ignore" tolerates absent names.drop= [namefornameinds.variablesifnamenotinkeep]
returnds.drop_vars(drop, errors="ignore")
# ===== end NEW helper =====# ===== existing code =====def_encode_ugrid(ds):
"""Encodes an unstructured grid represented under a ``Grid`` object as a ``xr.Dataset`` with an updated grid topology variable."""if"grid_topology"inds:
ds=ds.drop_vars(["grid_topology"])
# ... (rest of _encode_ugrid, unchanged) ...

[dylannelson]

@erogluorhan

4. Finally, ugrid (i.e. `_ugrid.py`) is only one of the several formats we support. What about other formats in our I/O, e.g. _mpas.py etc. Their grid files can come in a similar way. Maybe we will need to do this fix at the Grid level rather than particular I/O modules?

Sorry saw this after I was typing all the above content.

Edit: Here's an update, more focused on # 4 to start

This is a great question and I wish I had realized sooner. Trying to find a solution that was very high level, like in Grid.__innit__ seemed like a it made the most sense but kept causing problems downstream. Looking to something less high up made it easier, but I totally missed that it would only work for certain file types

So I had to look at 2 questions:

1. Does this leak happen for other file types?
   • If no, we are good, and it was a ugrid specific issue
   • If yes, and it's every other file, then we may have to find some function at a higher level that can handle all the leaks the same
2. Is there somewhere up, that can solve this problem by interacting with all the different file types in the same way, without disturbing other functionality
   • I've looked at this before and was struggling to find a place

Looking at # 1, I performed a similar check like before, read in a demo dataset, add time, save, reopen, check if there's time, and see if it crashes when there's a subset. Here's what I got

Reader	Demo grid tested (test/meshfiles/…)	How grid vars are populated	Leaks?
MPAS	mpas/QU/mesh.QU.1920km.151026.nc	DataArrays derived from in_ds	yes
ESMF	esmf/ne30/ne30pg3.grid.nc	xr.DataArray(<DataArray>, …)	yes
ICON	icon/R02B04/icon_grid_0010_R02B04_G.nc	in_ds["…"].T - 1 (DataArray)	yes
Exodus	exodus/outCSne8/outCSne8.g	DataArrays from the source	yes
SCRIP	scrip/ne30pg2/grid.nc	xr.DataArray(<numpy>, …)	no
GEOS-CS	geos-cs/c12/test-c12.native.nc4	in_ds["…"].values.ravel()	no
UGRID	ugrid/quad-hexagon/grid.nc	renames the input dataset in place	no
Structured / FESOM2	structured/outCSne30_vortex_structured.nc (FESOM2 untested)	crashes when re-opening	yes but worse

So basically from what I'm seeing right now, we have 3 different outcomes

• it's fine as is (good, 3 file types)
• it's the same problem as ugrid (bad)
• it introduces a problem that crashes when reading (worse than what happened with ugrid)

I'll keep looking for a solution for # 2, but worried at it's impact, especially since every file type seems to have a unique way to handle the data when opening.

[erogluorhan]

The table in your latest comment showing no leaks for UGRID confused me a bit (From your other words, it sounds like UGRID has the time leaking). Could you clarify?

[dylannelson]

@erogluorhan
For # 1 and # 2 and # 3, I think from what I'm seeing, _ds is used for a lot of functionality, and the added time being in there causes it to leak out into anything that can call _ds (which is a lot, shown below). The "Grid" and "_ds" can't be thought of as separate, because they are so heavily intertwined.

"the Grid" is essentially _ds

Almost every Grid property is a thin read of self._ds. Straight from uxarray/grid/grid.py (trimmed to the returns):

@propertydefn_face(self): returnself._ds.sizes["n_face"] # :878@propertydefface_node_connectivity(self): returnself._ds["face_node_connectivity"] # :1261@propertydefattrs(self): returnself._ds.attrs# :836@propertydefcoordinates(self): # only the NAMES the conventions recognize:return {cforcinSPHERICAL_COORD_NAMESifcinself._ds} | \
{cforcinCARTESIAN_COORD_NAMESifcinself._ds} # :799@propertydefdims(self):
return {dfordinDIM_NAMESifdinself._ds.dims} # :787

Line references (pinned to 17c20f91):

• n_face — grid.py:878
• face_node_connectivity — grid.py:1261
• attrs — grid.py:836
• coordinates — grid.py:807
• dims — grid.py:791

The same pattern appears to be true for others as well, like: n_node, n_edge, node_lon, every *_connectivity,
face_areas, bounds, … — all return self._ds[...].

This can be a problem because time can leak onto other calls, like what we saw with subset. When the subset pulls _ds["subgrid_face_indices"] out as an index array, time rides along and collides with the data's time.

[dylannelson]

@erogluorhan

The table in your latest comment showing no leaks for UGRID confused me a bit (From your other words, it sounds like UGRID has the time leaking). Could you clarify?

It's no leak for UGRID in my current version, which fixed it. main would instead show that yes, there is a leak in UGRID

[dylannelson]

and for

fixing at the Gridwithout altering _ds,

that's effectively what I had originally, via .values. (Version 1 in the big comment above)

The issue with this is that it fixes a symptom, without addressing the cause. We discussed this early on and it was turned down when we originally discussed it. It's basically a patch in the .isel indexer (_slice_from_grid) so that when the data is on the way out, we strip any extra data that may cause a clash. This works for our subset incident, but it won't permanently remove the leaked time and it won't guarantee it won't cause issues in other functions,

[erogluorhan]

and for

fixing at the Gridwithout altering _ds,

that's effectively what I had originally, via .values. (Version 1 in the big comment above)

The issue with this is that it fixes a symptom, without addressing the cause. We discussed this early on and it was turned down when we originally discussed it. It's basically a patch in the .isel indexer (_slice_from_grid) so that when the data is on the way out, we strip any extra data that may cause a clash. This works for our subset incident, but it won't permanently remove the leaked time and it won't guarantee it won't cause issues in other functions,

You're right; "without altering _ds" doesn't seem feasible or ideal, but I think we should discuss if the fix should be in the individual I/O modules (like _mpas etc.) or right before the Grid object creation, i.e. after the specific grid format reading. Thougts?

[rajeeja]

Other readers like MPAS, ESMF, Exodus, ICON, and FESOM2 still leave time on _ds — MPAS reproduces the same IndexError from #1444. Could the filter move to the end of the auto-detect branch in Grid.from_dataset, before the return cls(...) at grid.py:323, so it covers all formats at once? The keep-list also drops bounds and max_face_radius when re-reading a uxarray-written grid.

[dylannelson]

This new push should have everything we have wanted from the last few discussions
Early on we had the idea of moving the solution to grid.py in Grid.__init__ and doing the filtering there. The filtering would look at the dimensions and remove any that didn't belong. Unfortunately this didn't work, while also introducing downstream issues. It looked like this:

extra_dims=set(grid_ds.dims) -set(DIM_NAMES)
grid_ds=grid_ds.drop_dims(extra_dims)

The time in #1444 is a scalar (0-d) coordinate — it has no dimension, so it never appears in grid_ds.dims, extra_dims comes back empty, and drop_dims doesn't change it. I moved on from looking in grid.py at this point but now I'm back to it, and it seems to be working.

This solution is different in that it is dropping non-grid coordinates by name, which catches the scalar one that we've been having trouble with. Here's what this solution looks like instead:

grid_coord_names=set(ugrid.SPHERICAL_COORD_NAMES) |set(ugrid.CARTESIAN_COORD_NAMES)
stray= [cforcinds.coordsifcnotingrid_coord_names]
returnds.drop_vars(stray, errors="ignore")

Testing

Before:

After:

Notebook Preview:
18-new-fix-demo.html

Downstream impact

Performed a few downstream tests with claude code to test what affects it may have, here's the results of what it looked at:

__init__ runs on every construction path — from_dataset (all formats), from_topology/from_healpix/etc., the subset's sliced grids, and Grid.copy(). Static sweep:

Area	Finding	Verdict
Grid properties (n_face, node_lon, …)	read by name; recognized coords kept	safe
Reads of stray coords off grid._ds	grep '_ds.coords' / _ds["time"] → none	safe
Grid↔data linking (open_dataset)	by dimensions (_map_dims_to_ugrid), not coords	safe
Subset bookkeeping	subgrid_face/edge/node_indices are data vars	safe
Grid.copy()	re-runs __init__ on a clean _ds	safe (idempotent)
Encoders / to_xarray	no .coords use	safe
Every .coords read in the package	on the data objects / shapely — never grid._ds	safe
Tests	none assert coords on grid._ds	safe