Merge pull request #50 from IMMM-SFA/joss-review

Update documentation; add Jupyter notebook tutorial

.github/workflows/build_linux.yml

@@ -12,7 +12,7 @@ jobs:
 
       env:
         OS: ${{ matrix.os }}
-        PYTHON: '3.9'
+        PYTHON: '3.7'
 
       steps:
 
@@ -21,15 +21,14 @@ jobs:
         - name: Set up Python
           uses: actions/setup-python@master
           with:
-            python-version: 3.9
+            python-version: 3.7
 
         - name: Install dependencies
           run: |
             python -m pip install --upgrade pip
-            pip install -r requirements.txt
-            python setup.py install
+            pip install .
 
-        - name: Test and generate coverage report
+        - name: Test
           run: |
             pip install pytest
             pytest

.github/workflows/build_osx.yml

@@ -12,7 +12,7 @@ jobs:
 
       env:
         OS: ${{ matrix.os }}
-        PYTHON: '3.9'
+        PYTHON: '3.7'
 
       steps:
 
@@ -21,13 +21,12 @@ jobs:
         - name: Set up Python
           uses: actions/setup-python@master
           with:
-            python-version: 3.9
+            python-version: 3.7
 
         - name: Install dependencies
           run: |
             python -m pip install --upgrade pip
-            pip install -r requirements.txt
-            python setup.py install
+            pip install .
 
         - name: Test and generate coverage report
           run: |

.github/workflows/build_windows.yml

@@ -1,7 +1,6 @@
 name: windows
 
-#on: [push]
-on: workflow_dispatch
+on: [push]
 
 jobs:
   build:
@@ -13,7 +12,7 @@ jobs:
 
       env:
         OS: ${{ matrix.os }}
-        PYTHON: '3.9'
+        PYTHON: '3.7'
 
       steps:
 
@@ -22,23 +21,20 @@ jobs:
         - name: Set up Python
           uses: actions/setup-python@master
           with:
-            python-version: 3.9
+            python-version: 3.7
 
         - name: Install dependencies
           run: |
             python -m pip install --upgrade pip
-            curl -O https://download.lfd.uci.edu/pythonlibs/w4tscw6k/numpy-1.20.3+mkl-cp39-cp39-win_amd64.whl
-            pip install numpy-1.20.3+mkl-cp39-cp39-win_amd64.whl
-            curl -O https://download.lfd.uci.edu/pythonlibs/w4tscw6k/GDAL-3.3.0-cp39-cp39-win_amd64.whl
-            pip install GDAL-3.3.0-cp39-cp39-win_amd64.whl
-            curl -O https://download.lfd.uci.edu/pythonlibs/w4tscw6k/rasterio-1.2.4-cp39-cp39-win_amd64.whl
-            pip install rasterio-1.2.4-cp39-cp39-win_amd64.whl
-            pip install -r requirements.txt
-            python setup.py install
-            pip install numpy --upgrade
-            pip install python-benedict --upgrade
-
-        - name: Test and generate coverage report
+            curl -O https://download.lfd.uci.edu/pythonlibs/q4trcu4l/numpy-1.20.3+mkl-cp37-cp37m-win_amd64.whl
+            pip install numpy-1.20.3+mkl-cp37-cp37m-win_amd64.whl
+            curl -O https://download.lfd.uci.edu/pythonlibs/q4trcu4l/GDAL-3.3.0-cp37-cp37m-win_amd64.whl
+            pip install GDAL-3.3.0-cp37-cp37m-win_amd64.whl
+            curl -O https://download.lfd.uci.edu/pythonlibs/q4trcu4l/rasterio-1.2.4-cp37-cp37m-win_amd64.whl
+            pip install rasterio-1.2.4-cp37-cp37m-win_amd64.whl
+            pip install .
+
+        - name: Test
           run: |
             pip install pytest
             pytest

README.md

@@ -6,11 +6,16 @@
 
 ## getting started
 
