Add logic to reset the extra_reference_date after a network unwrapp…

…ing (#442)

* Add re-referencing for network unwrappers using `extra_reference_date`

* move out to separate helper function

* sort output paths

* delete

* add failing test

* for spurt, override the `interferogram_network` configuration

* add `IfgStackReader` for pre-formed nearest-3 interferograms

* add spurt to install

* fix optional type check

* revert spurt IfgReader

* final workaround to get spurt working

* work around for both spurt and for snaphu

* try adjusting deps

* fix disp reference test

.github/workflows/test-build-push.yml

@@ -46,12 +46,13 @@ jobs:
               - conda-forge
       - name: Install
         run: |
-          pip install --no-deps \
+          python -m pip install --no-deps \
             "snaphu>=0.4.0" \
             "opera-utils>=0.4.1" \
             git+https://github.com/isce-framework/tophu@main \
             git+https://github.com/isce-framework/whirlwind@40defb38d2d6deca2819934788ebbc57e418e32d
-          pip install --no-deps .
+          python -m pip install git+https://github.com/isce-framework/spurt@main
+          python -m pip install --no-deps .
       - name: Install test dependencies
         run: |
           micromamba install -f tests/requirements.txt -c conda-forge

src/dolphin/timeseries.py

@@ -1,6 +1,7 @@
 from __future__ import annotations
 
 import logging
+from datetime import datetime
 from enum import Enum
 from pathlib import Path
 from tempfile import NamedTemporaryFile
@@ -17,7 +18,7 @@
 from dolphin import DateOrDatetime, io, utils
 from dolphin._overviews import ImageType, create_overviews
 from dolphin._types import PathOrStr, ReferencePoint
-from dolphin.utils import flatten, format_dates, full_suffix
+from dolphin.utils import flatten, format_dates, full_suffix, get_nearest_date_idx
 from dolphin.workflows import CallFunc
 
 T = TypeVar("T")
@@ -55,6 +56,7 @@ def run(
     reference_point: tuple[int, int] = (-1, -1),
     wavelength: float | None = None,
     add_overviews: bool = True,
+    extra_reference_date: datetime | None = None,
 ) -> tuple[list[Path], ReferencePoint]:
     """Invert the unwrapped interferograms, estimate timeseries and phase velocity.
 
@@ -108,6 +110,9 @@ def run(
     add_overviews : bool, optional
         If True, creates overviews of the new velocity raster.
         Default is True.
+    extra_reference_date : datetime.datetime, optional
+        If provided, makes another set of interferograms referenced to this
+        for all dates later than it.
 
     Returns
     -------
@@ -203,7 +208,64 @@ def run(
             num_threads=num_threads,
         )
 
-    return inverted_phase_paths, ref_point
+    if extra_reference_date is None:
+        final_out = inverted_phase_paths
+    else:
+        final_out = _redo_reference(inverted_phase_paths, extra_reference_date)
+
+    return final_out, ref_point
+
+
+def _redo_reference(
+    inverted_phase_paths: Sequence[Path], extra_reference_date: datetime
+):
+    """Reset the reference date in `inverted_phase_paths`.
+
+    Affects all files whose secondary is after `extra_reference_date`.
+
+    E.g Given the (day 1, day 2), ..., (day 1, day N) pairs, outputs
+    (1, 2), (1, 3), ...(1, r), (r, r+1), ..., (r, N)
+    where r is the index of the `extra_reference_date`
+    """
+    output_path = inverted_phase_paths[0].parent
+    inverted_date_pairs: list[tuple[datetime, datetime]] = [
+        get_dates(p.stem)[:2] for p in inverted_phase_paths
+    ]
+    secondary_dates = [pair[1] for pair in inverted_date_pairs]
+    extra_ref_idx = get_nearest_date_idx(
+        secondary_dates, requested=extra_reference_date
+    )
+    ref_date = secondary_dates[extra_ref_idx]
+    logger.info(f"Re-referencing later timeseries files to {ref_date}")
+    extra_ref_img = inverted_phase_paths[extra_ref_idx]
+    ref = io.load_gdal(extra_ref_img, masked=True)
+
+    # Use a temp directory while re-referencing
+    extra_out_dir = inverted_phase_paths[0].parent / "extra"
+    extra_out_dir.mkdir(exist_ok=True)
+
+    for idx in range(extra_ref_idx + 1, len(inverted_date_pairs)):
+        # To create the interferogram (r, r+1), we subtract
+        # (1, r) from (1, r+1)
+        cur_img = inverted_phase_paths[idx]
+        new_stem = format_dates(ref_date, secondary_dates[idx])
+        cur_output_name = extra_out_dir / f"{new_stem}.tif"
+        cur = io.load_gdal(cur_img, masked=True)
+        new_out = cur - ref
+        io.write_arr(
+            arr=new_out, like_filename=extra_ref_img, output_name=cur_output_name
+        )
+
+    for idx, p in enumerate(inverted_phase_paths):
+        if idx <= extra_ref_idx:
+            p.rename(extra_out_dir / p.name)
+        else:
+            p.unlink()
+    # Finally, move them back in to the `timeseries/` folder
+    final_out = []
+    for p in extra_out_dir.glob("*.tif"):
+        final_out.append(p.rename(output_path / p.name))
+    return sorted(final_out)
 
 
 def _convert_and_reference(

src/dolphin/utils.py

@@ -690,3 +690,22 @@ def shutdown(self, wait: bool = True, cancel_futures: bool = True):  # noqa:D102
 
     def map(self, fn: Callable[P, T], *iterables, **kwargs):  # noqa: D102
         return map(fn, *iterables)
+
+
+def get_nearest_date_idx(
+    input_items: Sequence[datetime.datetime],
+    requested: datetime.datetime,
+) -> int:
+    """Find the index nearest to `requested` within `input_items`."""
+    sorted_inputs = sorted(input_items)
+    if not sorted_inputs[0] <= requested <= sorted_inputs[-1]:
+        msg = f"Requested {requested} falls outside of input range: "
+        msg += f"{sorted_inputs[0]}, {sorted_inputs[-1]}"
+        raise ValueError(msg)
+
+    nearest_idx = min(
+        range(len(input_items)),
+        key=lambda i: abs((input_items[i] - requested).total_seconds()),
+    )
+
+    return nearest_idx

src/dolphin/workflows/config/_displacement.py

@@ -27,7 +27,7 @@
     WorkflowBase,
     _read_file_list_or_glob,
 )
-from ._unwrap_options import UnwrapOptions
+from ._unwrap_options import UnwrapMethod, UnwrapOptions
 
 __all__ = [
     "DisplacementWorkflow",
@@ -246,3 +246,15 @@ def model_post_init(self, __context: Any) -> None:
         self.timeseries_options._velocity_file = (
             work_dir / self.timeseries_options._velocity_file
         )
+
+        # Modify interferogram options if using spurt for 3d unwrapping,
+        # which only does nearest-3 interferograms
+        if self.unwrap_options.unwrap_method == UnwrapMethod.SPURT:
+            logger.info(
+                "Using spurt: will form single reference interferograms, later convert"
+                " to nearest-3"
+            )
+            self.interferogram_network.reference_idx = 0
+            # Force all other network options to None
+            for attr in ["max_bandwidth", "max_temporal_baseline", "indexes"]:
+                setattr(self.interferogram_network, attr, None)

src/dolphin/workflows/displacement.py

@@ -271,6 +271,7 @@ def run(
             # num_threads=cfg.worker_settings....?
             wavelength=cfg.input_options.wavelength,
             add_overviews=cfg.output_options.add_overviews,
+            extra_reference_date=cfg.output_options.extra_reference_date,
         )
 
     else:

src/dolphin/workflows/wrapped_phase.py

@@ -11,6 +11,8 @@
 from dolphin import Bbox, Filename, interferogram, masking, ps
 from dolphin._log import log_runtime, setup_logging
 from dolphin.io import VRTStack
+from dolphin.utils import get_nearest_date_idx
+from dolphin.workflows import UnwrapMethod
 
 from . import InterferogramNetwork, sequential
 from .config import DisplacementWorkflow
@@ -139,7 +141,7 @@ def run(
 
     extra_reference_date = cfg.output_options.extra_reference_date
     if extra_reference_date:
-        new_compressed_slc_reference_idx = _get_nearest_idx(
+        new_compressed_slc_reference_idx = get_nearest_date_idx(
             [dtup[0] for dtup in input_dates], extra_reference_date
         )
     else:
@@ -216,18 +218,31 @@ def run(
         )
 
     logger.info(f"Creating virtual interferograms from {len(phase_linked_slcs)} files")
-    # TODO: with manual indexes, this may be split into 2 and redone
     reference_date = [
         get_dates(f, fmt=cfg.input_options.cslc_date_fmt)[0] for f in input_file_list
     ][cfg.phase_linking.output_reference_idx]
 
+    # TODO: remove this bad back to get around spurt's required input
+    # Reading direct nearest-3 ifgs is not working due to some slicing problem
+    # so we need to just give it single reference ifgs, all referenced to the beginning
+    # of the stack
+    extra_reference_date = (
+        # For spurt / networks of unwrappings, we ignore this this "changeover" date
+        # It will get applied in the `timeseries/` step
+        cfg.output_options.extra_reference_date
+        if (
+            cfg.unwrap_options.unwrap_method != UnwrapMethod.SPURT
+            or cfg.interferogram_network.max_bandwidth is not None
+        )
+        else None
+    )
     ifg_file_list: list[Path] = []
     ifg_file_list = create_ifgs(
         interferogram_network=ifg_network,
         phase_linked_slcs=phase_linked_slcs,
         contained_compressed_slcs=any(is_compressed),
         reference_date=reference_date,
-        extra_reference_date=cfg.output_options.extra_reference_date,
+        extra_reference_date=extra_reference_date,
     )
     return (
         ifg_file_list,
@@ -290,6 +305,8 @@ def create_ifgs(
     ifg_file_list: list[Path] = []
 
     secondary_dates = [get_dates(f)[0] for f in phase_linked_slcs]
+    # TODO: if we manually set an ifg network (i.e. not rely on spurt),
+    # we may still want to just pass it right to `Network`
     if not contained_compressed_slcs and extra_reference_date is None:
         # When no compressed SLCs/extra reference were passed in to the config,
         # we can directly pass options to `Network` and get the ifg list
@@ -326,7 +343,9 @@ def create_ifgs(
             for f in phase_linked_slcs
         ]
     else:
-        manual_reference_idx = _get_nearest_idx(secondary_dates, extra_reference_date)
+        manual_reference_idx = get_nearest_date_idx(
+            secondary_dates, extra_reference_date
+        )
 
         single_ref_ifgs = [
             interferogram.convert_pl_to_ifg(
@@ -474,22 +493,3 @@ def _get_mask(
         mask_filename = nodata_mask_file
 
     return mask_filename
-
-
-def _get_nearest_idx(
-    input_dates: Sequence[datetime.datetime],
-    selected_date: datetime.datetime,
-) -> int:
-    """Find the index nearest to `selected_date` within `input_dates`."""
-    sorted_inputs = sorted(input_dates)
-    if not sorted_inputs[0] <= selected_date <= sorted_inputs[-1]:
-        msg = f"Requested {selected_date} falls outside of input range: "
-        msg += f"{sorted_inputs[0]}, {sorted_inputs[-1]}"
-        raise ValueError(msg)
-
-    nearest_idx = min(
-        range(len(input_dates)),
-        key=lambda i: abs((input_dates[i] - selected_date).total_seconds()),
-    )
-
-    return nearest_idx

tests/test_workflows_displacement.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import importlib.util
 from pathlib import Path
 
 import pytest
@@ -188,14 +189,20 @@ def test_separate_workflow_runs(slc_file_list, tmp_path):
     assert [f.name for f in ifgs3_b] == [f.name for f in ifgs3]
 
 
-def test_displacement_run_extra_reference_date(opera_slc_files: list[Path], tmpdir):
+@pytest.mark.parametrize("unwrap_method", ["phass", "spurt"])
+def test_displacement_run_extra_reference_date(
+    opera_slc_files: list[Path], tmpdir, unwrap_method: str
+):
+    if unwrap_method == "spurt" and importlib.util.find_spec("spurt") is None:
+        pytest.skip(reason="spurt unwrapper not installed")
+
     with tmpdir.as_cwd():
         cfg = config.DisplacementWorkflow(
             # start_date = 20220101
             # shape = (4, 128, 128)
             # First one is COMPRESSED_
             output_options={"extra_reference_date": "2022-01-03"},
-            unwrap_options={"unwrap_method": "phass"},
+            unwrap_options={"unwrap_method": unwrap_method},
             cslc_file_list=opera_slc_files,
             input_options={"subdataset": "/data/VV"},
             phase_linking={
@@ -210,16 +217,28 @@ def test_displacement_run_extra_reference_date(opera_slc_files: list[Path], tmpd
 
         # The unwrapped/timeseries files should have a changeover to the new reference
         assert paths.unwrapped_paths is not None
-        unw_names = [pp.name for pp in paths.unwrapped_paths]
-        assert unw_names == [
-            "20220101_20220102.unw.tif",
-            "20220101_20220103.unw.tif",
-            "20220103_20220104.unw.tif",
-        ]
         assert paths.timeseries_paths is not None
+
         ts_names = [pp.name for pp in paths.timeseries_paths]
         assert ts_names == [
             "20220101_20220102.tif",
             "20220101_20220103.tif",
             "20220103_20220104.tif",
         ]
+
+        unw_names = [pp.name for pp in paths.unwrapped_paths]
+        if cfg.unwrap_options.unwrap_method == "spurt":
+            assert unw_names == [
+                "20220101_20220102.unw.tif",
+                "20220101_20220103.unw.tif",
+                "20220101_20220104.unw.tif",
+                "20220102_20220103.unw.tif",
+                "20220102_20220104.unw.tif",
+                "20220103_20220104.unw.tif",
+            ]
+        else:
+            assert unw_names == [
+                "20220101_20220102.unw.tif",
+                "20220101_20220103.unw.tif",
+                "20220103_20220104.unw.tif",
+            ]