From e776da6db4994c2d1336e343e56cfdd5aba745b6 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sat, 12 Oct 2024 00:01:14 -0400 Subject: [PATCH] Add logic to reset the `extra_reference_date` after a network unwrapping (#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 | 5 +- src/dolphin/phase_link/_core.jl | 19 ------ src/dolphin/timeseries.py | 66 ++++++++++++++++++- src/dolphin/utils.py | 19 ++++++ src/dolphin/workflows/config/_displacement.py | 14 +++- src/dolphin/workflows/displacement.py | 1 + src/dolphin/workflows/wrapped_phase.py | 46 ++++++------- tests/test_workflows_displacement.py | 35 +++++++--- 8 files changed, 150 insertions(+), 55 deletions(-) delete mode 100644 src/dolphin/phase_link/_core.jl diff --git a/.github/workflows/test-build-push.yml b/.github/workflows/test-build-push.yml index cd1804bb..9f98296b 100644 --- a/.github/workflows/test-build-push.yml +++ b/.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 diff --git a/src/dolphin/phase_link/_core.jl b/src/dolphin/phase_link/_core.jl deleted file mode 100644 index 619672f9..00000000 --- a/src/dolphin/phase_link/_core.jl +++ /dev/null @@ -1,19 +0,0 @@ -function emi(C) - Gam = abs.(C) - Gam_inv = inv(Gam) - A = Gam_inv .* C - - - σ = 0.99 - F = lu(A - σ * I) - Fmap = LinearMap{ComplexF64}((y, x) -> ldiv!(y, F, x), 25, ismutating=true) - - x = ones(eltype(A), size(A, 1)) - # lam_, v = IterativeSolvers.powm!(A, x; inverse=true, shift=σ) - # lam_, v = IterativeSolvers.powm!(Fmap, x; inverse=true, shift=σ) - # lam_, v = IterativeSolvers.powm!(Fmap, x; inverse=true, shift=σ) - # λ, x = invpowm!(Fmap, x) - # λ, x = invpowm!(Fmap, x) - λ, x = invpowm(Fmap, shift=σ) - return λ, v -end diff --git a/src/dolphin/timeseries.py b/src/dolphin/timeseries.py index 32469632..1266355c 100644 --- a/src/dolphin/timeseries.py +++ b/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( diff --git a/src/dolphin/utils.py b/src/dolphin/utils.py index a5cb9789..4505aba1 100644 --- a/src/dolphin/utils.py +++ b/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 diff --git a/src/dolphin/workflows/config/_displacement.py b/src/dolphin/workflows/config/_displacement.py index 067e1fe3..9c86330f 100644 --- a/src/dolphin/workflows/config/_displacement.py +++ b/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) diff --git a/src/dolphin/workflows/displacement.py b/src/dolphin/workflows/displacement.py index 998fa47b..936f3b28 100755 --- a/src/dolphin/workflows/displacement.py +++ b/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: diff --git a/src/dolphin/workflows/wrapped_phase.py b/src/dolphin/workflows/wrapped_phase.py index 5c992fd9..7ff71d0d 100644 --- a/src/dolphin/workflows/wrapped_phase.py +++ b/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 diff --git a/tests/test_workflows_displacement.py b/tests/test_workflows_displacement.py index fbd062fc..034e419e 100644 --- a/tests/test_workflows_displacement.py +++ b/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", + ]