Allowing hourly output. Allow height setting pass in

pointsnobal/c_snobal/snobal.pyx

@@ -272,7 +272,8 @@ def do_tstep_grid(input1, input2, output_rec, tstep_rec, mh, params, int first_s
             tstep_info[i].intervals = int(tstep_rec[i]['intervals'])
         if tstep_rec[i]['threshold'] is not None:
             tstep_info[i].threshold = tstep_rec[i]['threshold']
-        tstep_info[i].output = int(tstep_rec[i]['output'])
+        # NOTE: the 'output' flag is intentionally not read; snobal's built-in
+        # output paths are disabled here and outputting is driven from Python.
 
     # start1 = clock()
     cdef OUTPUT_REC_ARR output1_c
@@ -649,7 +650,8 @@ def do_tstep(input1, input2, output_rec, tstep_rec, mh, params, first_step=True)
             tstep_info[i].intervals = int(tstep_rec[i]['intervals'])
         if tstep_rec[i]['threshold'] is not None:
             tstep_info[i].threshold = tstep_rec[i]['threshold']
-        tstep_info[i].output = int(tstep_rec[i]['output'])
+        # NOTE: the 'output' flag is intentionally not read; snobal's built-in
+        # output paths are disabled here and outputting is driven from Python.
 
 
 

pointsnobal/cli.py

@@ -32,23 +32,24 @@ def main():
         "--output_file", type=str, default=None,
         help="Optional path to output file"
     )
+    parser.add_argument(
+        "--output_frequency", type=int, default=24,
+        help="Output frequency in hours"
+    )
     args = parser.parse_args()
 
     LOG.info(f"Reading in {args.filepath}")
     df_inputs = pd.read_csv(
         args.filepath,
         parse_dates=["datetime"], index_col="datetime"
     )
-    # Start and end dates
-    start_date = df_inputs.index.min()
-    end_date = df_inputs.index.max()
 
     output_file = args.output_file or "./pointsnobal_results.csv"
 
     # Run the model for that file
     LOG.info(f"Running pointsnobal...")
     df_out = run_model(
-        start_date, end_date, args.elevation, df_inputs
+        args.elevation, df_inputs, output_frequency_hours=args.output_frequency
     )
     LOG.info(f"Finished pointsnobal, outputting to {output_file}")
     df_out.to_csv(output_file)

pointsnobal/point_model.py

@@ -46,11 +46,17 @@
 
 def initialize_model(
         model_datetimes: pd.DatetimeIndex, elevation: float,
+        z_u: float = 5.0, z_t: float = 2.0, z_g: float = 0.3,
 ):
     """
     Args:
         model_datetimes: datetime index from the forcing data
         elevation: elevation in meters
+        z_u: wind speed measurement height in meters (relative to the snow
+            surface)
+        z_t: air temperature measurement height in meters (relative to the
+            snow surface)
+        z_g: soil temperature depth in meters
 
     Returns:
         output_record: output dictionary for start (mostly 0.0s)
@@ -61,37 +67,49 @@ def initialize_model(
     # initialize isnobal state
     LOG.info('Initializing snobal Model')
     freq = pd.infer_freq(model_datetimes)
-    # Convert to offset and calculate number of minutes
+    if freq is None:
+        raise ValueError(
+            "Could not infer a regular frequency from the input datetime "
+            "index. pointsnobal requires evenly-spaced input data."
+        )
+    # Seconds between input records (the data timestep)
     offset = pd.tseries.frequencies.to_offset(freq)
+    data_tstep = pd.Timedelta(offset).total_seconds()
+    # The normal/medium/small sub-timestep hierarchy below assumes the data
+    # step is a whole number of hours: the normal step is one hour and is run
+    # once per hour of the data step. Reject anything that would not divide
+    # cleanly rather than silently under-integrating the timestep.
+    if data_tstep % 3600 != 0:
+        raise ValueError(
+            "pointsnobal requires input data on a whole-hour timestep, got "
+            f"{data_tstep / 3600.0:.4g} hours between records."
+        )
+    normal_intervals = int(data_tstep // 3600)
+
+    # Only the keys consumed by the C layer (via the PARAMS struct) are kept
     constants = {
-            'time_step': offset.n * 60,
             'max_h2o_vol': 0.01,
-            'c': True,
-            'K': True,
-            'mass_threshold': 60,
-            'time_z': 0,
             'max_z_s_0': 0.25,
-            'z_u': 5.0,
-            'z_t': 2.0,
-            'z_g': 0.3,
+            'z_u': z_u,
+            'z_t': z_t,
+            'z_g': z_g,
             'relative_heights': True,
-            'max_density': 550,
-            'max_compact_density': 500,
-            'max_liquid_density': 500,
         }
 
-    # get the timestep info
+    # get the timestep info. Each data step is run as `normal_intervals`
+    # one-hour normal steps, which are adaptively subdivided into medium
+    # (15 min) and small (1 min) steps when the pack mass is below threshold.
     tstep_info = [
         {
-            'level': 0, 'output': 2, 'threshold': None,
-            'time_step': offset.n * 3600.0,  # Datatstep in seconds
+            'level': 0, 'threshold': None,
+            'time_step': data_tstep,  # data timestep in seconds
             'intervals': None
         },
-        {'level': 1, 'output': False, 'threshold': 60.0, 'time_step': 3600.0,
-        'intervals': 1},
-        {'level': 2, 'output': False, 'threshold': 10.0, 'time_step': 900.0,
+        {'level': 1, 'threshold': 60.0, 'time_step': 3600.0,
+        'intervals': normal_intervals},
+        {'level': 2, 'threshold': 10.0, 'time_step': 900.0,
         'intervals': 4},
-        {'level': 3, 'output': False, 'threshold': 1.0, 'time_step': 60.0,
+        {'level': 3, 'threshold': 1.0, 'time_step': 60.0,
         'intervals': 15}
     ]
     # get init params
@@ -185,27 +203,34 @@ def save_timsteps(
 
 
 def run_model(
-        start: pd.Timestamp, end: pd.Timestamp, elevation: float,
-        df_inputs: pd.DataFrame
+        elevation: float,
+        df_inputs: pd.DataFrame,
+        z_u: float = 5.0, z_t: float = 2.0, z_g: float = 0.3,
+        output_frequency_hours: int = 24,
 ) -> pd.DataFrame:
     """
     Run snobal with given input data
     Args:
