Hack coordinate names when opening grid file to avoid conflict

There is a 'theta' variable in grid files, which prevents naming a
dimension 'theta'.

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

xbout/geometries.py

@@ -400,7 +400,7 @@ def add_toroidal_geometry_coords(ds, *, coordinates=None, grid=None):
     ds = ds.rename(y=coordinates["y"])
 
     # TODO automatically make this coordinate 1D in simplified cases?
-    ds = ds.rename(psixy=coordinates["x"])
+    ds[coordinates["x"]] = ds["psixy"]
     ds = ds.set_coords(coordinates["x"])
     ds[coordinates["x"]].attrs["units"] = "Wb"
 

xbout/load.py

@@ -367,7 +367,12 @@ def attrs_remove_section(obj, section):
                 ds.metadata[v] = grid[v].values
 
     # Update coordinates to match particular geometry of grid
-    ds = geometries.apply_geometry(ds, geometry, grid=grid)
+    if input_type == "grid":
+        # Specify coordinate names to avoid name clash
+        coordinates = {"x": "psi_poloidal", "y": "y", "z": "zeta"}
+    else:
+        coordinates = None
+    ds = geometries.apply_geometry(ds, geometry, grid=grid, coordinates=coordinates)
 
     if remove_yboundaries:
         ds = ds.bout.remove_yboundaries()