Add wavelength attribute to workflow to convert rasters to meters (#…

…367)

* Simplify `print_yaml_schema`, V2 seems to have fixed problems

* add `input_options.wavelength`, auto-recognize OPERA CSLC-S1

* add `wavelength` as arguments in `timeseries` functions

* add `units` and `description` options to writing rasters

* add `constants` module

* Add `_convert_and_reference` to timeseries instead of just copying

* fix output units for full workflow, fix single-reference check

* fix the sign flp during `invert_unw_network`

src/dolphin/constants.py

@@ -0,0 +1,4 @@
+SPEED_OF_LIGHT = 299_792_458  # meters / second
+
+SENTINEL_1_FREQUENCY = 5.405e9  # Hz
+SENTINEL_1_WAVELENGTH = SPEED_OF_LIGHT / SENTINEL_1_FREQUENCY  # meters

src/dolphin/timeseries.py

@@ -1,7 +1,6 @@
 from __future__ import annotations
 
 import logging
-import shutil
 from enum import Enum
 from pathlib import Path
 from tempfile import NamedTemporaryFile
@@ -52,6 +51,7 @@ def run(
     block_shape: tuple[int, int] = (256, 256),
     num_threads: int = 4,
     reference_point: tuple[int, int] = (-1, -1),
+    wavelength: float | None = None,
 ) -> tuple[list[Path], ReferencePoint]:
     """Invert the unwrapped interferograms, estimate timeseries and phase velocity.
 
@@ -92,6 +92,12 @@ def run(
         Reference point (row, col) used if performing a time series inversion.
         If not provided, a point will be selected from a consistent connected
         component with low amplitude dispersion or high temporal coherence.
+    wavelength : float, optional
+        The wavelength of the radar signal, in meters.
+        If provided, the output rasters are in meters and meters / year for
+        the displacement and velocity rasters.
+        If not provided, the outputs are in radians.
+        See Notes for line of sight convention.
 
     Returns
     -------
@@ -102,6 +108,12 @@ def run(
         If passed as input, simply returned back as output.
         Otherwise, the result is the auto-selection from `select_reference_point`.
 
+    Notes
+    -----
+    When wavelength is provided, the output rasters are in meters and meters / year,
+    where positive values indicate motion *toward* from the radar (i.e. positive values
+    in both ascending and descending tracks imply uplift).
+
     """
     Path(output_dir).mkdir(exist_ok=True, parents=True)
 
