diff --git a/src/dolphin/_cli_unwrap.py b/src/dolphin/_cli_unwrap.py index ed0786d60..12cdc9e6a 100644 --- a/src/dolphin/_cli_unwrap.py +++ b/src/dolphin/_cli_unwrap.py @@ -77,6 +77,16 @@ def get_parser(subparser=None, subcommand_name="unwrap") -> argparse.ArgumentPar default=UnwrapMethod.SNAPHU.value, help="Choice of unwrapping algorithm to use.", ) + algorithm_opts.add_argument( + "--run-goldstein", + action="store_true", + help="Run Goldstein filter before unwrapping.", + ) + algorithm_opts.add_argument( + "--run-interpolation", + action="store_true", + help="Run interpolation before unwrapping.", + ) tophu_opts = parser.add_argument_group("Tophu options") # Add ability for downsampling/tiling with tophu diff --git a/src/dolphin/timeseries.py b/src/dolphin/timeseries.py index ff08720ad..b931ac010 100644 --- a/src/dolphin/timeseries.py +++ b/src/dolphin/timeseries.py @@ -35,6 +35,26 @@ class ReferencePointError(ValueError): pass +def argmin_index(arr: ArrayLike) -> tuple[int, ...]: + """Get the index tuple of the minimum value of the array. + + If multiple occurrences of the minimum value exist, returns + the index of the first such occurrence in the flattened array. + + Parameters + ---------- + arr : array_like + The input array. + + Returns + ------- + tuple of int + The index of the minimum value. + + """ + return np.unravel_index(np.argmin(arr), np.shape(arr)) + + @jit def weighted_lstsq_single( A: ArrayLike, @@ -607,28 +627,33 @@ def correlation_to_variance(correlation: ArrayLike, nlooks: int) -> Array: def select_reference_point( ccl_file_list: Sequence[PathOrStr], - amp_dispersion_file: PathOrStr, + condition_file: PathOrStr, output_dir: Path, + condition_func: Callable[[ArrayLike], tuple[int, ...]] = argmin_index, block_shape: tuple[int, int] = (512, 512), num_threads: int = 5, ) -> ReferencePoint: """Automatically select a reference point for a stack of unwrapped interferograms. - Uses the amplitude dispersion and connected component labels, the point is selected + Uses the condition file and connected component labels, the point is selected which 1. is within intersection of all nonzero connected component labels (always valid) - 2. has the lowest amplitude dispersion + 2. has the condition applied to condition file. for example: has the lowest + amplitude dispersion Parameters ---------- ccl_file_list : Sequence[PathOrStr] List of connected component label phase files. - amp_dispersion_file: PathOrStr - Amplitude dispersion, stitched and multilooked to the same size each raster - in `ccl_file_list` + condition_file: PathOrStr + A file with the same size as each raster, like amplitude dispersion or + temporal coherence in `ccl_file_list` output_dir: Path Path to store the computed "conncomp_intersection.tif" raster + condition_func: Callable[[ArrayLike, ]] + The function to apply to the condition file, + for example numpy.argmin which finds the pixel with lowest value block_shape: tuple[int, int] Size of blocks to read from while processing `ccl_file_list` Default = (512, 512) @@ -693,12 +718,12 @@ def intersect_conncomp(arr: np.ma.MaskedArray, axis: int) -> np.ndarray: # Create a mask of pixels with this label isin_largest_conncomp = label == largest_idx - amp_dispersion = io.load_gdal(amp_dispersion_file, masked=True) + condition_file_values = io.load_gdal(condition_file, masked=True) # Mask out where the conncomps aren't equal to the largest - amp_dispersion.mask = amp_dispersion.mask | (~isin_largest_conncomp) + condition_file_values.mask = condition_file_values.mask | (~isin_largest_conncomp) - # Pick the (unmasked) point with the lowest amplitude dispersion - ref_row, ref_col = np.unravel_index(np.argmin(amp_dispersion), amp_dispersion.shape) + # Pick the (unmasked) point with the condition applied to condition file + ref_row, ref_col = condition_func(condition_file_values) # Cast to `int` to avoid having `np.int64` types return ReferencePoint(int(ref_row), int(ref_col)) diff --git a/src/dolphin/unwrap/_snaphu_py.py b/src/dolphin/unwrap/_snaphu_py.py index ee223e775..3e6fe2c88 100644 --- a/src/dolphin/unwrap/_snaphu_py.py +++ b/src/dolphin/unwrap/_snaphu_py.py @@ -1,6 +1,9 @@ from __future__ import annotations from pathlib import Path +from typing import Optional + +from numpy.typing import ArrayLike from dolphin._log import get_log from dolphin._types import Filename @@ -21,13 +24,13 @@ def unwrap_snaphu_py( ntiles: tuple[int, int] = (1, 1), tile_overlap: tuple[int, int] = (0, 0), nproc: int = 1, - mask_file: Filename | None = None, + mask_file: Optional[Filename] = None, zero_where_masked: bool = False, - unw_nodata: float | None = DEFAULT_UNW_NODATA, - ccl_nodata: int | None = DEFAULT_CCL_NODATA, + unw_nodata: Optional[float] = DEFAULT_UNW_NODATA, + ccl_nodata: Optional[int] = DEFAULT_CCL_NODATA, init_method: str = "mst", cost: str = "smooth", - scratchdir: Filename | None = None, + scratchdir: Optional[Filename] = None, ) -> tuple[Path, Path]: """Unwrap an interferogram using at multiple scales using `tophu`. @@ -151,3 +154,74 @@ def unwrap_snaphu_py( return _zero_from_mask(unw_filename, cc_filename, mask_file) return Path(unw_filename), Path(cc_filename) + + +def grow_conncomp_snaphu( + unw_filename: Filename, + corr_filename: Filename, + nlooks: float, + mask: Optional[ArrayLike] = None, + ccl_nodata: Optional[int] = DEFAULT_CCL_NODATA, + cost: str = "smooth", + min_conncomp_frac: float = 0.0001, + scratchdir: Optional[Filename] = None, +) -> Path: + """Compute connected component labels using SNAPHU. + + Parameters + ---------- + unw_filename : Filename + Path to output unwrapped phase file. + corr_filename : Filename + Path to input correlation file. + nlooks : float + Effective number of looks used to form the input correlation data. + mask : Array, optional + binary byte mask array, by default None. + Assumes that 1s are valid pixels and 0s are invalid. + ccl_nodata : float, optional + Nodata value for the connected component labels. + cost : str + Statistical cost mode. + Default = "smooth" + min_conncomp_frac : float, optional + Minimum size of a single connected component, as a fraction of the total number + of pixels in the array. Defaults to 0.0001. + scratchdir : Filename, optional + If provided, uses a scratch directory to save the intermediate files + during unwrapping. + + Returns + ------- + conncomp_path : Filename + Path to output connected component label file. + + """ + import snaphu + + unw_suffix = full_suffix(unw_filename) + cc_filename = str(unw_filename).replace(unw_suffix, CONNCOMP_SUFFIX) + + with ( + snaphu.io.Raster(unw_filename) as unw, + snaphu.io.Raster(corr_filename) as corr, + snaphu.io.Raster.create( + cc_filename, + like=unw, + nodata=ccl_nodata, + dtype="u2", + **DEFAULT_TIFF_OPTIONS_RIO, + ) as conncomp, + ): + snaphu.grow_conncomps( + unw=unw, + corr=corr, + nlooks=nlooks, + mask=mask, + cost=cost, + min_conncomp_frac=min_conncomp_frac, + scratchdir=scratchdir, + conncomp=conncomp, + ) + + return Path(cc_filename) diff --git a/src/dolphin/unwrap/_unwrap.py b/src/dolphin/unwrap/_unwrap.py index ec3d3f43c..bbb459234 100644 --- a/src/dolphin/unwrap/_unwrap.py +++ b/src/dolphin/unwrap/_unwrap.py @@ -21,7 +21,7 @@ UNW_SUFFIX, UNW_SUFFIX_ZEROED, ) -from ._snaphu_py import unwrap_snaphu_py +from ._snaphu_py import grow_conncomp_snaphu, unwrap_snaphu_py from ._tophu import multiscale_unwrap from ._utils import create_combined_mask, set_nodata_values @@ -37,7 +37,7 @@ def run( output_path: Filename, *, nlooks: float = 5, - mask_file: Optional[Filename] = None, + mask_filename: Optional[Filename] = None, zero_where_masked: bool = False, unwrap_method: UnwrapMethod = UnwrapMethod.SNAPHU, init_method: str = "mst", @@ -69,7 +69,7 @@ def run( Path to output directory. nlooks : int, optional Effective number of looks used to form the input correlation data. - mask_file : Filename, optional + mask_filename : Filename, optional Path to binary byte mask file, by default None. Assumes that 1s are valid pixels and 0s are invalid. zero_where_masked : bool, optional @@ -160,8 +160,8 @@ def run( out_files.append(outf) logger.info(f"{len(out_files)} left to unwrap") - if mask_file: - mask_file = Path(mask_file).resolve() + if mask_filename: + mask_filename = Path(mask_filename).resolve() # This keeps it from spawning a new process for a single job. Executor = ThreadPoolExecutor if max_jobs > 1 else DummyProcessPoolExecutor @@ -176,7 +176,7 @@ def run( init_method=init_method, cost=cost, unwrap_method=unwrap_method, - mask_file=mask_file, + mask_filename=mask_filename, zero_where_masked=zero_where_masked, downsample_factor=downsample_factor, ntiles=ntiles, @@ -197,7 +197,7 @@ def run( # We're not passing all the unw files in, so we need to tally up below _unw_path, _cc_path = fut.result() - if zero_where_masked and mask_file is not None: + if zero_where_masked and mask_filename is not None: all_out_files = [ Path(str(outf).replace(UNW_SUFFIX, UNW_SUFFIX_ZEROED)) for outf in all_out_files @@ -219,7 +219,7 @@ def unwrap( corr_filename: Filename, unw_filename: Filename, nlooks: float, - mask_file: Optional[Filename] = None, + mask_filename: Optional[Filename] = None, zero_where_masked: bool = False, ntiles: Union[int, tuple[int, int]] = 1, tile_overlap: tuple[int, int] = (0, 0), @@ -250,7 +250,7 @@ def unwrap( Path to output unwrapped phase file. nlooks : float Effective number of looks used to form the input correlation data. - mask_file : Filename, optional + mask_filename : Filename, optional Path to binary byte mask file, by default None. Assumes that 1s are valid pixels and 0s are invalid. zero_where_masked : bool, optional @@ -306,9 +306,9 @@ def unwrap( Returns ------- - unw_path : Path + unw_path : Filename Path to output unwrapped phase file. - conncomp_path : Path + conncomp_path : Filename Path to output connected component label file. """ @@ -320,13 +320,13 @@ def unwrap( unwrap_method = UnwrapMethod(unwrap_method) # Check for a nodata mask - if io.get_raster_nodata(ifg_filename) is None or mask_file is None: + if io.get_raster_nodata(ifg_filename) is None or mask_filename is None: # With no marked `nodata`, just use the passed in mask - combined_mask_file = mask_file + combined_mask_file = mask_filename else: combined_mask_file = Path(ifg_filename).with_suffix(".mask.tif") create_combined_mask( - mask_filename=mask_file, + mask_filename=mask_filename, image_filename=ifg_filename, output_filename=combined_mask_file, ) @@ -395,6 +395,7 @@ def unwrap( f"Masking pixels with correlation below {interpolation_cor_threshold}" ) coherent_pixel_mask = corr[:] >= interpolation_cor_threshold + logger.info(f"Interpolating {pre_interp_ifg_filename} -> {interp_ifg_filename}") modified_ifg = interpolate( ifg=ifg, @@ -402,6 +403,7 @@ def unwrap( weight_cutoff=interpolation_cor_threshold, max_radius=max_radius, ) + logger.info(f"Writing interpolated output to {interp_ifg_filename}") io.write_arr( arr=modified_ifg, @@ -481,4 +483,19 @@ def unwrap( options=opts, ) + # Regrow connected components after phase modification + corr = io.load_gdal(corr_filename) + mask = corr[:] > 0 + # TODO decide whether we want to have the + # 'min_conncomp_frac' option in the config + conncomp_path = grow_conncomp_snaphu( + unw_filename=unw_filename, + corr_filename=corr_filename, + nlooks=nlooks, + mask=mask, + ccl_nodata=ccl_nodata, + cost=cost, + scratchdir=scratchdir, + ) + return unw_path, conncomp_path diff --git a/src/dolphin/workflows/_cli_config.py b/src/dolphin/workflows/_cli_config.py index 051d66d46..9139846e8 100644 --- a/src/dolphin/workflows/_cli_config.py +++ b/src/dolphin/workflows/_cli_config.py @@ -39,6 +39,8 @@ def create_config( no_unwrap: bool = False, no_inversion: bool = False, n_parallel_unwrap: int = 1, + run_goldstein: bool = False, + run_interpolation: bool = False, unwrap_method: UnwrapMethod = UnwrapMethod.SNAPHU, troposphere_files: Optional[list[str]] = None, tropo_date_fmt: str = "%Y%m%d", @@ -104,6 +106,8 @@ def create_config( "n_parallel_jobs": n_parallel_unwrap, "run_unwrap": not no_unwrap, "zero_where_masked": zero_where_masked, + "run_goldstein": run_goldstein, + "run_interpolation": run_interpolation, }, timeseries_options={ "run_inversion": not no_inversion, @@ -274,6 +278,16 @@ def get_parser(subparser=None, subcommand_name="run"): default=1, help="Number of interferograms to unwrap in parallel.", ) + unwrap_group.add_argument( + "--run-goldstein", + action="store_true", + help="Run Goldstein filter before unwrapping.", + ) + unwrap_group.add_argument( + "--run-interpolation", + action="store_true", + help="Run interpolation before unwrapping.", + ) # Correction options correction_group = parser.add_argument_group("Correction options") diff --git a/src/dolphin/workflows/unwrapping.py b/src/dolphin/workflows/unwrapping.py index f9c7a3796..2c099480f 100644 --- a/src/dolphin/workflows/unwrapping.py +++ b/src/dolphin/workflows/unwrapping.py @@ -78,7 +78,7 @@ def run( cor_filenames=cor_file_list, output_path=output_path, nlooks=nlooks, - mask_file=output_mask, + mask_filename=output_mask, zero_where_masked=unwrap_options.zero_where_masked, max_jobs=unwrap_options.n_parallel_jobs, ntiles=unwrap_options.ntiles,