Add BoundedRegion to keep track of boundary's child sections

This inherits the behavior from `GriddedRegion` but is instead created
from an instance of the new `Section` class in `sectionate` (see related
[commit](hdrake/sectionate@5bfce02)).

The benefit of this is that transports computed normal to the region's
boundary can then be decomposed into the constituent parts because
we retain information about the child sections.

regionate/region.py

@@ -1,12 +1,17 @@
 import numpy as np
+import xarray as xr
 
+import sectionate as sec
 from sectionate import is_section_counterclockwise
 from .utilities import *
 from .grid_conform import (
     get_region_boundary_grid_indices,
     mask_from_grid_boundaries
 )
 
+import os
+from pathlib import Path
+
 class Region:
     """
     A named polygonal region defined by a list or array of geographical coordinates.
@@ -101,7 +106,14 @@ def remove_duplicate_points(self, closeness_threshold=5.e3):
         closeness_threshold : float
             A short distance within which points are deemed to be identical. Default: 5.e3.
         """
-        self.lons, self.lats = unique_lonlat(self.lons, self.lats, closeness_threshold=closeness_threshold)
+        self.lons, self.lats = unique_lonlat(
+            self.lons,
+            self.lats,
+            closeness_threshold=closeness_threshold
+        )
+
+    def __repr__(self):
+        return f"""{str(type(self))[8:-2]}("{self.name}")"""
 
 class GriddedRegion(Region):
     """
@@ -115,7 +127,7 @@ def __init__(
         grid,
         positive_in=True,
         mask=None,
-        ij=None
+        ij=None,
         ):
         """
         Create a Region object (named `name`) from arrays of (`lons`, `lats`) and an ocean model `grid`. 