-        start: start date
-        end: end date
         elevation: elevation in meters for the point
         df_inputs: hourly input pd.Dataframe
+        z_u: wind speed measurement height in meters (relative to the snow
+            surface)
+        z_t: air temperature measurement height in meters (relative to the
+            snow surface)
+        z_g: soil temperature depth in meters
+        output_frequency_hours: how often (in hours) to record an output row.
+            Defaults to 24 (daily). Use 1 for hourly output. Accumulated
+            outputs (e.g. SWI) are summed over each output interval, so at
+            hourly output they represent hourly totals rather than daily.
 
     Returns:
-        Dataframe of daily outputs indexed on datetime
+        Dataframe of outputs indexed on datetime
     """
-    # Returns '=1H'
-    frequency = pd.infer_freq(df_inputs.index)
     # Variable for storing the outputs
     output_list = []
     # Get the variables for snobal
     output_record, tstep_info, constants, model_datetimes = initialize_model(
-        df_inputs.index, elevation)
+        df_inputs.index, elevation, z_u=z_u, z_t=z_t, z_g=z_g)
 
     # Tracking how often we output
     output_record['current_time'] = 1.0 * np.zeros(
@@ -244,8 +269,8 @@ def run_model(
         # copy the second inputs to now be the starting inputs
         input1 = copy.deepcopy(input2)
 
-        # output at the frequency and the last time step
-        if (j * (data_tstep / 3600.0) % 24 == 0) \
+        # output at the requested frequency and the last time step
+        if (j * (data_tstep / 3600.0) % output_frequency_hours == 0) \
                 or (j == len(model_datetimes) - 1):
             LOG.debug('Outputting {}'.format(tstep))
 

tests/test_run_snobal.py

@@ -2,13 +2,82 @@
 import pytest
 from pathlib import Path
 
-from pointsnobal.point_model import run_model
+from pointsnobal.point_model import initialize_model, run_model
+
+
+class TestInitializeModel:
+    """
+    Unit tests for how the snobal timestep hierarchy and parameters are
+    constructed. These do not run the model, so they are fast and do not
+    depend on the compiled extension producing specific values.
+    """
+
+    @staticmethod
+    def _index(freq, periods=48):
+        return pd.date_range("2022-10-01", periods=periods, freq=freq)
+
+    @pytest.mark.parametrize("freq, expected_intervals", [
+        ("1H", 1),
+        ("3H", 3),
+        ("6H", 6),
+    ])
+    def test_normal_intervals_track_frequency(self, freq, expected_intervals):
+        # The NORMAL (1 hr) level must run once per hour of the data step so
+        # that multi-hour input integrates the whole step, not just one hour.
+        _, tstep_info, _, _ = initialize_model(self._index(freq), 2103.0)
+        assert tstep_info[1]["intervals"] == expected_intervals
+
+    @pytest.mark.parametrize("freq, expected_seconds", [
+        ("1H", 3600.0),
+        ("6H", 21600.0),
+    ])
+    def test_data_timestep_in_seconds(self, freq, expected_seconds):
+        _, tstep_info, _, _ = initialize_model(self._index(freq), 2103.0)
+        assert tstep_info[0]["time_step"] == expected_seconds
+
+    def test_constants_only_hold_c_params_keys(self):
+        # constants is passed to the C layer as both `mh` and `params`; only
+        # these six keys are read from it (via the PARAMS struct).
+        _, _, constants, _ = initialize_model(self._index("1H"), 2103.0)
+        assert set(constants) == {
+            "z_u", "z_t", "z_g",
+            "relative_heights", "max_h2o_vol", "max_z_s_0",
+        }
+
+    def test_measurement_heights_passed_through(self):
+        _, _, constants, _ = initialize_model(
+            self._index("1H"), 2103.0, z_u=7.5, z_t=3.0, z_g=0.5
+        )
+        assert constants["z_u"] == 7.5
+        assert constants["z_t"] == 3.0
+        assert constants["z_g"] == 0.5
+
+    def test_measurement_height_defaults(self):
+        _, _, constants, _ = initialize_model(self._index("1H"), 2103.0)
+        assert (constants["z_u"], constants["z_t"], constants["z_g"]) == (
+            5.0, 2.0, 0.3
+        )
+
+    def test_tstep_info_has_no_output_flag(self):
+        # snobal's built-in output paths are disabled; the flag is not used.
+        _, tstep_info, _, _ = initialize_model(self._index("1H"), 2103.0)
+        assert all("output" not in level for level in tstep_info)
+
+    def test_sub_hourly_frequency_rejected(self):
+        with pytest.raises(ValueError, match="whole-hour timestep"):
+            initialize_model(self._index("30min"), 2103.0)
+
+    def test_irregular_index_rejected(self):
+        index = self._index("1H").delete(5)  # break the regular spacing
+        with pytest.raises(ValueError, match="regular frequency"):
+            initialize_model(index, 2103.0)
 
 
 class TestRunSnobal:
     TEST_FILE = Path(__file__).parent.joinpath(
         "data/inputs_csl_2023.csv"
     )
+    ELEVATION = 2103.0
 
     @pytest.fixture(scope="class")
     def test_data(self):
@@ -17,15 +86,76 @@ def test_data(self):
             parse_dates=["datetime"], index_col="datetime"
         )
 
-    def test_run_snobal(self, test_data):
-        start_date = test_data.index.min()
-        end_date = test_data.index.max()
+    @pytest.fixture(scope="class")
+    def daily_result(self, test_data):
+        return run_model(self.ELEVATION, test_data)
+
+    def test_run_snobal(self, daily_result):
+        # 6H input, daily output. The specific_mass gold reflects the corrected
+        # sub-timestep integration (one normal step per data hour); a loose
+        # tolerance keeps it robust across compiler/architecture float drift.
+        assert len(daily_result) == 302
+        assert daily_result["specific_mass"].values[200] == pytest.approx(
+            1548.35, rel=1e-3
+        )
+
+    def test_snowpack_accumulates_and_melts_out(self, daily_result):
+        mass = daily_result["specific_mass"]
+        assert (mass >= 0).all()
+        assert mass.max() > 500  # a real snowpack built up
+        assert mass.iloc[-1] == pytest.approx(0.0, abs=1e-6)  # melted out
 
-        # Run the model for that file
-        result = run_model(
-            start_date, end_date, 2103.0, test_data
+    def test_measurement_heights_change_results(self, test_data, daily_result):
+        # Different measurement heights must actually reach the model and
+        # change the turbulent/soil fluxes -> a different snowpack.
+        altered = run_model(
+            self.ELEVATION, test_data, z_u=2.0, z_t=1.0, z_g=0.1
         )
-        assert len(result) == 302
-        assert result["specific_mass"].values[200] == pytest.approx(
-            1623.4878530514477
+        # Lower/closer sensors change the turbulent + soil fluxes, giving a
+        # different modeled snowpack than the default heights.
+        assert altered["specific_mass"].values[200] == pytest.approx(
+            1476.46, rel=1e-3
+        )
+        assert not daily_result["specific_mass"].equals(
+            altered["specific_mass"]
+        )
+
+    def test_hourly_output_yields_more_rows(self, test_data, daily_result):
+        sub_daily = run_model(
+            self.ELEVATION, test_data, output_frequency_hours=1
         )
+        # 6H input -> the finest available cadence is every data step
+        assert len(sub_daily) == 1208
+
+    def test_output_cadence_preserves_instantaneous_state(
+        self, test_data, daily_result
+    ):
+        # The output cadence only changes which timesteps are recorded, not
+        # the physics: instantaneous state must match at shared timestamps.
+        sub_daily = run_model(
+            self.ELEVATION, test_data, output_frequency_hours=1
+        )
+        shared = daily_result.index.intersection(sub_daily.index)
+        assert len(shared) == len(daily_result)
+        for column in ["thickness", "specific_mass", "snow_density",
+                       "temp_snowcover"]:
+            diff = (
+                daily_result.loc[shared, column]
+                - sub_daily.loc[shared, column]
+            ).abs().max()
+            assert diff == pytest.approx(0.0, abs=1e-9)
+
+    def test_swi_conserved_across_output_cadence(self, test_data, daily_result):
+        # SWI is summed over each output interval, so the season total is
+        # independent of how finely it is reported.
+        sub_daily = run_model(
+            self.ELEVATION, test_data, output_frequency_hours=1
+        )
+        assert sub_daily["SWI"].sum() == pytest.approx(
+            daily_result["SWI"].sum(), rel=1e-9
+        )
+
+    def test_run_model_rejects_irregular_input(self, test_data):
+        broken = test_data.drop(test_data.index[10])
+        with pytest.raises(ValueError):
+            run_model(self.ELEVATION, broken)