Merge pull request #217 from boutproject/document-workaround-load-conflicts

Document a workaround for conflicts in open_boutdataset(); some tidying

Author: johnomotani

xbout/boutdataarray.py

@@ -189,7 +189,7 @@ def from_field_aligned(self):
         return result
 
     @property
-    def regions(self):
+    def _regions(self):
         if "regions" not in self.data.attrs:
             raise ValueError(
                 "Called a method requiring regions, but these have not been created. "
@@ -279,7 +279,7 @@ def interpolate_parallel(
                 self.interpolate_parallel(
                     region, n=n, toroidal_points=toroidal_points, method=method
                 ).bout.to_dataset()
-                for region in self.data.regions
+                for region in self._regions
             ]
 
             # 'region' is not the same for all parts, and should not exist in the result,
@@ -301,7 +301,7 @@ def interpolate_parallel(
 
         # Select a particular 'region' and interpolate to higher parallel resolution
         da = self.data
-        region = da.regions[region]
+        region = da.bout._regions[region]
         tcoord = da.metadata["bout_tdim"]
         xcoord = da.metadata["bout_xdim"]
         ycoord = da.metadata["bout_ydim"]
@@ -409,9 +409,9 @@ def remove_yboundaries(self, return_dataset=False, remove_extra_upper=False):
 
         ycoord = self.data.metadata["bout_ydim"]
         parts = []
-        for region in self.data.regions:
+        for region in self._regions:
             part = self.data.bout.from_region(region, with_guards=0)
-            part_region = list(part.regions.values())[0]
+            part_region = list(part.bout._regions.values())[0]
             if part_region.connection_lower_y is None:
                 part = part.isel({ycoord: slice(myg, None)})
             if part_region.connection_upper_y is None:
@@ -525,7 +525,7 @@ def ddy(self, region=None):
         if region is None:
             # Call the single-region version of this method for each region, and combine
             # the results together
-            parts = [self.ddy(r).to_dataset() for r in self.data.regions]
+            parts = [self.ddy(r).to_dataset() for r in self._regions]
 
             # 'region' is not the same for all parts, and should not exist in the
             # result, so delete before merging
@@ -537,7 +537,7 @@ def ddy(self, region=None):
             result = xr.combine_by_coords(parts)[f"d({name})/dy"]
 
             # regions get mixed up during the split and combine_by_coords, so reset them
-            result.attrs["regions"] = self.data.regions
+            result.attrs["regions"] = self._regions
 
             return result
 
@@ -568,7 +568,7 @@ def ddy(self, region=None):
         result = (da.shift({ycoord: -1}) - da.shift({ycoord: 1})) / (2.0 * dy)
 
         # Remove any y-guard cells
-        region_object = da.regions[region]
+        region_object = da.bout._regions[region]
         if region_object.connection_lower_y is None:
             ylower = None
         else:
@@ -1000,13 +1000,26 @@ def interpolate_from_unstructured(
 
     # BOUT-specific plotting functionality: methods that plot on a poloidal (R-Z) plane
     def contour(self, ax=None, **kwargs):
+        """
+        Contour-plot a radial-poloidal slice on the R-Z plane
+        """
         return plotfuncs.plot2d_wrapper(self.data, xr.plot.contour, ax=ax, **kwargs)
 
     def contourf(self, ax=None, **kwargs):
+        """
+        Filled-contour-plot a radial-poloidal slice on the R-Z plane
+        """
         return plotfuncs.plot2d_wrapper(self.data, xr.plot.contourf, ax=ax, **kwargs)
 
     def pcolormesh(self, ax=None, **kwargs):
+        """
+        Colour-plot a radial-poloidal slice on the R-Z plane
+        """
         return plotfuncs.plot2d_wrapper(self.data, xr.plot.pcolormesh, ax=ax, **kwargs)
 
-    def regions(self, ax=None, **kwargs):
-        return plotfuncs.regions(self.data, ax=ax, **kwargs)
+    def plot_regions(self, ax=None, **kwargs):
+        """
+        Plot the regions into which xBOUT splits radial-poloidal arrays to handle
+        tokamak topology.
+        """
+        return plotfuncs.plot_regions(self.data, ax=ax, **kwargs)

xbout/boutdataset.py

@@ -143,7 +143,7 @@ def from_region(self, name, with_guards=None):
         return _from_region(self.data, name, with_guards)
 
     @property
-    def regions(self):
+    def _regions(self):
         if "regions" not in self.data.attrs:
             raise ValueError(
                 "Called a method requiring regions, but these have not been created. "

xbout/load.py

@@ -78,6 +78,21 @@ def open_boutdataset(
     `run_name` are ignored. `geometry` is treated specially, and can be passed when
     reloading a Dataset (along with `gridfilepath` if needed).
 
+    Troubleshooting
+    ---------------
+    Variable conflicts: sometimes, for example when loading data from multiple restarts,
+    some variables may have conflicts (e.g. a source term was changed between some of
+    the restarts, but the source term is saved as time-independent, without a
+    t-dimension). In this case one workaround is to pass a list of variable names to the
+    keyword argument `drop_vars` to ignore the variables with conflicts, e.g. if `"S1"`
+    and `"S2"` have conflicts
+    ```
+    ds = open_boutdataset("data*/boutdata.nc", drop_vars=["S1", "S2"])
+    ```
+    will open a Dataset which is missing `"S1"` and `"S2"`.\
+    [`drop_vars` is an argument of `xarray.open_dataset()` that is passed down through
+    `kwargs`.]
+
     Parameters
     ----------
     datapath : str or (list or tuple of xr.Dataset), optional
@@ -124,7 +139,7 @@ def open_boutdataset(
     info : bool or "terse", optional
     kwargs : optional
         Keyword arguments are passed down to `xarray.open_mfdataset`, which in
-        turn extra kwargs down to `xarray.open_dataset`.
+        turn passes extra kwargs down to `xarray.open_dataset`.
 
     Returns
     -------
@@ -327,9 +342,6 @@ def collect(
     info=True,
     prefix="BOUT.dmp",
 ):
-
-    from os.path import join
-
     """
 
     Extract the data pertaining to a specified variable in a BOUT++ data set
@@ -369,6 +381,7 @@ def collect(
     ds : numpy.ndarray
 
     """
+    from os.path import join
 
     datapath = join(path, prefix + "*.nc")
 

xbout/plotting/plotfuncs.py

@@ -24,7 +24,7 @@
     matplotlib.rcParams["pcolor.shading"] = "auto"
 
 
-def regions(da, ax=None, **kwargs):
+def plot_regions(da, ax=None, **kwargs):
     """
     Plots each logical plotting region as a different color for debugging.
 

xbout/plotting/utils.py

@@ -64,7 +64,10 @@ def plot_separatrix(da, sep_pos, ax, radial_coord="x"):
 
 def _decompose_regions(da):
 
-    return {region: da.bout.from_region(region, with_guards=1) for region in da.regions}
+    return {
+        region: da.bout.from_region(region, with_guards=1)
+        for region in da.bout._regions
+    }
 
 
 def _is_core_only(da):
@@ -90,7 +93,7 @@ def plot_separatrices(da, ax, *, x="R", y="Z"):
     ycoord = da0.metadata["bout_ydim"]
 
     for da_region in da_regions.values():
-        inner = list(da_region.regions.values())[0].connection_inner_x
+        inner = list(da_region.bout._regions.values())[0].connection_inner_x
         if inner in da_regions:
             da_inner = da_regions[inner]
 
@@ -136,14 +139,14 @@ def plot_targets(da, ax, *, x="R", y="Z", hatching=True):
         y_boundary_guards = 0
 
     for da_region in da_regions.values():
-        if list(da_region.regions.values())[0].connection_lower_y is None:
+        if list(da_region.bout._regions.values())[0].connection_lower_y is None:
             # lower target exists
             x_target = da_region.coords[x].isel(**{ycoord: y_boundary_guards})
             y_target = da_region.coords[y].isel(**{ycoord: y_boundary_guards})
             [line] = ax.plot(x_target, y_target, "k-", linewidth=2)
             if hatching:
                 _add_hatching(line, ax)
-        if list(da_region.regions.values())[0].connection_upper_y is None:
+        if list(da_region.bout._regions.values())[0].connection_upper_y is None:
             # upper target exists
             x_target = da_region.coords[x].isel(**{ycoord: -y_boundary_guards - 1})
             y_target = da_region.coords[y].isel(**{ycoord: -y_boundary_guards - 1})

xbout/region.py

@@ -1221,11 +1221,11 @@ def _concat_inner_guards(da, da_global, mxg):
     if mxg <= 0:
         return da
 
-    if len(da.regions) > 1:
+    if len(da.bout._regions) > 1:
         raise ValueError("da passed should have only one region")
-    region = list(da.regions.values())[0]
+    region = list(da.bout._regions.values())[0]
 
-    if region.connection_inner_x not in da_global.regions:
+    if region.connection_inner_x not in da_global.bout._regions:
         # No connection, or plotting restricted set of regions not including this
         # connection
         return da
@@ -1236,9 +1236,9 @@ def _concat_inner_guards(da, da_global, mxg):
     ycoord = da_global.metadata["bout_ydim"]
 
     da_inner = da_global.bout.from_region(region.connection_inner_x, with_guards=0)
-    if len(da_inner.regions) > 1:
+    if len(da_inner.bout._regions) > 1:
         raise ValueError("da_inner should have only one region")
-    region_inner = list(da_inner.regions.values())[0]
+    region_inner = list(da_inner.bout._regions.values())[0]
 
     if (
         myg_da > 0
@@ -1276,7 +1276,7 @@ def _concat_inner_guards(da, da_global, mxg):
         da_inner_lower = da_inner_lower.isel(
             **{xcoord: slice(-mxg, None), ycoord: slice(-myg_da, None)}
         )
-        save_regions = da_inner.regions
+        save_regions = da_inner.bout._regions
         da_inner = xr.concat((da_inner_lower, da_inner), ycoord, join="exact")
         # xr.concat takes attributes from the first variable, but we need da_inner's
         # regions
@@ -1314,7 +1314,7 @@ def _concat_inner_guards(da, da_global, mxg):
         da_inner[xcoord].data[...] = new_xcoord.data
         da_inner[ycoord].data[...] = new_ycoord.data
 
-    save_regions = da.regions
+    save_regions = da.bout._regions
     da = xr.concat((da_inner, da), xcoord, join="exact")
     # xr.concat takes attributes from the first variable (for xarray>=0.15.0, keeps attrs
     # that are the same in all objects for xarray<0.15.0)
@@ -1332,11 +1332,11 @@ def _concat_outer_guards(da, da_global, mxg):
     if mxg <= 0:
         return da
 
-    if len(da.regions) > 1:
+    if len(da.bout._regions) > 1:
         raise ValueError("da passed should have only one region")
-    region = list(da.regions.values())[0]
+    region = list(da.bout._regions.values())[0]
 
-    if region.connection_outer_x not in da_global.regions:
+    if region.connection_outer_x not in da_global.bout._regions:
         # No connection, or plotting restricted set of regions not including this
         # connection
         return da
@@ -1347,9 +1347,9 @@ def _concat_outer_guards(da, da_global, mxg):
     ycoord = da_global.metadata["bout_ydim"]
 
     da_outer = da_global.bout.from_region(region.connection_outer_x, with_guards=0)
-    if len(da_outer.regions) > 1:
+    if len(da_outer.bout._regions) > 1:
         raise ValueError("da_outer should have only one region")
-    region_outer = list(da_outer.regions.values())[0]
+    region_outer = list(da_outer.bout._regions.values())[0]
 
     if (
         myg_da > 0
@@ -1387,7 +1387,7 @@ def _concat_outer_guards(da, da_global, mxg):
         da_outer_lower = da_outer_lower.isel(
             **{xcoord: slice(-mxg, None), ycoord: slice(-myg_da, None)}
         )
-        save_regions = da_outer.regions
+        save_regions = da_outer.bout._regions
         da_outer = xr.concat((da_outer_lower, da_outer), ycoord, join="exact")
         # xr.concat takes attributes from the first variable, but we need da_outer's
         # regions
@@ -1425,7 +1425,7 @@ def _concat_outer_guards(da, da_global, mxg):
         da_outer[xcoord].data[...] = new_xcoord.data
         da_outer[ycoord].data[...] = new_ycoord.data
 
-    save_regions = da.regions
+    save_regions = da.bout._regions
     da = xr.concat((da, da_outer), xcoord, join="exact")
     # xarray<0.15.0 only keeps attrs that are the same on all variables passed to concat
     da.attrs["regions"] = save_regions
@@ -1442,11 +1442,11 @@ def _concat_lower_guards(da, da_global, mxg, myg):
     if myg <= 0:
         return da
 
-    if len(da.regions) > 1:
+    if len(da.bout._regions) > 1:
         raise ValueError("da passed should have only one region")
-    region = list(da.regions.values())[0]
+    region = list(da.bout._regions.values())[0]
 
-    if region.connection_lower_y not in da_global.regions:
+    if region.connection_lower_y not in da_global.bout._regions:
         # No connection, or plotting restricted set of regions not including this
         # connection
         return da
@@ -1525,7 +1525,7 @@ def _concat_lower_guards(da, da_global, mxg, myg):
         da_lower[xcoord].data[...] = new_xcoord.data
         da_lower[ycoord].data[...] = new_ycoord.data
 
-    save_regions = da.regions
+    save_regions = da.bout._regions
     da = xr.concat((da_lower, da), ycoord, join="exact")
     # xr.concat takes attributes from the first variable (for xarray>=0.15.0, keeps attrs
     # that are the same in all objects for xarray<0.15.0)
@@ -1543,11 +1543,11 @@ def _concat_upper_guards(da, da_global, mxg, myg):
     if myg <= 0:
         return da
 
-    if len(da.regions) > 1:
+    if len(da.bout._regions) > 1:
         raise ValueError("da passed should have only one region")
-    region = list(da.regions.values())[0]
+    region = list(da.bout._regions.values())[0]
 
-    if region.connection_upper_y not in da_global.regions:
+    if region.connection_upper_y not in da_global.bout._regions:
         # No connection, or plotting restricted set of regions not including this
         # connection
         return da
@@ -1625,7 +1625,7 @@ def _concat_upper_guards(da, da_global, mxg, myg):
         da_upper[xcoord].data[...] = new_xcoord.data
         da_upper[ycoord].data[...] = new_ycoord.data
 
-    save_regions = da.regions
+    save_regions = da.bout._regions
     da = xr.concat((da, da_upper), ycoord, join="exact")
     # xarray<0.15.0 only keeps attrs that are the same on all variables passed to concat
     da.attrs["regions"] = save_regions
@@ -1637,7 +1637,7 @@ def _from_region(ds_or_da, name, with_guards):
     # ensure we do not modify the input
     ds_or_da = ds_or_da.copy(deep=True)
 
-    region = ds_or_da.regions[name]
+    region = ds_or_da.bout._regions[name]
     xcoord = ds_or_da.metadata["bout_xdim"]
     ycoord = ds_or_da.metadata["bout_ydim"]
 
@@ -1669,7 +1669,7 @@ def _from_region(ds_or_da, name, with_guards):
 
     # If the result (which only has a single region) is passed to from_region a
     # second time, don't want to slice anything.
-    single_region = list(result.regions.values())[0]
+    single_region = list(result.bout._regions.values())[0]
     single_region.xinner_ind = None
     single_region.xouter_ind = None
     single_region.ylower_ind = None