Write the PCRaster header in the time series outputs

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 <https://github.com/LabSid-USP/RUBEM/issues/347>`__).
 - The impervious area interception ``i_imp`` must lie between 1 and 3 mm, the range of the
   published formulation (`#327 <https://github.com/LabSid-USP/RUBEM/issues/327>`__);
   the model overview and the user guide state the same range.

doc/source/fileformats.rst

@@ -763,6 +763,36 @@ Resulting maps of Accumulated Total Runoff [:raw-html:`m<sup>3</sup>s<sup>-1</su
   - Columns = :ref:`clone columns<fileformats:Mask of Catchment (Clone) raster>`;
   - Cell Size = :ref:`clone cell size<fileformats:Mask of Catchment (Clone) raster>`.
 
+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 <fileformats:Stations (samples) raster>`
+   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
 ^^^^^^^^^^^^^^^^^^^^^^^^^
 

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:

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;<id>``
+    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(

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))

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))