@@ -120,10 +132,24 @@ def run(
     sar_dates = sorted(set(utils.flatten(ifg_date_pairs)))
     # if we did single-reference interferograms, for `n` sar dates, we will only have
     # `n-1` interferograms. Any more than n-1 ifgs means we need to invert
-    needs_inversion = len(unwrapped_paths) > len(sar_dates) - 1
+    is_single_reference = (len(unwrapped_paths) == len(sar_dates) - 1) and all(
+        pair[0] == ifg_date_pairs[0][0] for pair in ifg_date_pairs
+    )
+
     # check if we even need to invert, or if it was single reference
     inverted_phase_paths: list[Path] = []
-    if needs_inversion:
+    if is_single_reference:
+        logger.info(
+            "Skipping inversion step: only single reference interferograms exist."
+        )
+        # Copy over the unwrapped paths to `timeseries/`
+        inverted_phase_paths = _convert_and_reference(
+            unwrapped_paths,
+            output_dir=output_dir,
+            reference_point=ref_point,
+            wavelength=wavelength,
+        )
+    else:
         logger.info("Selecting a reference point for unwrapped interferograms")
 
         logger.info("Inverting network of %s unwrapped ifgs", len(unwrapped_paths))
@@ -133,23 +159,9 @@ def run(
             output_dir=output_dir,
             block_shape=block_shape,
             num_threads=num_threads,
+            wavelength=wavelength,
             method=method,
         )
-    else:
-        logger.info(
-            "Skipping inversion step: only single reference interferograms exist."
-        )
-        # Copy over the unwrapped paths to `timeseries/`
-        for p in unwrapped_paths:
-            # if it ends in `.unw.tif`, change to just `.tif` for consistency
-            # with the case where we run an inversion
-            cur_name = Path(p).name
-            unw_suffix = full_suffix(p)
-            target_name = str(cur_name).replace(unw_suffix, ".tif")
-            target = Path(output_dir) / target_name
-            if not target.exists():  # Check to prevent overwriting
-                shutil.copy(p, target)
-            inverted_phase_paths.append(target)
 
     if run_velocity:
         #  We can't pass the correlations after an inversion- the numbers don't match
@@ -175,6 +187,51 @@ def run(
     return inverted_phase_paths, ref_point
 
 
+def _convert_and_reference(
+    unwrapped_paths: Sequence[PathOrStr],
+    *,
+    output_dir: PathOrStr,
+    reference_point: ReferencePoint,
+    wavelength: float | None = None,
+) -> list[Path]:
+    if wavelength is not None:
+        # Positive values are motion towards the radar
+        constant = -1 * (wavelength / (4 * np.pi))
+        units = "meters"
+    else:
+        constant = -1
+        units = "radians"
+
+    ref_row, ref_col = reference_point
+    out_paths: list[Path] = []
+    for p in unwrapped_paths:
+        # if it ends in `.unw.tif`, change to just `.tif` for consistency
+        # with the case where we run an inversion
+        cur_name = Path(p).name
+        unw_suffix = full_suffix(p)
+        target_name = str(cur_name).replace(unw_suffix, ".tif")
+        target = Path(output_dir) / target_name
+        out_paths.append(target)
+
+        if target.exists():  # Check to prevent overwriting
+            continue
+
+        arr_radians = io.load_gdal(p)
+        # Reference to the
+        ref_value = arr_radians[ref_row, ref_col]
+        if np.isnan(ref_value):
+            logger.warning(
+                "{ref_point!r} is NaN for {p} . Skipping reference subtraction."
+            )
+        else:
+            arr_radians -= ref_value
+        io.write_arr(
+            arr=arr_radians * constant, output_name=target, units=units, like_filename=p
+        )
+
+    return out_paths
+
+
 def argmin_index(arr: ArrayLike) -> tuple[int, ...]:
     """Get the index tuple of the minimum value of the array.
 
@@ -576,7 +633,6 @@ def read_and_fit(
     readers = [unw_reader]
     if cor_reader is not None:
         readers.append(cor_reader)
-
     writer = io.BackgroundRasterWriter(output_file, like_filename=unw_file_list[0])
     io.process_blocks(
         readers=readers,
@@ -590,6 +646,8 @@ def read_and_fit(
     if add_overviews:
         logger.info("Creating overviews for velocity image")
         create_overviews([output_file])
+    if units := io.get_raster_units(unw_file_list[0]):
+        io.set_raster_units(output_file, units=units)
     logger.info("Completed create_velocity")
 
 
@@ -664,6 +722,7 @@ def invert_unw_network(
     cor_threshold: float = 0.2,
     n_cor_looks: int = 1,
     ifg_date_pairs: Sequence[Sequence[DateOrDatetime]] | None = None,
+    wavelength: float | None = None,
     method: InversionMethod = InversionMethod.L2,
     block_shape: tuple[int, int] = (256, 256),
     num_threads: int = 4,
@@ -681,15 +740,6 @@ def invert_unw_network(
         from all other points when solving.
     output_dir : PathOrStr
         The directory to save the output files
-    ifg_date_pairs : Sequence[Sequence[DateOrDatetime]], optional
-        List of date pairs to use for the inversion. If not provided, will be
-        parsed from filenames in `unw_file_list`.
-    method : str, choices = "L1", "L2"
-        Inversion method to use when solving Ax = b.
-        Default is L2, which uses least squares to solve Ax = b (faster).
-        "L1" minimizes |Ax - b|_1 at each pixel.
-    block_shape : tuple[int, int], optional
-        The shape of the blocks to process in parallel
     cor_file_list : Sequence[PathOrStr], optional
         List of correlation files to use for weighting the inversion
     cor_threshold : float, optional
@@ -699,6 +749,20 @@ def invert_unw_network(
         The number of looks used to form the input correlation data, used
         to convert correlation to phase variance.
         Default is 1.
+    ifg_date_pairs : Sequence[Sequence[DateOrDatetime]], optional
+        List of date pairs to use for the inversion. If not provided, will be
+        parsed from filenames in `unw_file_list`.
+    method : str, choices = "L1", "L2"
+        Inversion method to use when solving Ax = b.
+        Default is L2, which uses least squares to solve Ax = b (faster).
+        "L1" minimizes |Ax - b|_1 at each pixel.
+    wavelength : float, optional
+        The wavelength of the radar signal, in meters.
+        If provided, the output rasters are in meters.
+        If not provided, the outputs are in radians.
+    block_shape : tuple[int, int], optional
+        The shape of the blocks to process in parallel.
+        Default is (256, 256).
     num_threads : int
         The parallel blocks to process at once.
         Default is 4.
@@ -747,6 +811,14 @@ def invert_unw_network(
     ref_row, ref_col = reference
     ref_data = unw_reader[:, ref_row, ref_col].reshape(-1, 1, 1)
 
+    if wavelength is not None:
+        # Positive values are motion towards the radar
+        constant = -1 * (wavelength / (4 * np.pi))
+        units = "meters"
+    else:
+        constant = -1
+        units = "radians"
+
     def read_and_solve(
         readers: Sequence[io.StackReader], rows: slice, cols: slice
     ) -> tuple[slice, slice, np.ndarray]:
@@ -771,7 +843,8 @@ def read_and_solve(
             phases = invert_stack_l1(A, stack)[0]
         else:
             phases = invert_stack(A, stack, weights)[0]
-        return np.asarray(phases), rows, cols
+        # Convert to meters, with LOS convention:
+        return constant * np.asarray(phases), rows, cols
 
     if cor_file_list is not None:
         cor_reader = io.VRTStack(
@@ -783,7 +856,9 @@ def read_and_solve(
         readers = [unw_reader]
         logger.info("Using unweighted unw inversion")
 
-    writer = io.BackgroundStackWriter(out_paths, like_filename=unw_file_list[0])
+    writer = io.BackgroundStackWriter(
+        out_paths, like_filename=unw_file_list[0], units=units
+    )
 
     io.process_blocks(
         readers=readers,
@@ -795,7 +870,7 @@ def read_and_solve(
     writer.notify_finished()
 
     if add_overviews:
-        logger.info("Creating overviews for unwrapped images")
+        logger.info("Creating overviews for timeseries images")
         create_overviews(out_paths, image_type=ImageType.UNWRAPPED)
 
     logger.info("Completed invert_unw_network")

src/dolphin/unwrap/_unwrap.py

@@ -122,6 +122,8 @@ def run(
             options=unwrap_options.spurt_options,
             scratchdir=scratchdir,
         )
+        for f in unw_paths:
+            io.set_raster_units(f, "radians")
         return unw_paths, conncomp_paths
 
     ifg_suffixes = [full_suffix(f) for f in ifg_filenames]
@@ -190,6 +192,8 @@ def run(
             Path(str(outf).replace(UNW_SUFFIX, CONNCOMP_SUFFIX))
             for outf in all_out_files
         ]
+    for f in all_out_files:
+        io.set_raster_units(f, "radians")
     return all_out_files, conncomp_files
 
 

src/dolphin/workflows/config/_common.py

@@ -258,6 +258,14 @@ class InputOptions(BaseModel, extra="forbid"):
         "%Y%m%d",
         description="Format of dates contained in CSLC filenames",
     )
+    wavelength: Optional[float] = Field(
+        None,
+        description=(
+            "Radar wavelength (in meters) of the transmitted data. used to convert the"
+            " units in the rasters in `timeseries/` to from radians to meters. If None"
+            " and sensor is not recognized, outputs remain in radians."
+        ),
+    )
 
 
 class OutputOptions(BaseModel, extra="forbid"):

src/dolphin/workflows/config/_displacement.py

@@ -4,7 +4,7 @@
 from pathlib import Path
 from typing import Annotated, Any, Optional
 
-from opera_utils import get_dates, sort_files_by_date
+from opera_utils import get_burst_id, get_dates, sort_files_by_date
 from pydantic import (
     BaseModel,
     ConfigDict,
@@ -14,6 +14,7 @@
     model_validator,
 )
 
+from dolphin import constants
 from dolphin._types import TropoModel, TropoType
 
 from ._common import (
@@ -187,6 +188,16 @@ def _check_input_files_exist(self) -> DisplacementWorkflow:
     def model_post_init(self, __context: Any) -> None:
         """After validation, set up properties for use during workflow run."""
         super().model_post_init(__context)
+
+        if self.input_options.wavelength is None and self.cslc_file_list:
+            # Try to infer the wavelength from filenames
+            try:
+                get_burst_id(self.cslc_file_list[-1])
+                # The Burst ID was recognized for OPERA-S1 SLCs: use S1 wavelength
+                self.input_options.wavelength = constants.SENTINEL_1_WAVELENGTH
+            except ValueError:
+                pass
+
         # Ensure outputs from workflow steps are within work directory.
         if not self.keep_paths_relative:
             # Resolve all CSLC paths:

src/dolphin/workflows/config/_yaml_model.py

@@ -1,9 +1,7 @@
 import json
 import sys
 import textwrap
-import warnings
 from io import StringIO
-from itertools import repeat
 from typing import Optional, TextIO, Union
 
 from pydantic import BaseModel
@@ -106,32 +104,9 @@ def print_yaml_schema(
             Number of spaces to indent per level.
 
         """
-        full_dict = cls._construct_empty()
-        # UserWarning: Pydantic serializer warnings:
-        # New V2 warning, but seems harmless for just printing the schema
-        with warnings.catch_warnings():
-            warnings.filterwarnings("ignore", category=UserWarning)
-            cls.model_construct(**full_dict).to_yaml(
-                output_path, with_comments=True, indent_per_level=indent_per_level
-            )
-
-    @classmethod
-    def _construct_empty(cls):
-        """Construct a model with all fields filled in.
-
-        Uses defaults, or `None` for required fields.
-
-        The .construct is a pydantic method to disable validation
-        https://docs.pydantic.dev/usage/models/#creating-models-without-validation
-        But Fields without a default value don't get filled in
-        so first, we manually make a dict of all the fields with None values
-        then we update it with the default-filled values
-        """
-        all_none_vals = dict(
-            zip(cls.model_json_schema()["properties"].keys(), repeat(None))
+        cls.model_construct().to_yaml(
+            output_path, with_comments=True, indent_per_level=indent_per_level
         )
-        all_none_vals.update(cls.model_construct().model_dump())
-        return all_none_vals
 
     def _to_yaml_obj(self, by_alias: bool = True) -> CommentedMap:
         # Make the YAML object to add comments to

src/dolphin/workflows/displacement.py

@@ -265,6 +265,7 @@ def run(
             num_threads=ts_opts.num_parallel_blocks,
             # TODO: do i care to configure block shape, or num threads from somewhere?
             # num_threads=cfg.worker_settings....?
+            wavelength=cfg.input_options.wavelength,
         )
 
     else:

tests/test_io_writers.py

@@ -6,7 +6,7 @@
 import rasterio as rio
 from rasterio.errors import NotGeoreferencedWarning
 
-from dolphin.io._core import load_gdal
+from dolphin.io._core import get_raster_units, load_gdal
 from dolphin.io._writers import (
     BackgroundBlockWriter,
     BackgroundRasterWriter,
@@ -150,3 +150,11 @@ def test_setitem(self, output_file_list, slc_file_list):
             assert np.allclose(
                 load_gdal(output_file_list[i], rows=rows, cols=cols), data[i]
             )
+
+    def test_file_kwargs(self, output_file_list, slc_file_list):
+        w = BackgroundStackWriter(
+            output_file_list, like_filename=slc_file_list[0], units="custom units"
+        )
+        w.close()
+        for f in output_file_list:
+            assert get_raster_units(f) == "custom units"