diff --git a/src/dolphin/phase_link/_compress.py b/src/dolphin/phase_link/_compress.py index e020ca80..2df1e38f 100644 --- a/src/dolphin/phase_link/_compress.py +++ b/src/dolphin/phase_link/_compress.py @@ -9,6 +9,7 @@ def compress( slc_stack: ArrayLike, pl_cpx_phase: ArrayLike, + first_real_slc_idx: int = 0, slc_mean: ArrayLike | None = None, reference_idx: int | None = None, ): @@ -21,6 +22,10 @@ def compress( pl_cpx_phase : ArrayLike The estimated complex phase from phase linking. shape = (nslc, rows // strides.y, cols // strides.x) + first_real_slc_idx : int + Index within `slc_stack` where the "real" SLCs start. + Indexes before this (which are assumed to be already Compressed SLCs) are + excluded from the dot product during compression slc_mean : ArrayLike, optional The mean SLC magnitude, shape (rows, cols), to use as output pixel magnitudes. If None, the mean is computed from the input SLC stack. @@ -39,18 +44,20 @@ def compress( pl_referenced = pl_cpx_phase * pl_cpx_phase[reference_idx][None, :, :].conj() else: pl_referenced = pl_cpx_phase + + # Slice away the compressed layers *after* doing the reference. + pl_referenced = pl_referenced[first_real_slc_idx:, :, :] + slcs = slc_stack[first_real_slc_idx:] # If the output is downsampled, we need to make `pl_cpx_phase` the same shape # as the output - pl_estimate_upsampled = upsample_nearest(pl_referenced, slc_stack.shape[1:]) + pl_estimate_upsampled = upsample_nearest(pl_referenced, slcs.shape[1:]) # For each pixel, project the SLCs onto the (normalized) estimated phase # by performing a pixel-wise complex dot product with warnings.catch_warnings(): warnings.filterwarnings("ignore", message="invalid value encountered") - phase = np.angle( - np.nansum(slc_stack * np.conjugate(pl_estimate_upsampled), axis=0) - ) + phase = np.angle(np.nansum(slcs * np.conjugate(pl_estimate_upsampled), axis=0)) if slc_mean is None: - slc_mean = np.mean(np.abs(slc_stack), axis=0) + slc_mean = np.mean(np.abs(slcs), axis=0) # If the phase is invalid, set the mean to NaN slc_mean[phase == 0] = np.nan return slc_mean * np.exp(1j * phase) diff --git a/src/dolphin/stack.py b/src/dolphin/stack.py index 273a7c91..c192550e 100755 --- a/src/dolphin/stack.py +++ b/src/dolphin/stack.py @@ -5,6 +5,7 @@ import logging import warnings from datetime import date, datetime +from enum import Enum from os import fspath from pathlib import Path from typing import Optional, Sequence @@ -20,6 +21,14 @@ logger = logging.getLogger(__name__) +class CompressedSlcPlan(str, Enum): + """Plan for creating Compressed SLCs during phase linking.""" + + ALWAYS_FIRST = "always_first" + FIRST_PER_MINISTACK = "first_per_ministack" + LAST_PER_MINISTACK = "last_per_ministack" + + class BaseStack(BaseModel): """Base class for mini- and full stack classes.""" @@ -342,6 +351,11 @@ class MiniStackInfo(BaseStack): ), ) + def __rich_repr__(self): + yield from super().__rich_repr__() + yield "output_reference_idx", self.output_reference_idx + yield "compressed_reference_idx", self.compressed_reference_idx + @property def output_reference_date(self): """Date of the reference phase of the stack.""" @@ -391,6 +405,7 @@ class MiniStackPlanner(BaseStack): 0, description="Index of the SLC to use as reference during phase linking", ) + compressed_slc_plan: CompressedSlcPlan = CompressedSlcPlan.ALWAYS_FIRST def plan( self, ministack_size: int, compressed_idx: int | None = None @@ -448,12 +463,20 @@ def plan( ) cur_output_folder = self.output_folder / new_date_str - if compressed_idx is None: - # TODO: To change to Ansari-style ministack references, change this - # so that a new `output_reference_idx` gets made - compressed_reference_idx = self.output_reference_idx - else: + if compressed_idx is not None: compressed_reference_idx = compressed_idx + elif self.compressed_slc_plan == CompressedSlcPlan.ALWAYS_FIRST: + # Simplest operational version: CompSLCs have same base phase, + # but different "residual" added on + compressed_reference_idx = 0 + elif self.compressed_slc_plan == CompressedSlcPlan.FIRST_PER_MINISTACK: + # Like Ansari, 2017 paper: each ministack is "self contained" + compressed_reference_idx = num_ccslc + # Ansari, 2017 also had output_reference_idx = num_ccslcs, and + # use the "Datum Adjustment" step to get outputs relative to day 0 + elif self.compressed_slc_plan == CompressedSlcPlan.LAST_PER_MINISTACK: + # Alternative that allows sequential interferograms across ministacks + compressed_reference_idx = -1 cur_ministack = MiniStackInfo( file_list=combined_files, @@ -469,3 +492,9 @@ def plan( compressed_slc_infos.append(cur_comp_slc) return output_ministacks + + def __rich_repr__(self): + yield from super().__rich_repr__() + yield "max_num_compressed", self.max_num_compressed + yield "output_reference_idx", self.output_reference_idx + yield "compressed_slc_plan", self.compressed_slc_plan diff --git a/src/dolphin/workflows/config/_common.py b/src/dolphin/workflows/config/_common.py index 21d1bbb8..f83ff0c5 100644 --- a/src/dolphin/workflows/config/_common.py +++ b/src/dolphin/workflows/config/_common.py @@ -19,6 +19,7 @@ from dolphin import __version__ as _dolphin_version from dolphin._types import Bbox from dolphin.io import DEFAULT_HDF5_OPTIONS, DEFAULT_TIFF_OPTIONS +from dolphin.stack import CompressedSlcPlan from ._enums import ShpMethod from ._yaml_model import YamlModel @@ -129,6 +130,7 @@ class PhaseLinkingOptions(BaseModel, extra="forbid"): "`n` interferograms. `baseline_line` must be positive." ), ) + compressed_slc_plan: CompressedSlcPlan = CompressedSlcPlan.ALWAYS_FIRST class InterferogramNetwork(BaseModel, extra="forbid"): diff --git a/src/dolphin/workflows/sequential.py b/src/dolphin/workflows/sequential.py index 26fc2354..4ffd9402 100644 --- a/src/dolphin/workflows/sequential.py +++ b/src/dolphin/workflows/sequential.py @@ -19,7 +19,7 @@ from dolphin._types import Filename from dolphin.io import VRTStack from dolphin.similarity import create_similarities -from dolphin.stack import MiniStackPlanner +from dolphin.stack import CompressedSlcPlan, MiniStackPlanner from .config import ShpMethod from .single import run_wrapped_phase_single @@ -45,6 +45,7 @@ def run_wrapped_phase_sequential( shp_nslc: Optional[int] = None, use_evd: bool = False, beta: float = 0.00, + compressed_slc_plan: CompressedSlcPlan = CompressedSlcPlan.ALWAYS_FIRST, max_num_compressed: int = 100, output_reference_idx: int = 0, new_compressed_reference_idx: int | None = None, @@ -69,6 +70,7 @@ def run_wrapped_phase_sequential( output_folder=output_folder, max_num_compressed=max_num_compressed, output_reference_idx=output_reference_idx, + compressed_slc_plan=compressed_slc_plan, ) ministacks = ministack_planner.plan( ministack_size, compressed_idx=new_compressed_reference_idx diff --git a/src/dolphin/workflows/single.py b/src/dolphin/workflows/single.py index fde676a0..23adff15 100644 --- a/src/dolphin/workflows/single.py +++ b/src/dolphin/workflows/single.py @@ -237,18 +237,14 @@ def run_wrapped_phase_single( # Get the mean to set as pixel magnitudes abs_stack = np.abs(cur_data[first_real_slc_idx:, in_trim_rows, in_trim_cols]) cur_data_mean, cur_amp_dispersion, _ = calc_ps_block(abs_stack) - # NOTE: the `ministack.compressed_reference_idx` was set relative to *all* - # input images, real and compressed. - # Since we pass in just the first real idx, we have to subtract that off - ref_index_relative = ministack.compressed_reference_idx - first_real_slc_idx cur_comp_slc = compress( # Get the inner portion of the full-res SLC data - cur_data[first_real_slc_idx:, in_trim_rows, in_trim_cols], - pl_output.cpx_phase[first_real_slc_idx:, out_trim_rows, out_trim_cols], + cur_data[:, in_trim_rows, in_trim_cols], + pl_output.cpx_phase[:, out_trim_rows, out_trim_cols], + first_real_slc_idx=first_real_slc_idx, slc_mean=cur_data_mean, - reference_idx=ref_index_relative, + reference_idx=ministack.compressed_reference_idx, ) - # TODO: truncate # ### Save results ### diff --git a/tests/test_show_versions.py b/tests/test_show_versions.py index a2ba7c38..a61cf5ba 100644 --- a/tests/test_show_versions.py +++ b/tests/test_show_versions.py @@ -15,7 +15,6 @@ def test_get_deps_info(): assert "osgeo.gdal" in deps_info assert "numpy" in deps_info assert "pydantic" in deps_info - assert "ruamel_yaml" in deps_info assert "h5py" in deps_info