Merge pull request #233 from amoodie/output_coords

Change default configuration of netcdf file output coordinates

docs/source/examples/custom_saving.rst

@@ -17,11 +17,11 @@ For example, ``self._save_fig_list['active_layer'] = ['active_layer']`` will pro
 When adding variables or metadata to be initialized and subsequently saved in the output netCDF, the key-value pair relationship is as follows.
 The key added to ``self._save_var_list`` is the name of the variable as it will be recorded in the netCDF file, this *does not* have to correspond to the name of an attribute in the model.
 To add a variable to the metadata, a key must be added to ``self._save_var_list['meta']``.
-The expected value for a given key is a list containing strings indicating the model attribute to be saved, its units, the variable type, and lastly the variable dimensions (e.g., ``['active_layer', 'fraction', 'f4', ('total_time', 'length', 'width')]`` for the active layer).
+The expected value for a given key is a list containing strings indicating the model attribute to be saved, its units, the variable type, and lastly the variable dimensions (e.g., ``['active_layer', 'fraction', 'f4', ('time', 'x', 'y')]`` for the active layer).
 
 .. important::
 
-    The dimensions of the custom variable being specified must match *exactly* with one of the three standard dimensions: length, width, total_time.
+    The dimensions of the custom variable being specified must match *exactly* with one of the three standard dimensions: `x`, `y`, `time`.
     Use of an invalid dimension will result in an error.
 
 An example of using the hook and creating a model subclass to customize the figures, gridded variables, and metadata being saved is provided below.
@@ -49,8 +49,8 @@ An example of using the hook and creating a model subclass to customize the figu
     ...
     ...         # save the active layer grid each save_dt w/ a short name
     ...         self._save_var_list['actlay'] = ['active_layer', 'fraction',
-    ...                                          'f4', ('total_time',
-    ...                                                 'length', 'width')]
+    ...                                          'f4', ('time',
+    ...                                                 'x', 'y')]
     ...
     ...         # save number of water parcels w/ a long name
     ...         self._save_var_list['meta']['water_parcels'] = ['Np_water',
@@ -83,4 +83,4 @@ For simplicity we will just check that the appropriate parameters were added to
     {'active_layer': ['active_layer']}
 
     >>> print(mdl._save_var_list)
-    {'meta': {'water_parcels': ['Np_water', 'parcels', 'i8', ()]}, 'actlay': ['active_layer', 'fraction', 'f4', ('total_time', 'length', 'width')]}
+    {'meta': {'water_parcels': ['Np_water', 'parcels', 'i8', ()]}, 'actlay': ['active_layer', 'fraction', 'f4', ('time', 'x', 'y')]}

docs/source/examples/variable_velocity.rst

@@ -94,12 +94,12 @@ We define a model subclass to handle the changing boundary condition:
 and then simply run with:
 
 .. plot::
-    :context: close-figs
+    :context:
 
     # we create the model here, just to be sure it works (for good docs)
     with pyDeltaRCM.shared_tools._docs_temp_directory() as output_dir:
         mdl = ChangingVelocityModel(
-            end_time=end_time,
+            end_time=86400*100,
             out_dir=output_dir)
 
 .. code::

docs/source/info/outputfile.rst

@@ -9,23 +9,31 @@ Gridded Variables
 =================
 
 In any given run, the saving parameters "save_<var>_grids" control whether or
-not that 2-D grid variable (e.g. velocity) is saved to the netCDF4 file.
-In the netCDF4 file, a 3-D array with the dimensions
-*time* x *length* x *width* is created for each 2-D grid variable that is set
-to be saved. The appropriate units for these variables are stored as well,
-such as "meters per second" for the *velocity* grid.
+not that 2-D grid variable (e.g. velocity) is saved to the netCDF4 file. In
+the netCDF4 file, a 3-D array with the dimensions `time` :math:`\times`
+`x` :math:`\times` `y` is created for each 2-D grid variable that is set to
+be saved. Note that `x` is the *downstream* coordinate, rather than the
+Cartesian `x` when displaying the grid. The appropriate units for all
+variables are stored: for example "meters per second" for the *velocity*
+grid.
+
+.. note::
+
+   The format of the output netCDF file coordinate changed in `v2.1.0`. The
+   old format is documented
+   in :attr:`~pyDeltaRCM.model.DeltaModel.legacy_netcdf`, and that input
+   parameter `legacy_netcdf` can be used to create on output netcdf file with
+   the old coordinate configuration.
 
 
 Grid Coordinates
 ================
 
-To save the model information associated with the domain itself, variables
-associated with the grid are saved as well. These are the meshed 2-D grids
-corresponding to the distance of each cell from the boundary in the "Width"
-dimension of the domain, *x* in meters. As well as the distance away from the
-boundary of each cell in the "Length" dimension, as *y* in meters. Similarly, a
-*time* variable is stored which is a 1-D array (vector) holding the model time
-values in seconds, associated with each set of saved output data.
+Grid coordinates are specified in the variables `time`, `x`, and `y` in the output netCDF4 file.
+These arrays are 1D arrays, which specify the location of each cell in the domain in *dimensional* coordinates (e.g., meters).
+In the downstream direction,  the distance of each cell from the inlet boundary is specified in `x` in meters.
+Similarly, the cross-domain distance is specified in `y` in meters.
+Lastly, the `time` variable is stored as a 1D array with model `time` in seconds.
 
 
 Model Metadata

docs/source/info/yamlparameters.rst

@@ -114,6 +114,8 @@ Output Settings
 
 :attr:`pyDeltaRCM.model.DeltaModel.clobber_netcdf`
 
+:attr:`pyDeltaRCM.model.DeltaModel.legacy_netcdf`
+
 
 Reduced-Complexity Routing Parameters
 =====================================

pyDeltaRCM/_version.py

@@ -4,4 +4,4 @@ def __version__():
     Private version declaration, gets assigned to pyDeltaRCM.__version__
     during import
     """
-    return '2.0.3'
+    return '2.1.0'

pyDeltaRCM/default.yml

@@ -181,3 +181,6 @@ sand_frac_bc:
 clobber_netcdf:
   type: 'bool'
   default: False
+legacy_netcdf:
+  type: 'bool'
+  default: False

pyDeltaRCM/init_tools.py

@@ -176,6 +176,10 @@ def import_files(self, kwargs_dict={}):
         self._input_file_vars = input_file_vars
         self.out_dir = self._input_file_vars['out_dir']
         self.verbose = self._input_file_vars['verbose']
+        if self._input_file_vars['legacy_netcdf']:
+            self._netcdf_coords = ('total_time', 'length', 'width')
+        else:
+            self._netcdf_coords = ('time', 'x', 'y')
 
     def process_input_to_model(self):
         """Process input file to model variables.
@@ -275,6 +279,7 @@ def create_other_variables(self):
                                 self._save_velocity_components)
         if self._save_any_grids:  # always save metadata if saving grids
             self._save_metadata = True
