Fix loading of grid files

[johnomotani]

There is (now) a 'theta' variable in grid files, which prevents naming a dimension 'theta', so need to choose a different dimension name when opening a grid file.

Also keep the 'psixy' variable in the Dataset - create a new variable referecing it rather than renaming it to make the 'psi_poloidal' coordinate.

Along with https://github.com/boutproject/hypnotoad/pull/187, fixes #306.

The examples/plot_grid.py script should work again.

[johnomotani]

Test failures seem to be numpy/xarray incompatibility issues. Maybe would be fixed by #303?

[mrhardman]

I have tried checking out bout and #307, rebuilding, and running the scripts featured in #306 . I still see the error

Traceback (most recent call last):
  File "/*/lib/python3.13/site-packages/xarray/core/dataset.py", line 1317, in _construct_dataarray
    variable = self._variables[name]
               ~~~~~~~~~~~~~~~^^^^^^
KeyError: 't_array'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/*/plot_test.py", line 11, in <module>
    grid = open_boutdataset(gridfilepath, geometry="toroidal")
  File "/*/xBOUT/xbout/load.py", line 279, in open_boutdataset
    ds, remove_yboundaries = _auto_open_mfboutdataset(
                             ~~~~~~~~~~~~~~~~~~~~~~~~^
        datapath=datapath,
        ^^^^^^^^^^^^^^^^^^
    ...<4 lines>...
        **kwargs,
        ^^^^^^^^^
    )
    ^
  File "/*/xBOUT/xbout/load.py", line 741, in _auto_open_mfboutdataset
    _, unique_indices = unique(ds["t_array"], return_index=True)
                               ~~^^^^^^^^^^^
  File "/*/lib/python3.13/site-packages/xarray/core/dataset.py", line 1410, in __getitem__
    return self._construct_dataarray(key)
           ~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^
  File "/*/lib/python3.13/site-packages/xarray/core/dataset.py", line 1319, in _construct_dataarray
    _, name, variable = _get_virtual_variable(self._variables, name, self.dims)
                        ~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/*/lib/python3.13/site-packages/xarray/core/dataset.py", line 175, in _get_virtual_variable
    raise KeyError(key)
KeyError: 't_array'

What did I miss?

[johnomotani]

If you ncdump -h your grid file, does it have a 't' dimension? With the updated https://github.com/boutproject/hypnotoad/pull/187 it shouldn't.

[mrhardman]

If you ncdump -h your grid file, does it have a 't' dimension? With the updated boutproject/hypnotoad#187 it shouldn't.

Here's my check:

netcdf grid {
dimensions:
        t = UNLIMITED ; // (5 currently)

Which is odd, given that I have an editable hypnotoad install on your branch.

hypnotoad$ git log
commit b23d387b93d3a15107205a572866322ecf7c13d7 (HEAD -> leg-only-grid, origin/leg-only-grid)
Author: John Omotani <[email protected]>
Date:   Thu Nov 7 16:44:51 2024 +0000

    Make "single_region" option work for non-orthogonal grids

[johnomotani]

You need to pull the latest version, with the changes I added yesterday.

[mrhardman]

You need to pull the latest version, with the changes I added yesterday.

My mistake. Now I have the following:

/*/lib/python3.13/site-packages/xarray/core/concat.py:527: FutureWarning: unique with argument that is not not a Series, Index, ExtensionArray, or np.ndarray is deprecated and will raise in a future version.
  common_dims = tuple(pd.unique([d for v in vars for d in v.dims]))
Traceback (most recent call last):
  File "/home/mrhardman/hermes-3-work/hypnotoad_experiment/plot_test.py", line 18, in <module>
    grid["psi_poloidal"].bout.regions()
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AttributeError: 'BoutDataArrayAccessor' object has no attribute 'regions'

from snippet

grid = open_boutdataset(gridfilepath, geometry="toroidal")
grid["psi_poloidal"].bout.contourf()
grid["psi_poloidal"].bout.contour()
grid["psi_poloidal"].bout.pcolormesh()
grid["psi_poloidal"].bout.pcolormesh(shading="gouraud")
grid["psi_poloidal"].bout.regions()

[johnomotani]

Ah, that was just a mistake in the example plot script that I've forgotten to update. The line should have been

grid["psi_poloidal"].bout.plot_regions()

[mrhardman]

Plots are now generating from the example script on my local machine. Thanks for tracking down the issues!