-Ensure you have Python >= 3.9 available (consider using a [virtual environment](https://github.com/pyenv/pyenv), see the docs [here](https://mosartwmpy.readthedocs.io/en/latest/virtualenv.html) for a brief tutorial), then install `mosartwmpy` with:
+Ensure you have Python >= 3.7 available (consider using a [virtual environment](https://github.com/pyenv/pyenv), see the docs [here](https://mosartwmpy.readthedocs.io/en/latest/virtualenv.html) for a brief tutorial), then install `mosartwmpy` with:
 ```shell
 pip install mosartwmpy
 ```
 
+Alternatively, install via conda with:
+```shell
+conda install mosartwmpy
+```
+
 Download a sample input dataset spanning May 1981 by running the following and selecting option `1` for "tutorial". This will download and unpack the inputs to your current directory. Optionally specify a path to download and extract to instead of the current directory.
 
 ```shell
@@ -108,6 +113,10 @@ plt.show()
 
 ```
 
+## model coupling
+
+A common use case for `mosartwmpy` is to run coupled with output from the Community Land Model (CLM). To see an example of how to drive `mosartwmpy` with runoff from a coupled model, check out the [Jupyter notebook tutorial](https://github.com/IMMM-SFA/mosartwmpy/blob/main/notebooks/tutorial.ipynb)!
+
 ## model input
 
 Several input files in NetCDF format are required to successfully run a simulation, which are not shipped with this repository due to their large size. The grid files, reservoir files, and a small range of runoff and demand input files can be obtained using the download utility by running `python -m mosartwmpy.download` and choosing option 1 for "sample_input". Currently, all input files are assumed to be at the same resolution (for the sample files this is 1/8 degree over the CONUS). Below is a summary of the various input files:
@@ -197,31 +206,6 @@ Several input files in NetCDF format are required to successfully run a simulati
 </tbody>
 </table>
 
-Alternatively, certain model inputs can be set using the BMI interface. This can be useful for coupling `mosartwmpy` with other models. If setting an input that would typically be read from a file, be sure to disable the `read_from_file` configuration value for that input. For example:
-```python
-import numpy as np
-from mosartwmpy import Model
-
-mosart_wm = Model()
-mosart_wm.initialize()
-
-# get a list of model input variables
-mosart_wm.get_input_var_names()
-
-# disable the runoff read_from_file
-mosart_wm.config['runoff.read_from_file'] = False
-
-# set the runoff values manually (i.e. from another model's output)
-surface_runoff = np.empty(mosart_wm.get_grid_size())
-surface_runoff[:] = # <values from coupled model>
-mosart_wm.set_value('surface_runoff_flux', surface_runoff)
-
-# advance one timestep
-mosart_wm.update()
-
-# continue coupling...
-```
-
 ## model output
 
 By default, key model variables are output on a monthly basis at a daily averaged resolution to `./output/<simulation name>/<simulation name>_<year>_<month>.nc`. See the configuration file for examples of how to modify the outputs, and the `./mosartwmpy/state/state.py` file for state variable names.

mosartwmpy/config_defaults.yaml

@@ -229,11 +229,11 @@ water_management:
         latitude: lat
         # streamflow field name [m3/s]
         streamflow: Qmon
-        # streamflow time resolution - 'month' or 'epiweek'
+        # streamflow time resolution - only 'month' supported so far
         streamflow_time_resolution: month
         # demand field name [m3/s]
         demand: demand
-        # demand time resolution - 'month' or 'epiweek'
+        # demand time resolution - only 'month' supported so far
         demand_time_resolution: month
         # grid id to reservoir id mapping field
         grid_to_reservoir: gridID_from_Dam

mosartwmpy/grid/grid.py

@@ -404,11 +404,11 @@ def to_files(self, path: str) -> None:
             for df in dfs:
                 names.append(f'{df["key"]}.df.nc')
                 paths.append(f'{tmpdir}/{names[-1]}')
-                df['frame'].to_xarray().to_netcdf(paths[-1], engine='h5netcdf')
+                df['frame'].to_xarray().to_netcdf(paths[-1], engine='scipy')
             for ds in xrs:
                 names.append(f'{ds["key"]}.xr.nc')
                 paths.append(f'{tmpdir}/{names[-1]}')
-                ds['data_array'].to_netcdf(paths[-1], engine='h5netcdf')
+                ds['data_array'].to_netcdf(paths[-1], engine='scipy')
             with ZipFile(path, 'w', compression=ZIP_DEFLATED, compresslevel=9) as zip:
                 for i, filename in enumerate(paths):
                     zip.write(filename, names[i])
@@ -441,13 +441,13 @@ def from_files(path: Path) -> 'Grid':
                             setattr(grid, key, npdf[key].values)
                     if filename.endswith('df.nc'):
                         key = filename.split('.')[0]
-                        ds = xr.open_dataset(file, engine='h5netcdf')
+                        ds = xr.open_dataset(file, engine='scipy')
                         df = ds.to_dataframe()
                         setattr(grid, key, df)
                         ds.close()
                     if filename.endswith('xr.nc'):
                         key = filename.split('.')[0]
-                        ds = xr.open_dataarray(file, engine='h5netcdf').load()
+                        ds = xr.open_dataarray(file, engine='scipy').load()
                         setattr(grid, key, ds)
                         ds.close()
 

mosartwmpy/model.py

@@ -10,7 +10,6 @@
 from benedict import benedict
 from bmipy import Bmi
 from datetime import datetime, time, timedelta
-from epiweeks import Week
 from pathlib import Path
 from pathvalidate import sanitize_filename
 from timeit import default_timer as timer
@@ -161,22 +160,27 @@ def initialize(self, config_file_path: str = None, grid: Grid = None, state: Sta
             raise e
 
         logging.debug(f'Initialization completed in {pretty_timer(timer() - t)}.')
+        logging.info(f'Done.')
 
     def update(self) -> None:
         t = timer()
         step = datetime.fromtimestamp(self.get_current_time()).isoformat(" ")
         # perform one timestep
         logging.debug(f'Beginning timestep {step}...')
         try:
-            # read runoff
             if self.config.get('runoff.read_from_file', False):
+                # read runoff from file
                 logging.debug(f'Reading runoff input from file.')
                 load_runoff(self.state, self.grid, self.config, self.current_time)
+            else:
+                # convert provided runoff from mm/s to m3/s
+                self.state.hillslope_surface_runoff = 0.001 * self.grid.land_fraction * self.grid.area * self.state.hillslope_surface_runoff
+                self.state.hillslope_subsurface_runoff = 0.001 * self.grid.land_fraction * self.grid.area * self.state.hillslope_subsurface_runoff
+                self.state.hillslope_wetland_runoff = 0.001 * self.grid.land_fraction * self.grid.area * self.state.hillslope_wetland_runoff
             # read demand
             if self.config.get('water_management.enabled', False):
                 if self.config.get('water_management.demand.read_from_file', False):
                     # only read new demand and compute new release if it's the very start of simulation or new time period
-                    # TODO this currently assumes monthly demand input
                     if self.current_time == datetime.combine(self.config.get('simulation.start_date'), time.min) or self.current_time == datetime(self.current_time.year, self.current_time.month, 1):
                         logging.debug(f'Reading demand rate input from file.')
                         # load the demand from file
@@ -187,8 +191,6 @@ def update(self) -> None:
                 self.state.grid_cell_supply[:] = 0
                 self.state.grid_cell_unmet_demand[:] = 0
                 # get streamflow for this time period
-                # TODO this is still written assuming monthly, but here's the epiweek for when that is relevant
-                epiweek = Week.fromdate(self.current_time).week
                 month = self.current_time.month
                 streamflow_time_name = self.config.get('water_management.reservoirs.streamflow_time_resolution')
                 self.state.reservoir_streamflow[:] = self.grid.reservoir_streamflow_schedule.sel({streamflow_time_name: month}).values
@@ -206,10 +208,16 @@ def update(self) -> None:
         except Exception as e:
             logging.exception('Failed to write output or restart file; see below for stacktrace.')
             raise e
-        # clear runoff input arrays
-        self.state.hillslope_surface_runoff[:] = 0
-        self.state.hillslope_subsurface_runoff[:] = 0
-        self.state.hillslope_wetland_runoff[:] = 0
+        if self.config.get('runoff.read_from_file', False):
+            # clear runoff input arrays
+            self.state.hillslope_surface_runoff[:] = 0
+            self.state.hillslope_subsurface_runoff[:] = 0
+            self.state.hillslope_wetland_runoff[:] = 0
+        else:
+            # convert back to mm/s
+            self.state.hillslope_surface_runoff = self.state.hillslope_surface_runoff * 1000.0 / self.grid.land_fraction / self.grid.area
+            self.state.hillslope_subsurface_runoff = self.state.hillslope_subsurface_runoff * 1000.0 / self.grid.land_fraction / self.grid.area
+            self.state.hillslope_wetland_runoff = self.state.hillslope_wetland_runoff * 1000.0 / self.grid.land_fraction / self.grid.area
 
     def update_until(self, time: float) -> None:
         # make sure that requested end time is after now

mosartwmpy/reservoirs/grid.py

@@ -57,7 +57,7 @@ def load_reservoirs(self, config: Benedict, parameters: Parameters) -> None:
     # index by grid cell
     self.reservoir_to_grid_mapping = self.reservoir_to_grid_mapping.set_index('grid_cell_id')
 
-    # prepare the month or epiweek based reservoir schedules mapped to the domain
+    # prepare the month based reservoir schedules mapped to the domain
     prepare_reservoir_schedule(self, config, parameters, reservoirs)
 
     reservoirs.close()
@@ -75,8 +75,6 @@ def prepare_reservoir_schedule(self, config: Benedict, parameters: Parameters, r
     # the reservoir streamflow and demand are specified by the time resolution and reservoir id
     # so let's remap those to the actual mosart domain for ease of use
 
-    # TODO i had wanted to convert these all to epiweeks no matter what format provided, but we don't know what year all the data came from
-
     # streamflow flux
     streamflow_time_name = config.get('water_management.reservoirs.streamflow_time_resolution')
     streamflow = reservoirs[config.get('water_management.reservoirs.streamflow')]
@@ -92,7 +90,7 @@ def prepare_reservoir_schedule(self, config: Benedict, parameters: Parameters, r
         else:
             schedule = concat([schedule, sched], dim=streamflow_time_name)
     self.reservoir_streamflow_schedule = schedule.assign_coords(
-        # if monthly, convert to 1 based month index (instead of starting from 0)
+        # convert to 1 based month index (instead of starting from 0)
         {streamflow_time_name: (streamflow_time_name, schedule[streamflow_time_name].values + (1 if streamflow_time_name == 'month' else 0))}
     ).streamflow
 
@@ -110,12 +108,11 @@ def prepare_reservoir_schedule(self, config: Benedict, parameters: Parameters, r
         else:
             schedule = concat([schedule, sched], dim=demand_time_name)
     self.reservoir_demand_schedule = schedule.assign_coords(
-        # if monthly, convert to 1 based month index (instead of starting from 0)
+        # convert to 1 based month index (instead of starting from 0)
         {demand_time_name: (demand_time_name, schedule[demand_time_name].values + (1 if demand_time_name == 'month' else 0))}
     ).demand
 
     # initialize prerelease based on long term mean flow and demand (Biemans 2011)
-    # TODO this assumes demand and flow use the same timescale :(
     flow_avg = self.reservoir_streamflow_schedule.mean(dim=streamflow_time_name)
     demand_avg = self.reservoir_demand_schedule.mean(dim=demand_time_name)
     prerelease = (1.0 * self.reservoir_streamflow_schedule)

mosartwmpy/reservoirs/reservoirs.py

@@ -1,20 +1,16 @@
 import numpy as np
 
 from datetime import datetime
-from epiweeks import Week
 from benedict.dicts import benedict as Benedict
 
 from mosartwmpy.config.parameters import Parameters
 from mosartwmpy.state.state import State
 from mosartwmpy.grid.grid import Grid
 
-# TODO in fortran mosart there is a StorCalibFlag that affects how storage targets are calculated -- code so far is written assuming that it == 0
 
 def reservoir_release(state, grid, config, parameters, current_time):
     # compute release from reservoirs
-
-    # TODO so much logic was dependent on monthly, so still assuming monthly for now, but here's the epiweek for when that is relevant
-    epiweek = Week.fromdate(current_time).week
+
     month = current_time.month
 
     # if it's the start of the operational year for the reservoir, set it's start of op year storage to the current storage
@@ -30,9 +26,7 @@ def reservoir_release(state, grid, config, parameters, current_time):
 
 def regulation_release(state, grid, config, parameters, current_time):
     # compute the expected monthly release based on Biemans (2011)
-
-    # TODO this is still written assuming monthly, but here's the epiweek for when that is relevant
-    epiweek = Week.fromdate(current_time).week
+
     month = current_time.month
     streamflow_time_name = config.get('water_management.reservoirs.streamflow_time_resolution')
 
@@ -77,8 +71,6 @@ def storage_targets(state: State, grid: Grid, config: Benedict, parameters: Para
 
     # TODO the logic here is really hard to follow... can it be simplified or made more readable?
 
-    # TODO this is still written assuming monthly, but here's the epiweek for when that is relevant
-    epiweek = Week.fromdate(current_time).week
     month = current_time.month
     streamflow_time_name = config.get('water_management.reservoirs.streamflow_time_resolution')
 
@@ -96,7 +88,7 @@ def storage_targets(state: State, grid: Grid, config: Benedict, parameters: Para
     )
     drop = 0 * state.reservoir_month_flood_control_start
     n_month = 0 * drop
-    for m in np.arange(1,13): # TODO assumes monthly
+    for m in np.arange(1,13):
         m_and_condition = (m >= state.reservoir_month_flood_control_start) & (m < state.reservoir_month_flood_control_end)
         m_or_condition = (m >= state.reservoir_month_flood_control_start) | (m < state.reservoir_month_flood_control_end)
         drop = np.where(
@@ -128,25 +120,6 @@ def storage_targets(state: State, grid: Grid, config: Benedict, parameters: Para
         (month >= state.reservoir_month_flood_control_end) |
         (month < state.reservoir_month_start_operations)
     )
-    # TODO this logic exists in fortran mosart but isn't used...
-    # fill = 0 * drop
-    # n_month = 0 * drop
-    # for m in np.arange(1,13): # TODO assumes monthly
-    #     m_condition = (m >= self.state.reservoir_month_flood_control_end.values) &
-    #         (self.reservoir_streamflow_schedule.sel({streamflow_time_name: m}).values > self.reservoir_streamflow_schedule.mean(dim=streamflow_time_name).values) & (
-    #             (first_condition & (m <= self.state.reservoir_month_start_operations)) |
-    #             (second_condition & (m <= 12))
-    #         )
-    #     fill = np.where(
-    #         m_condition,
-    #         fill + np.abs(self.reservoir_streamflow_schedule.mean(dim=streamflow_time_name).values - self.reservoir_streamflow_schedule.sel({streamflow_time_name: m}).values),
-    #         fill
-    #     )
-    #     n_month = np.where(
-    #         m_condition,
-    #         n_month + 1,
-    #         n_month
-    #     )
     state.reservoir_release = np.where(
         (state.reservoir_release> grid.reservoir_streamflow_schedule.mean(dim=streamflow_time_name).values) & (first_condition | second_condition),
         grid.reservoir_streamflow_schedule.mean(dim=streamflow_time_name).values,