+
         self._is_finalized = False
 
     def set_constants(self):
@@ -430,12 +435,15 @@ def create_domain(self):
         # resolve any boundary conditions
         self._hb = self.hb or self.h0  # basin depth
 
-        # ---- empty arrays ----
+        # ---- coordinates ----
+        self.xc = np.arange(0, self.L) * self._dx
+        self.yc = np.arange(0, self.W) * self._dx
         self.x, self.y = np.meshgrid(np.arange(0, self.W),
                                      np.arange(0, self.L))
         self.X, self.Y = np.meshgrid(np.arange(0, self.W+1)*self._dx,
                                      np.arange(0, self.L+1)*self._dx)
 
+        # ---- empty arrays ----
         self.cell_type = np.zeros((self.L, self.W), dtype=np.int64)
         self.eta = np.zeros((self.L, self.W), dtype=np.float32)
         self.eta0 = np.copy(self.eta)  # establish eta0 copy
@@ -593,25 +601,39 @@ def init_output_file(self):
             self.output_netcdf.description = 'Output from pyDeltaRCM'
             self.output_netcdf.history = ('Created '
                                           + time_lib.ctime(time_lib.time()))
-            self.output_netcdf.source = 'pyDeltaRCM'
+            self.output_netcdf.source = 'pyDeltaRCM v{ver}'.format(
+                ver=self.__pyDeltaRCM_version__)
 
-            # create master dimensions
-            self.output_netcdf.createDimension('length', self.L)
-            self.output_netcdf.createDimension('width', self.W)
-            self.output_netcdf.createDimension('total_time', None)
+            # create master dimensions (pulls from `self._netcdf_coords`)
+            self.output_netcdf.createDimension(self._netcdf_coords[1], self.L)
+            self.output_netcdf.createDimension(self._netcdf_coords[2], self.W)
+            self.output_netcdf.createDimension(self._netcdf_coords[0], None)
 
             # create master coordinates (as netCDF variables)
-            x = self.output_netcdf.createVariable(
-                'x', 'f4', ('length', 'width'))
-            y = self.output_netcdf.createVariable(
-                'y', 'f4', ('length', 'width'))
-            time = self.output_netcdf.createVariable('time', 'f4',
-                                                     ('total_time',))
-            x.units = 'meters'
-            y.units = 'meters'
+            time = self.output_netcdf.createVariable(
+                'time', 'f4', (self._netcdf_coords[0],))
             time.units = 'second'
