From c7937ea0ed270c714fcf10dffb375f76c2c7f7ec Mon Sep 17 00:00:00 2001 From: Micah Sandusky Date: Mon, 31 Aug 2026 13:59:11 -0600 Subject: [PATCH] Allowing hourly output. Allow height setting pass in --- pointsnobal/c_snobal/snobal.pyx | 6 +- pointsnobal/cli.py | 9 +- pointsnobal/point_model.py | 83 ++++++++++++------ tests/test_run_snobal.py | 150 +++++++++++++++++++++++++++++--- 4 files changed, 203 insertions(+), 45 deletions(-) diff --git a/pointsnobal/c_snobal/snobal.pyx b/pointsnobal/c_snobal/snobal.pyx index 6f8399e..9e582bc 100644 --- a/pointsnobal/c_snobal/snobal.pyx +++ b/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. diff --git a/pointsnobal/cli.py b/pointsnobal/cli.py index cc81a91..c4f4670 100644 --- a/pointsnobal/cli.py +++ b/pointsnobal/cli.py @@ -32,6 +32,10 @@ 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}") @@ -39,16 +43,13 @@ def main(): 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) diff --git a/pointsnobal/point_model.py b/pointsnobal/point_model.py index da4fdf6..d1181af 100644 --- a/pointsnobal/point_model.py +++ b/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)) diff --git a/tests/test_run_snobal.py b/tests/test_run_snobal.py index c4aae9b..4121fe6 100644 --- a/tests/test_run_snobal.py +++ b/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)