diff --git a/CHANGELOG.md b/CHANGELOG.md index 1bb5d0c8..c44dc97d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,8 +1,26 @@ # Changelog -## [Unreleased](https://github.com/isce-framework/dolphin/compare/v0.30.0...main) +## [Unreleased](https://github.com/isce-framework/dolphin/compare/v0.35.1...main) -## [0.35.0](https://github.com/isce-framework/dolphin/compare/v0.34.0...v0.35.0) - 2025-12-09 +### Changed + +- `timeseries.py`: Auto reference point selection + - Picks the center of mass instead of arbitrary `argmax` result + - Rename `condition_file` to `quality_file` + +### Changed +- Removed `condition` parameter in `timeeries` reference point functions +- Removed `CallFunc` enum + +## [0.35.1](https://github.com/isce-framework/dolphin/compare/v0.35.0...v0.35.1) - 2025-01-15 + +### Fixed + +- `filtering.py` Fix in_bounds_pixels masking, set default to 25 km +- Set `output_reference_idx` separately from `compressed_reference_idx` during phase linking setup + + +## [0.35.0](https://github.com/isce-framework/dolphin/compare/v0.34.0...v0.35.0) - 2025-01-09 ### Added diff --git a/src/dolphin/_cli_timeseries.py b/src/dolphin/_cli_timeseries.py index 567904e2..f220de72 100644 --- a/src/dolphin/_cli_timeseries.py +++ b/src/dolphin/_cli_timeseries.py @@ -5,7 +5,6 @@ from dolphin._log import setup_logging from dolphin.timeseries import InversionMethod -from dolphin.workflows import CallFunc if TYPE_CHECKING: _SubparserType = argparse._SubParsersAction[argparse.ArgumentParser] @@ -57,21 +56,12 @@ def get_parser(subparser=None, subcommand_name="timeseries") -> argparse.Argumen ), ) parser.add_argument( - "--condition-file", + "--quality-file", help=( "A file with the same size as each raster, like amplitude dispersion or " "temporal coherence to find reference point" ), ) - parser.add_argument( - "--condition", - type=CallFunc, - default=CallFunc.MIN, - help=( - "A condition to apply to condition file to find the reference point. " - "Options are [min, max]. default=min" - ), - ) parser.add_argument( "--method", type=InversionMethod, diff --git a/src/dolphin/timeseries.py b/src/dolphin/timeseries.py index 2a6820f3..e561d856 100644 --- a/src/dolphin/timeseries.py +++ b/src/dolphin/timeseries.py @@ -19,7 +19,6 @@ from dolphin._overviews import ImageType, create_overviews from dolphin._types import PathOrStr, ReferencePoint from dolphin.utils import flatten, format_dates, full_suffix, get_nearest_date_idx -from dolphin.workflows import CallFunc T = TypeVar("T") DateOrDatetime = datetime | date @@ -42,8 +41,7 @@ class ReferencePointError(ValueError): def run( unwrapped_paths: Sequence[PathOrStr], conncomp_paths: Sequence[PathOrStr] | None, - condition_file: PathOrStr, - condition: CallFunc, + quality_file: PathOrStr, output_dir: PathOrStr, method: InversionMethod = InversionMethod.L1, reference_candidate_threshold: float = 0.95, @@ -67,13 +65,9 @@ def run( Sequence unwrapped interferograms to invert. conncomp_paths : Sequence[Path] Sequence connected component files, one per file in `unwrapped_paths` - condition_file: PathOrStr + quality_file: PathOrStr A file with the same size as each raster, like amplitude dispersion or temporal coherence - condition: CallFunc - The function to apply to the condition file, - for example numpy.argmin which finds the pixel with lowest value - the options are [min, max] output_dir : Path Path to the output directory. method : str, choices = "L1", "L2" @@ -81,12 +75,11 @@ def run( Default is L2, which uses least squares to solve Ax = b (faster). "L1" minimizes |Ax - b|_1 at each pixel. reference_candidate_threshold: float - The threshold for the condition function to be considered a candidate + The threshold for the quality metric to be considered a candidate reference point pixel. - If the condition function is CallFunc.MAX, then only pixels with values - in `condition_file` greater than `reference_candidate_threshold` will be - considered a candidate. - Default = 0.95 + Only pixels with values in `quality_file` greater than + `reference_candidate_threshold` will be considered a candidate. + Default is 0.95. run_velocity : bool Whether to run velocity estimation on the inverted phase series corr_paths : Sequence[Path], optional @@ -147,9 +140,8 @@ def run( if reference_point == (-1, -1): logger.info("Selecting a reference point for unwrapped interferograms") ref_point = select_reference_point( - condition_file=condition_file, + quality_file=quality_file, output_dir=Path(output_dir), - condition=condition, candidate_threshold=reference_candidate_threshold, ccl_file_list=conncomp_paths, ) @@ -1176,9 +1168,8 @@ def correlation_to_variance(correlation: ArrayLike, nlooks: int) -> Array: def select_reference_point( *, - condition_file: PathOrStr, + quality_file: PathOrStr, output_dir: Path, - condition: CallFunc = CallFunc.MAX, candidate_threshold: float = 0.95, ccl_file_list: Sequence[PathOrStr] | None = None, block_shape: tuple[int, int] = (256, 256), @@ -1186,30 +1177,27 @@ def select_reference_point( ) -> ReferencePoint: """Automatically select a reference point for a stack of unwrapped interferograms. - Uses the condition file and (optionally) connected component labels. + Uses the quality file and (optionally) connected component labels. The point is selected which - 1. has the condition applied to condition file. for example: has the lowest - amplitude dispersion - 2. (optionally) is within intersection of all nonzero connected component labels + 1. (optionally) is within intersection of all nonzero connected component labels + 2. Has value in `quality_file` above the threshold `candidate_threshold` + + Among all points which meet this, the centroid selected using the function + `scipy.ndimage.center_of_mass`. Parameters ---------- - condition_file: PathOrStr - A file with the same size as each raster, like amplitude dispersion or - temporal coherence in `ccl_file_list` + quality_file: PathOrStr + A file with the same size as each raster in `ccl_file_list` containing a quality + metric, such as temporal coherence. output_dir: Path Path to store the computed "conncomp_intersection.tif" raster - condition: CallFunc - The function to apply to the condition file, for example CallFunc.MIN - which finds the pixel with lowest value. - Default = CallFunc.MAX candidate_threshold: float - The threshold for the condition function to be considered a candidate + The threshold for the quality metric function to be considered a candidate reference point pixel. - If the condition function is CallFunc.MAX, then only pixels with values - in `condition_file` greater than `candidate_threshold` will be considered - a candidate. + Only pixels with values in `quality_file` greater than `candidate_threshold` are + considered a candidate. Default = 0.95 ccl_file_list : Sequence[PathOrStr] List of connected component label phase files. @@ -1239,10 +1227,10 @@ def select_reference_point( return ref_point logger.info("Selecting reference point") - condition_file_values = io.load_gdal(condition_file, masked=True) + quality_file_values = io.load_gdal(quality_file, masked=True) # Start with all points as valid candidates - isin_largest_conncomp = np.ones(condition_file_values.shape, dtype=bool) + isin_largest_conncomp = np.ones(quality_file_values.shape, dtype=bool) if ccl_file_list: try: isin_largest_conncomp = _get_largest_conncomp_mask( @@ -1253,14 +1241,11 @@ def select_reference_point( ) except ReferencePointError: msg = "Unable to find find a connected component intersection." - msg += f"Proceeding using only {condition_file = }" + msg += f"Proceeding using only {quality_file = }" logger.warning(msg, exc_info=True) # Find pixels meeting the threshold criteria - if condition == CallFunc.MAX: - is_candidate = condition_file_values > candidate_threshold - else: # CallFunc.MIN - is_candidate = condition_file_values < candidate_threshold + is_candidate = quality_file_values > candidate_threshold # Restrict candidates to the largest connected component region is_candidate &= isin_largest_conncomp @@ -1274,8 +1259,7 @@ def select_reference_point( f"No pixels above threshold={candidate_threshold}. Choosing best among" " available." ) - condition_func = argmax_index if condition == CallFunc.MAX else argmin_index - ref_row, ref_col = condition_func(condition_file_values) + ref_row, ref_col = argmax_index(quality_file_values) else: # Find the largest region of connected candidate pixels label_counts = np.bincount(labeled.ravel()) diff --git a/src/dolphin/workflows/config/_enums.py b/src/dolphin/workflows/config/_enums.py index 40500287..efb2514d 100644 --- a/src/dolphin/workflows/config/_enums.py +++ b/src/dolphin/workflows/config/_enums.py @@ -1,7 +1,6 @@ from enum import Enum __all__ = [ - "CallFunc", "ShpMethod", "UnwrapMethod", ] @@ -25,10 +24,3 @@ class UnwrapMethod(str, Enum): PHASS = "phass" SPURT = "spurt" WHIRLWIND = "whirlwind" - - -class CallFunc(str, Enum): - """Call function for the timeseries method to find reference point.""" - - MIN = "min" - MAX = "max" diff --git a/src/dolphin/workflows/displacement.py b/src/dolphin/workflows/displacement.py index 14fd5648..097ff753 100755 --- a/src/dolphin/workflows/displacement.py +++ b/src/dolphin/workflows/displacement.py @@ -16,7 +16,6 @@ from dolphin import __version__, io, timeseries, utils from dolphin._log import log_runtime, setup_logging from dolphin.timeseries import ReferencePoint -from dolphin.workflows import CallFunc from . import stitching_bursts, unwrapping, wrapped_phase from ._utils import _create_burst_cfg, _remove_dir_if_empty, parse_ionosphere_files @@ -260,8 +259,10 @@ def run( unwrapped_paths=unwrapped_paths, conncomp_paths=conncomp_paths, corr_paths=stitched_paths.interferometric_corr_paths, - condition_file=stitched_paths.temp_coh_file, - condition=CallFunc.MAX, + quality_file=stitched_paths.temp_coh_file, + # TODO: Right now aew don't have the option to pick a different candidate + # or quality file. + reference_candidate_threshold=0.95, output_dir=ts_opts._directory, method=timeseries.InversionMethod(ts_opts.method), run_velocity=ts_opts.run_velocity, diff --git a/tests/test_timeseries.py b/tests/test_timeseries.py index 5766db07..4c503ab6 100644 --- a/tests/test_timeseries.py +++ b/tests/test_timeseries.py @@ -7,7 +7,6 @@ from numpy.linalg import lstsq as lstsq_numpy from dolphin import io, timeseries -from dolphin.timeseries import CallFunc from dolphin.utils import format_dates NUM_DATES = 10 @@ -283,9 +282,8 @@ def test_all_ones_center(self, tmp_path): io.write_arr(arr=arr, output_name=coh_file) ref_point = timeseries.select_reference_point( - condition_file=coh_file, + quality_file=coh_file, output_dir=tmp_path, - condition=CallFunc.MAX, candidate_threshold=0.95, # everything is above 0.95 ccl_file_list=None, ) @@ -303,9 +301,8 @@ def test_half_03_half_099(self, tmp_path): io.write_arr(arr=arr, output_name=coh_file) ref_point = timeseries.select_reference_point( - condition_file=coh_file, + quality_file=coh_file, output_dir=tmp_path, - condition=CallFunc.MAX, candidate_threshold=0.95, ccl_file_list=None, ) @@ -342,9 +339,8 @@ def test_with_conncomp(self, tmp_path): ref_point = timeseries.select_reference_point( ccl_file_list=[ccl_file1, ccl_file2], - condition_file=coh_file, + quality_file=coh_file, output_dir=tmp_path, - condition=CallFunc.MAX, candidate_threshold=0.95, )