-            x[:] = self.x
-            y[:] = self.y
+
+            if self._legacy_netcdf:
+                # old format is 2d array x and y
+                x = self.output_netcdf.createVariable(
+                    'x', 'f4', self._netcdf_coords[1:])
+                y = self.output_netcdf.createVariable(
+                    'y', 'f4', self._netcdf_coords[1:])
+                x[:] = self.x
+                y[:] = self.y
+
+            else:
+                # new output format is 1d x and y
+                x = self.output_netcdf.createVariable(
+                    'x', 'f4', ('x'))
+                y = self.output_netcdf.createVariable(
+                    'y', 'f4', ('y'))
+                x[:] = self.xc
+                y[:] = self.yc
+
+            x.units = 'meter'
+            y.units = 'meter'
 
             # set up variables for output data grids
             def _create_grid_variable(varname, varunits,
@@ -692,25 +714,28 @@ def init_metadata_list(self):
         self._save_var_list['meta']['hb'] = ['hb', 'meters', 'f4', ()]
         self._save_var_list['meta']['cell_type'] = ['cell_type',
                                                     'type', 'i8',
-                                                    ('length', 'width')]
+                                                    self._netcdf_coords[1:]]
         # subsidence metadata
         if self._toggle_subsidence:
             self._save_var_list['meta']['start_subsidence'] = [
                 'start_subsidence', 'seconds', 'i8', ()
             ]
             self._save_var_list['meta']['sigma'] = [
                 'sigma', 'meters per timestep', 'f4',
-                ('length', 'width')
+                self._netcdf_coords[1:]
             ]
         # time-varying metadata
         self._save_var_list['meta']['H_SL'] = [None, 'meters', 'f4',
-                                               ('total_time')]
+                                               (self._netcdf_coords[0])]
         self._save_var_list['meta']['f_bedload'] = [None, 'fraction',
-                                                    'f4', ('total_time')]
+                                                    'f4',
+                                                    (self._netcdf_coords[0])]
         self._save_var_list['meta']['C0_percent'] = [None, 'percent',
-                                                     'f4', ('total_time')]
+                                                     'f4',
+                                                     (self._netcdf_coords[0])]
         self._save_var_list['meta']['u0'] = [None, 'meters per second',
-                                             'f4', ('total_time')]
+                                             'f4',
+                                             (self._netcdf_coords[0])]
 
     def load_checkpoint(self, defer_output=False):
         """Load the checkpoint from the .npz file.

pyDeltaRCM/model.py

@@ -809,7 +809,7 @@ def save_eta_grids(self, save_eta_grids):
         if (save_eta_grids is True) and \
           ('eta' not in self._save_var_list.keys()):
             self._save_var_list['eta'] = ['eta', 'meters', 'f4',
-                                          ('total_time', 'length', 'width')]
+                                          self._netcdf_coords]
         elif ((save_eta_grids is False) and
               ('eta' in self._save_var_list.keys())):
             del self._save_var_list['eta']
@@ -827,7 +827,7 @@ def save_stage_grids(self, save_stage_grids):
         if (save_stage_grids is True) and \
           ('stage' not in self._save_var_list.keys()):
             self._save_var_list['stage'] = ['stage', 'meters', 'f4',
-                                            ('total_time', 'length', 'width')]
+                                            self._netcdf_coords]
         elif ((save_stage_grids is False) and
               ('stage' in self._save_var_list.keys())):
             del self._save_var_list['stage']
@@ -845,7 +845,7 @@ def save_depth_grids(self, save_depth_grids):
         if (save_depth_grids is True) and \
           ('depth' not in self._save_var_list.keys()):
             self._save_var_list['depth'] = ['depth', 'meters', 'f4',
-                                            ('total_time', 'length', 'width')]
+                                            self._netcdf_coords]
         elif ((save_depth_grids is False) and
               ('depth' in self._save_var_list.keys())):
             del self._save_var_list['depth']
@@ -865,8 +865,7 @@ def save_discharge_grids(self, save_discharge_grids):
             self._save_var_list['discharge'] = ['qw',
                                                 'cubic meters per second',
                                                 'f4',
-                                                ('total_time', 'length',
-                                                 'width')]
+                                                self._netcdf_coords]
         elif ((save_discharge_grids is False) and
               ('discharge' in self._save_var_list.keys())):
             del self._save_var_list['discharge']
@@ -884,8 +883,7 @@ def save_velocity_grids(self, save_velocity_grids):
         if (save_velocity_grids is True) and \
           ('velocity' not in self._save_var_list.keys()):
             self._save_var_list['velocity'] = ['uw', 'meters per second', 'f4',
-                                               ('total_time', 'length',
-                                                'width')]
+                                               self._netcdf_coords]
         elif ((save_velocity_grids is False) and
               ('velocity' in self._save_var_list.keys())):
             del self._save_var_list['velocity']
@@ -904,8 +902,7 @@ def save_sedflux_grids(self, save_sedflux_grids):
           ('sedflux' not in self._save_var_list.keys()):
             self._save_var_list['sedflux'] = ['qs', 'cubic meters per second',
                                               'f4',
-                                              ('total_time', 'length',
-                                               'width')]
+                                              self._netcdf_coords]
         elif ((save_sedflux_grids is False) and
               ('sedflux' in self._save_var_list.keys())):
             del self._save_var_list['sedflux']
@@ -924,8 +921,7 @@ def save_sandfrac_grids(self, save_sandfrac_grids):
         if (save_sandfrac_grids is True) and \
           ('sandfrac' not in self._save_var_list.keys()):
             self._save_var_list['sandfrac'] = ['sand_frac', 'fraction', 'f4',
-                                               ('total_time', 'length',
-                                                'width')]
+                                               self._netcdf_coords]
         elif ((save_sandfrac_grids is False) and
               ('sandfrac' in self._save_var_list.keys())):
             del self._save_var_list['sandfrac']
@@ -944,11 +940,11 @@ def save_discharge_components(self, save_discharge_components):
             if ('discharge_x' not in self._save_var_list.keys()):
                 self._save_var_list['discharge_x'] = [
                     'qx', 'cubic meters per second', 'f4',
-                    ('total_time', 'length', 'width')]
+                    self._netcdf_coords]
             if ('discharge_y' not in self._save_var_list.keys()):
                 self._save_var_list['discharge_y'] = [
                     'qy', 'cubic meters per second', 'f4',
-                    ('total_time', 'length', 'width')]
+                    self._netcdf_coords]
         elif (save_discharge_components is False):
             if ('discharge_x' in self._save_var_list.keys()):
                 del self._save_var_list['discharge_x']
@@ -969,11 +965,11 @@ def save_velocity_components(self, save_velocity_components):
             if ('velocity_x' not in self._save_var_list.keys()):
                 self._save_var_list['velocity_x'] = [
                     'ux', 'meters per second', 'f4',
-                    ('total_time', 'length', 'width')]
+                    self._netcdf_coords]
             if ('velocity_y' not in self._save_var_list.keys()):
                 self._save_var_list['velocity_y'] = [
                     'uy', 'meters per second', 'f4',
-                    ('total_time', 'length', 'width')]
+                    self._netcdf_coords]
         elif (save_velocity_components is False):
             if ('velocity_x' in self._save_var_list.keys()):
                 del self._save_var_list['velocity_x']
@@ -1263,6 +1259,47 @@ def clobber_netcdf(self):
     def clobber_netcdf(self, clobber_netcdf):
         self._clobber_netcdf = clobber_netcdf
 
+    @property
+    def legacy_netcdf(self):
+        """Enable output in legacy netCDF format.
+
+        .. note:: new in `v2.1.0`.
+
+        Default behavior, legacy_netcdf: False, is for the model to use the
+        new `v2.1.0` output netCDF format. The updated format is configured
+        to match the input expected by `xarray`, which eases interaction with
+        model outputs. The change in format is from inconsistently named
+        dimensions and *coordinate variables*, to homogeneous definitions.
+        Also, the legacy format specified the variables `x` and `y` as 2d
+        grids, whereas the updated format uses 1d coordinate arrays.
+        
+        .. important::
+
+            There are no changes to the dimensionality of data such as bed
+            elevation or velocity, only the metadata specifying the location of
+            the data are changed.
+
+        +-------------+-------------------+---------------------------------+
+        |             | default           | legacy                          |
+        +=============+===================+=================================+
+        | dimensions  | `time`, `x`, `y`  | `total_time`, `length`, `width` |
+        +-------------+-------------------+---------------------------------+
+        | variables   | `time`, `x`, `y`  | `time`, `y`, `x`; x, y as 2D    |
+        +-------------+-------------------+---------------------------------+
+        | data        | `t-x-y` array     | `t-y-x` array                   |
+        +-------------+-------------------+---------------------------------+
+
+        .. hint::
+
+            If you are beginning a new project, use `legacy_netcdf == False`,
+            and update scripts accordingly.
+        """
+        return self._legacy_netcdf
+
+    @legacy_netcdf.setter
+    def legacy_netcdf(self, legacy_netcdf):
+        self._legacy_netcdf = legacy_netcdf
+
     @property
     def time(self):
         """Elapsed model time in seconds.