@@ -149,6 +161,7 @@ def __init__(
         >>> region = reg.Region("Bermuda Triangle", lons, lats, grid)
         """
         self.grid = grid
+        self.save = {}
 
         if len(lons)>=3 and len(lats)>=3 and ij is None:
             self.initiate_from_boundary(
@@ -186,7 +199,7 @@ def initiate_from_boundary(
         lons,
         lats,
         positive_in=True,
-        mask=None,
+        mask=None
         ):
         """
         TO DO
@@ -196,7 +209,7 @@ def initiate_from_boundary(
             get_region_boundary_grid_indices(
                 lons.copy(),
                 lats.copy(),
-                self.grid,
+                self.grid
             )
         )
         if mask is None:
@@ -205,4 +218,141 @@ def initiate_from_boundary(
                 self.lats,
                 self.grid
             )
-        self.mask = mask.astype(bool) ^ (not positive_in)
+        self.mask = mask.astype(bool) ^ (not positive_in)
+
+    def to_gr(self, path):
+        """Save the GriddedRegion object in a .gr format file directory
+
+        There are two key files within each .gr directory:
+          - a `grid.nc` file that contains information about the coordinates
+          requires to create an `xgcm.Grid` instance
+          - a `boundary.nc` file that contains the region's tracer cell mask
+          and the coordinates and indices of the corner cells that define its
+          boundary.
+
+        To do:
+          - Subdirectory for child boundary information
+
+        Arguments
+        ---------
+        path [str] -- path to directory where the .gr file directory
+            should be saved. The filename will be [GriddedRegion.name].gr
+
+        Example
+        -------
+        >>> gridded_region.save("../data/")
+        """
+        gr_path = f"{path}/{self.name.replace(" ","_")}.gr/"
+        Path(gr_path).mkdir(parents=True, exist_ok=True)
+
+        grid_path = f"{gr_path}/grid.nc"
+        parent_grid = f"{path}/../grid.nc"
+        if os.path.isfile(parent_grid):
+            os.symlink(parent_grid, grid_path)
+        else:
+            grid = self.grid
+            grid._ds.drop_vars([v for v in grid._ds.data_vars]).to_netcdf(grid_path)
+
+        # Write boundary information to NetCDF file
+        vertex = xr.DataArray(np.arange(0, self.i.size), dims=("vertex",))
+        face = xr.DataArray(np.arange(0.5, self.i.size-1), dims=("face",))
+        ds = xr.Dataset({}, coords={"vertex": vertex, "face": face})
+        for v in ["lons", "lats", "i", "j"]:
+            var = getattr(self,v)
+            if v in ["lons", "lats"]:
+                var = loop(var)
+            ds[v] = xr.DataArray(var, dims=("vertex",))
+        for v in ["lons_uv", "lats_uv"]:
+            ds[v] = xr.DataArray(getattr(self,v), dims=("face",))
+        ds["mask"] = self.mask
+        ds.to_netcdf(f"{gr_path}/region.nc")
+
+        for (k,v) in self.save.items():
+            v.to_netcdf(f"{gr_path}/{k}.nc")
+
+        # Write boundary information for each child section to NetCDF file
+        child_path = f"{gr_path}/children/"
+        Path(child_path).mkdir(parents=True, exist_ok=True)
+        for child in self.children.values():
+            child_name = child.name.replace(" ","_")
+            sec_path = f"{gr_path}/children/{child_name}.sec"
+            Path(sec_path).mkdir(parents=True, exist_ok=True)
+
+            vertex = xr.DataArray(np.arange(0, child.i.size), dims=("vertex",))
+            ds = xr.Dataset({}, coords={"vertex": vertex})
+            for k in ["lons", "lats", "i", "j"]:
+                var = getattr(child,k)
+                ds[k] = xr.DataArray(var, dims=("vertex",))
+            ds.to_netcdf(f"{sec_path}/section.nc")
+            for k, ds_save in child.save.items():
+                ds_save.to_netcdf(f"{sec_path}/{k}.nc")
+
+class BoundedRegion(GriddedRegion):
+    def __init__(self, section, grid, **kwargs):
+        super().__init__(
+            section.name,
+            section.lons,
+            section.lats,
+            grid,
+            **kwargs
+        )
+        self.children = {}
+        for child_name, child in section.children.items():
+            i, j, lons, lats = sec.grid_section(grid, child.lons, child.lats)
+            coords = sec.coords_from_lonlat(lons, lats)
+            child_section = sec.Section(
+                child_name,
+                coords,
+                children={},
+                parents={}
+            )
+            child_section.i = i
+            child_section.j = j
+            self.children[child_name] = child_section
+
+def open_gr(path, ds_to_grid):
+
+    ds_grid = xr.open_dataset(f"{path}/grid.nc")
+    grid = ds_to_grid(ds_grid)
+    ds = xr.open_dataset(f"{path}/region.nc")
+
+    name = path.split("/")[-1][:-3].replace("_", " ")
+    region = GriddedRegion(
+        name,
+        ds.lons.values,
+        ds.lats.values,
+        grid,
+        mask = ds.mask,
+        ij = (ds.i.values, ds.j.values)
+    )
+    gr_files = [
+        f for f in os.listdir(f"{path}/")
+        if (".nc" in f) and (f not in ["grid.nc", "region.nc"])
+    ]
+    for file in gr_files:
+        v = file.split(".")[0]
+        region.save[v] = xr.open_dataset(f"{path}/{file}")
+
+    region.children = {}
+    children_path = f"{path}/children/"
+    child_paths = [
+        f"{children_path}{file}"
+        for file in os.listdir(children_path)
+    ]
+    for child_path in child_paths:
+        child_name = child_path.split("/")[-1][:-4].replace("_", " ")
+        ds = xr.open_dataset(f"{child_path}/section.nc")
+
+        section = sec.Section(
+            child_name,
+            sec.coords_from_lonlat(ds.lons.values, ds.lats.values)
+        )
+
+        section.save = {}
+        for file in [f for f in os.listdir(child_path) if f != "section.nc"]:
+            v = file.split(".")[0]
+            section.save[v] = xr.open_dataset(f"{child_path}/{file}")
+
+        region.children[child_name] = section
+
+    return region

regionate/regions.py

@@ -1,10 +1,14 @@
 import numpy as np
+import xarray as xr
 
-from .region import Region, GriddedRegion
+from .region import Region, GriddedRegion, BoundedRegion, open_gr
 from .boundaries import grid_boundaries_from_mask
 from .overlaps import *
 from .utilities import *
 
+import os
+from pathlib import Path
+
 class Regions():
     """
     A dictionary of polygonal regions defined by a list or array of geographical coordinates.
@@ -93,7 +97,7 @@ def __init__(
 
         PARAMETERS
         ----------
-        region_dict : dictionary mapping region `name` (str) to `Region` or `GriddedRegion` instance
+        region_dict : dictionary mapping region `name` (str) to `Region`, `GriddedRegion`, or `BoundedRegion` instance
         name : str or None (default: None)
             Overarching name of the collection of regions
 
@@ -103,14 +107,31 @@ def __init__(
         """
         self.grid = grid
 
-        try:
-            if all([type(v) in [Region, GriddedRegion] for v in region_dict.values()]):
-                super().__init__(region_dict, name=name)
-            else:
-                raise NameError("Values in `region_dict` dictionary must be of instances of `Region` or `GriddedRegion`.")
-        except:
-            raise NameError("Must provide valid `region_dict` dictionary to initialize.")
-
+        super().__init__(region_dict, name=name)
+        # try:
+        #     if all([type(v) in [Region, GriddedRegion, BoundedRegion] for v in region_dict.values()]):
+        #         super().__init__(region_dict, name=name)
+        #     else:
+        #         raise NameError("""Values in `region_dict` dictionary must be instances of
+        #         `Region`, `GriddedRegion`, or `BoundedRegion`.""")
+        # except:
+        #     raise NameError("Must provide valid `region_dict` dictionary to initialize.")
+
+    def to_grs(self, path):
+
+        # Create .grs file directory
+        grs_path = f"{path}{self.name.replace(" ","_")}.grs"
+        Path(grs_path).mkdir(parents=True, exist_ok=True)
+
+        # Write grid dataset (without variables) to NetCDF file
+        grid = self.grid
+        grid_path = f"{grs_path}/grid.nc"
+        grid._ds.drop_vars([v for v in grid._ds.data_vars]).to_netcdf(grid_path)
+
+        # Write boundary information for each region to separate .gr file directory
+        for region_name, region in self.region_dict.items():
+            region.to_gr(f"{grs_path}/")
+
 class MaskRegions(GriddedRegions):
     """
     A dictionary of polygonal regions that exactly conform to the velocity faces bounding a mask in a C-grid ocean model.
@@ -159,6 +180,17 @@ def __init__(
         }
         super().__init__(region_dict, grid, name=name)
 
+def open_grs(path, ds_to_grid):
+    name = path.split("/")[-1][:-4]
+    grid = ds_to_grid(xr.open_dataset(f"{path}/grid.nc"))
+    grs_files = [f for f in os.listdir(f"{path}/") if f!="grid.nc"]
+    region_dict = {}
+    for gr_file in grs_files:
+        gr_name = gr_file[:-3]
+        region_dict[gr_name] = open_gr(f"{path}/{gr_file}", ds_to_grid)
+
+    return GriddedRegions(region_dict, grid, name=name)
+
 def sorted_tuple(s):
     """Sort tuples (by integer)"""
     return tuple(sorted(s, key=int))

regionate/version.py

@@ -1,3 +1,3 @@
 """regionate: version information"""
 
-__version__ = "0.1.1"
+__version__ = "0.2.0"