diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index c7ddaca4..ca4b485e 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -95,6 +95,12 @@ Added Changed ``````` +- The station time series the model keeps as PCRaster ``.tss`` files carry the + PCRaster header (title, number of columns, ``timestep`` line and one line per + station id), so PCRaster's own tools read them; the CSV conversion reads the + header, refuses a file without it or whose ids differ from the configured + stations, and converts the data rows only, leaving the CSV tables unchanged + (`#347 `__). - The impervious area interception ``i_imp`` must lie between 1 and 3 mm, the range of the published formulation (`#327 `__); the model overview and the user guide state the same range. diff --git a/doc/source/fileformats.rst b/doc/source/fileformats.rst index 4c8f0fe0..17aa3ce7 100644 --- a/doc/source/fileformats.rst +++ b/doc/source/fileformats.rst @@ -763,6 +763,36 @@ Resulting maps of Accumulated Total Runoff [:raw-html:`m3s-1`; - Cell Size = :ref:`clone cell size`. +Time series of the sampling stations +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The time series of the sampling stations are written by PCRaster as +:file:`*.tss` files and can be kept as such, converted to Comma-Separated +Values (CSV) :file:`*.csv` files, or both, according to the time series formats +of the results specification. + +A kept :file:`*.tss` file carries the PCRaster header: + +#. the title line, ``timeseries scalar``; +#. the number of columns, the time step column included; +#. the line ``timestep``; +#. one line per sampling station, with its identifier: the identifier of the + station in the :ref:`stations map ` + for the ``point`` and ``subcatchment`` aggregations, or the renumbered zone + column (``1`` to the number of zones, whose correspondence with the zone + identifiers is written to :file:`zones_mapping.csv`) for the ``zones`` + aggregation. + +The data rows follow the header, one per time step: the step number and one +value per station, in the order of the identifiers above. Because the header is +present, these files are read by PCRaster's own tools, Aguila included. + +The CSV of each variable is derived from the data rows of the corresponding +:file:`*.tss` file: its first line names the columns (``0`` for the time step, +then the station identifiers) and each following line is one data row of the +time series, with the values as PCRaster wrote them. The tables described below +are these CSV files. + Total Interception table ^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/rubem/_dynamic_model.py b/rubem/_dynamic_model.py index 77f8cf7c..fb90be07 100644 --- a/rubem/_dynamic_model.py +++ b/rubem/_dynamic_model.py @@ -590,6 +590,9 @@ def __initial_setup_timeoutput_timeseries(self): """Initial setup of timeoutput timeseries. Initialize Tss report at sample locations or pits for each enabled output variable. + Every file is written with the PCRaster header (title, number of + columns, ``timestep`` line and one line per station id), which is the + form PCRaster's own tools read. """ point_map = self.config.output_variables.aggregation == "point" and not is_geotiff( self.config.raster_files.sample_locations @@ -600,7 +603,7 @@ def __initial_setup_timeoutput_timeseries(self): str(Path(self.config.output_directory.path) / var.table_filename_prefix), self, id_map, - noHeader=True, + noHeader=False, ) self.sample_time_series_dict[var.id] = tss_file.sample if point_map: diff --git a/rubem/file/_file_conversions.py b/rubem/file/_file_conversions.py index e55c9e20..793dba9f 100644 --- a/rubem/file/_file_conversions.py +++ b/rubem/file/_file_conversions.py @@ -159,6 +159,75 @@ def _undo_installs(installs: list[tuple[str, str | None]]) -> None: logger.error("Error while restoring file %s. %s", dst_file_path, e) +def _is_numeric(token: str) -> bool: + """Whether ``token`` reads as a number, which a title line never does.""" + try: + float(token) + except ValueError: + return False + return True + + +def _split_header(lines: list[str], tss_file: str, cols_names: list[str]) -> list[str]: + """Return the data rows of ``lines``, after checking the PCRaster header. + + The header the model writes is the one :class:`TimeoutputTimeseries` + produces: a title line (``timeseries scalar``), the number of columns (the + time step column included), a ``timestep`` line and one line per station + id. The ids are checked against ``cols_names`` so that a file written for + other stations is never converted under the wrong column names. + + :raises ValueError: If the file carries no header, or if its header does + not describe the given column names. + """ + stripped = [line.strip() for line in lines] + headerless = ( + f"The time series file {tss_file} has no PCRaster header; the model writes " + "headed time series files." + ) + if not stripped or not stripped[0] or _is_numeric(stripped[0].split()[0]): + logger.error("The time series file %s has no PCRaster header.", tss_file) + raise ValueError(headerless) + try: + number_of_columns = int(stripped[1]) + except (IndexError, ValueError): + logger.error("The time series file %s has no PCRaster header.", tss_file) + raise ValueError(headerless) from None + # A header always announces at least the time step column, so a count + # below one is not a column count and the second line is not a header one. + if number_of_columns < 1: + logger.error("The time series file %s has no PCRaster header.", tss_file) + raise ValueError(headerless) + + ids = stripped[3 : 3 + number_of_columns - 1] + # A data row where an id line is expected means the header stops short of + # the announced column count. + if len(stripped) < 3 or len(ids) != number_of_columns - 1 or any(" " in id_ for id_ in ids): + logger.error("The time series file %s has a truncated PCRaster header.", tss_file) + raise ValueError( + f"The PCRaster header of the time series file {tss_file} is truncated: it " + f"announces {number_of_columns} column(s) but does not list " + f"{number_of_columns - 1} station id(s)." + ) + if number_of_columns - 1 != len(cols_names): + logger.error( + "The number of columns in the header of %s is different from the number of " + "column names.", + tss_file, + ) + raise ValueError( + f"The PCRaster header of the time series file {tss_file} announces " + f"{number_of_columns} column(s), the column names give {len(cols_names) + 1}." + ) + if ids != list(cols_names): + logger.error("The station ids in the header of %s are not the configured ones.", tss_file) + raise ValueError( + f"The station ids in the PCRaster header of the time series file {tss_file} " + f"are {ids}, the column names give {list(cols_names)}." + ) + return lines[3 + number_of_columns - 1 :] + + def tss2csv(tss_files, cols_names: list[str], should_delete_src_tss: bool = True) -> None: """Convert the given PCRaster Time Series (``*.tss``) files to ``*.csv``. @@ -176,15 +245,26 @@ def tss2csv(tss_files, cols_names: list[str], should_delete_src_tss: bool = True backup names are allocated so that unrelated files next to a destination are never overwritten. - :raises ValueError: If the column names are empty, if a source has no data - rows, or if a source has a different number of columns. + Every source must carry the PCRaster header the model writes (a title + line, the number of columns, the ``timestep`` line and one line per + station id): the header is checked against the column names given here and + only the data rows below it are converted, so the CSV keeps its ``0;`` + header and one row per time step. + + :raises ValueError: If the column names are empty, if a source carries no + PCRaster header, if its header does not describe the given column + names, if a source has no data rows, or if a source has a different + number of columns. :raises IsADirectoryError: If a destination path is an existing directory. :param tss_files: Paths of the ``.tss`` files to convert. Each item accepts anything :func:`rubem._paths.as_path` does, ``bytes`` included. :type tss_files: list - :param cols_names: List of strings of aliases for the column names. + :param cols_names: Station ids of the columns, in the order the time + series holds them. They name the columns of the CSV and are checked + against the station ids the PCRaster header of each source lists, so + they are those ids and not free aliases. :type cols_names: list[str] :param should_delete_src_tss: Remove the source files after conversion, defaults to ``True``. @@ -210,7 +290,8 @@ def tss2csv(tss_files, cols_names: list[str], should_delete_src_tss: bool = True with Path(tss_file).open(encoding="utf8") as f: lines = f.readlines() - data = [line.split() for line in lines if line.strip()] + rows = _split_header(lines, tss_file, cols_names) + data = [line.split() for line in rows if line.strip()] if not data: logger.error("The time series file %s has no data rows.", tss_file) raise ValueError( diff --git a/tests/unit/core/test_core.py b/tests/unit/core/test_core.py index ef347023..6ea22f65 100644 --- a/tests/unit/core/test_core.py +++ b/tests/unit/core/test_core.py @@ -78,6 +78,16 @@ def test_outputs_are_finite_on_valid_cells(self, tmp_path): assert valid.size assert np.isfinite(valid).all(), f"{name} has non-finite cells" + @pytest.mark.unit + def test_the_csv_tables_carry_only_the_column_names_and_the_data_rows(self, tmp_path): + """The PCRaster header of the sources never reaches the converted table.""" + run_model(str(tmp_path)) + lines = (tmp_path / "out" / "tss_arn.csv").read_text(encoding="utf8").splitlines() + + assert lines[0] == "0;1;2" + assert len(lines) == 3 + assert [line.split(";")[0] for line in lines[1:]] == ["1", "2"] + @pytest.mark.unit def test_time_series_cover_every_station_and_step(self, tmp_path): run_model(str(tmp_path)) diff --git a/tests/unit/core/test_dynamic_model.py b/tests/unit/core/test_dynamic_model.py index d28d7887..1c77a7e9 100644 --- a/tests/unit/core/test_dynamic_model.py +++ b/tests/unit/core/test_dynamic_model.py @@ -176,6 +176,28 @@ def test_stale_tss_files_in_the_output_directory_are_left_alone(self, tmp_path): assert stale.exists() assert not (tmp_path / "out" / "stale.csv").exists() + @pytest.mark.unit + def test_a_kept_tss_file_carries_the_pcraster_header(self, tmp_path): + """The kept time series are in the form PCRaster's own tools read.""" + from rubem.configuration.model_configuration_file import ModelConfigurationFile + from rubem.configuration.model_configuration_file_v1 import ModelConfigurationFileV1 + + legacy = ModelConfigurationFile.model_validate(write_synthetic_dataset(str(tmp_path))) + document = ModelConfigurationFileV1.from_legacy(legacy).to_dict() + document["model_simulation_output"]["time_series_samples"]["formats"] = ["PCRasterTSS"] + + DynamicFrameworkWrapper.load(ModelConfiguration(document)).run() + + tss = tmp_path / "out" / "tss_arn.tss" + lines = tss.read_text(encoding="utf8").splitlines() + # Title, number of columns (the time step column included), the + # ``timestep`` line and one line per station id of the dataset. + assert lines[:4] == ["timeseries scalar", "3", "timestep", "1"] + assert lines[3:5] == ["1", "2"] + data = [line.split() for line in lines[5:] if line.strip()] + assert [row[0] for row in data] == ["1", "2"] + assert all(len(row) == 3 for row in data) + @pytest.mark.unit def test_missing_first_ndvi_step_raises_a_clear_error(self, tmp_path): config = write_synthetic_dataset(str(tmp_path)) diff --git a/tests/unit/file/test_file_conversions.py b/tests/unit/file/test_file_conversions.py index 4b1c03f2..c2266f49 100644 --- a/tests/unit/file/test_file_conversions.py +++ b/tests/unit/file/test_file_conversions.py @@ -7,6 +7,16 @@ from rubem.file._file_conversions import _reserve_stage, tss2csv +def _read_exact(path): + """Return the text of ``path`` without newline translation. + + ``Path.read_text`` accepts ``newline`` only from Python 3.13 on; the CSV + writer emits CRLF and the comparison must see it. + """ + with path.open(encoding="utf8", newline="") as file: + return file.read() + + class FakeUUID: """A stand-in for ``uuid.uuid4()`` results with a fixed ``hex`` value.""" @@ -14,7 +24,28 @@ def __init__(self, hex_value): self.hex = hex_value -def write_tss(path, rows): +def write_tss(path, rows, ids=None, number_of_columns=None, title="timeseries scalar"): + """Write a ``.tss`` file in the headed form the model produces. + + ``ids`` defaults to the station ids the first row has values for, and + ``number_of_columns`` to the column count that header describes; both are + given explicitly by the tests that check a header disagreeing with the + column names. + """ + if ids is None: + ids = [str(column) for column in range(1, len(rows[0]))] + if number_of_columns is None: + number_of_columns = len(ids) + 1 + with open(path, "w", encoding="utf8") as f: + f.write(f"{title}\n{number_of_columns}\ntimestep\n") + for station in ids: + f.write(f"{station}\n") + for row in rows: + f.write(" ".join(str(value) for value in row) + "\n") + + +def write_headerless_tss(path, rows): + """Write the headerless form the model used to produce.""" with open(path, "w", encoding="utf8") as f: for row in rows: f.write(" ".join(str(value) for value in row) + "\n") @@ -118,12 +149,120 @@ def test_a_later_row_with_a_divergent_column_count_is_refused(self, tmp_path): assert not list(tmp_path.glob("*.csv")) assert not list(tmp_path.glob("*.tmp")) + @pytest.mark.unit + def test_a_headerless_source_is_refused(self, tmp_path): + """The model writes headed files; a headerless one is not converted.""" + tss = tmp_path / "tss_itp.tss" + write_headerless_tss(tss, [(1, 10.5), (2, 11.0)]) + + with pytest.raises(ValueError, match="has no PCRaster header"): + tss2csv([tss], ["1"]) + + assert tss.exists() + assert not list(tmp_path.glob("*.csv")) + assert not list(tmp_path.glob("*.tmp")) + + @pytest.mark.unit + def test_a_source_whose_second_line_is_not_an_integer_is_refused(self, tmp_path): + tss = tmp_path / "tss_itp.tss" + tss.write_text("timeseries scalar\nnot a number\ntimestep\n1\n 1 10.5\n", encoding="utf8") + + with pytest.raises(ValueError, match="has no PCRaster header"): + tss2csv([tss], ["1"]) + + assert not list(tmp_path.glob("*.csv")) + + @pytest.mark.unit + def test_a_source_announcing_no_column_is_refused(self, tmp_path): + """A count below one is not a column count: the time step is a column.""" + tss = tmp_path / "tss_itp.tss" + tss.write_text("timeseries scalar\n-3\ntimestep\n1\n 1 10.5\n", encoding="utf8") + + with pytest.raises(ValueError, match="has no PCRaster header"): + tss2csv([tss], ["1"]) + + assert tss.exists() + assert not list(tmp_path.glob("*.csv")) + + @pytest.mark.unit + def test_a_zero_byte_source_is_refused_and_keeps_the_previous_csv(self, tmp_path): + """A run cut short before any line must not replace usable results.""" + first = tmp_path / "tss_a.tss" + second = tmp_path / "tss_b.tss" + write_tss(first, [(1, 10.5)]) + second.write_text("", encoding="utf8") + (tmp_path / "tss_a.csv").write_text("previous a\n", encoding="utf8") + (tmp_path / "tss_b.csv").write_text("previous b\n", encoding="utf8") + + with pytest.raises(ValueError, match="has no PCRaster header"): + tss2csv([first, second], ["1"]) + + assert (tmp_path / "tss_a.csv").read_text(encoding="utf8") == "previous a\n" + assert (tmp_path / "tss_b.csv").read_text(encoding="utf8") == "previous b\n" + assert sorted(path.name for path in tmp_path.iterdir()) == [ + "tss_a.csv", + "tss_a.tss", + "tss_b.csv", + "tss_b.tss", + ] + + @pytest.mark.unit + def test_a_header_with_other_station_ids_is_refused(self, tmp_path): + tss = tmp_path / "tss_itp.tss" + write_tss(tss, [(1, 10.5, 20.5)], ids=["3", "4"]) + + with pytest.raises(ValueError, match=r"are \['3', '4'\].*give \['1', '2'\]"): + tss2csv([tss], ["1", "2"]) + + assert tss.exists() + assert not list(tmp_path.glob("*.csv")) + + @pytest.mark.unit + def test_a_header_with_another_column_count_is_refused(self, tmp_path): + tss = tmp_path / "tss_itp.tss" + write_tss(tss, [(1, 10.5, 20.5)], ids=["1", "2"]) + + with pytest.raises(ValueError, match=r"announces 3 column\(s\), the column names give 2"): + tss2csv([tss], ["1"]) + + assert tss.exists() + assert not list(tmp_path.glob("*.csv")) + + @pytest.mark.unit + def test_a_truncated_header_is_refused(self, tmp_path): + """A header that stops short of the column count it announces.""" + ending = tmp_path / "tss_ending.tss" + with_rows = tmp_path / "tss_with_rows.tss" + write_tss(ending, [], ids=["1"], number_of_columns=3) + write_tss(with_rows, [(1, 10.5, 20.5)], ids=["1"], number_of_columns=3) + + for tss in (ending, with_rows): + with pytest.raises(ValueError, match="is truncated"): + tss2csv([tss], ["1", "2"]) + + assert tss.exists() + assert not list(tmp_path.glob("*.csv")) + + @pytest.mark.unit + def test_the_csv_of_a_headed_file_is_the_one_the_headerless_form_produced(self, tmp_path): + """Only the data rows are converted, so the CSV text is byte for byte + the one the previous headerless sources produced from the same rows.""" + rows = [(1, 10.5, 20.5), (2, 11.0, 21.0)] + tss = tmp_path / "tss_itp.tss" + write_tss(tss, rows) + + tss2csv([tss], ["1", "2"]) + + produced = _read_exact(tmp_path / "tss_itp.csv") + assert produced == "0;1;2\r\n1;10.5;20.5\r\n2;11.0;21.0\r\n" + @pytest.mark.unit def test_failure_leaves_sources_and_no_partial_outputs(self, tmp_path): good = tmp_path / "tss_good.tss" bad = tmp_path / "tss_bad.tss" write_tss(good, [(1, 10.5)]) - write_tss(bad, [(1, 10.5, 99.0)]) + # The header agrees with the column names; the data row does not. + write_tss(bad, [(1, 10.5, 99.0)], ids=["1"]) with pytest.raises(ValueError): tss2csv([good, bad], ["1"]) @@ -271,7 +410,7 @@ def test_empty_source_fails_and_keeps_the_previous_csv(self, tmp_path): first = tmp_path / "tss_a.tss" second = tmp_path / "tss_b.tss" write_tss(first, [(1, 10.5)]) - second.write_text("", encoding="utf8") + write_tss(second, [], ids=["1"]) (tmp_path / "tss_a.csv").write_text("previous a\n", encoding="utf8") (tmp_path / "tss_b.csv").write_text("previous b\n", encoding="utf8")