From 46f5074d20b97b7f338f4c27f87b35350e3d20e8 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Thu, 24 Oct 2024 14:00:02 -0400 Subject: [PATCH 1/2] Update `dolphin timeseries` cli for new options --- src/dolphin/_cli_timeseries.py | 82 +++++++++++++++++++++++++++------- src/dolphin/timeseries.py | 6 +-- 2 files changed, 69 insertions(+), 19 deletions(-) diff --git a/src/dolphin/_cli_timeseries.py b/src/dolphin/_cli_timeseries.py index 4a09384a..70eddec4 100644 --- a/src/dolphin/_cli_timeseries.py +++ b/src/dolphin/_cli_timeseries.py @@ -1,7 +1,9 @@ import argparse +from datetime import datetime from pathlib import Path from typing import TYPE_CHECKING, Any +from dolphin.timeseries import InversionMethod from dolphin.workflows import CallFunc if TYPE_CHECKING: @@ -15,19 +17,17 @@ def get_parser(subparser=None, subcommand_name="timeseries") -> argparse.Argumen metadata = { "description": "Create a configuration file for a displacement workflow.", "formatter_class": argparse.ArgumentDefaultsHelpFormatter, - # https://docs.python.org/3/library/argparse.html#fromfile-prefix-chars "fromfile_prefix_chars": "@", } if subparser: - # Used by the subparser to make a nested command line interface parser = subparser.add_parser(subcommand_name, **metadata) else: parser = argparse.ArgumentParser(**metadata) # type: ignore[arg-type] - # parser._action_groups.pop() parser.add_argument( "-o", "--output-dir", + type=Path, default=Path(), help="Path to output directory to store results", ) @@ -58,30 +58,64 @@ def get_parser(subparser=None, subcommand_name="timeseries") -> argparse.Argumen parser.add_argument( "--condition-file", help=( - "A file with the same size as each raster, like amplitude dispersion or" - "temporal coherence to find reference point. default: amplitude dispersion" + "A file with the same size as each raster, like amplitude dispersion or " + "temporal coherence to find reference point" ), + required=True, ) parser.add_argument( "--condition", type=CallFunc, default=CallFunc.MIN, help=( - "A condition to apply to condition file to find the reference point" + "A condition to apply to condition file to find the reference point. " "Options are [min, max]. default=min" ), ) parser.add_argument( - "--num-threads", - type=int, - default=5, - help="Number of threads for the inversion", + "--method", + type=InversionMethod, + choices=list(InversionMethod), + default=InversionMethod.L1, + help=( + "Inversion method to use when solving Ax = b. L2 uses least squares" + " (faster), L1 minimizes |Ax - b|_1" + ), ) parser.add_argument( "--run-velocity", action="store_true", help="Run the velocity estimation from the phase time series", ) + parser.add_argument( + "--weight-velocity-by-corr", + action="store_true", + help=( + "Flag to indicate whether the velocity fitting should use correlation as" + " weights" + ), + ) + parser.add_argument( + "--correlation-threshold", + type=range_limited_float_type, + default=0.0, + metavar="[0-1]", + help="Pixels with correlation below this value will be masked out", + ) + parser.add_argument( + "--block-shape", + type=int, + nargs=2, + default=(256, 256), + metavar=("HEIGHT", "WIDTH"), + help="The shape of the blocks to process in parallel", + ) + parser.add_argument( + "--num-threads", + type=int, + default=4, + help="The parallel blocks to process at once", + ) parser.add_argument( "--reference-point", type=int, @@ -91,15 +125,31 @@ def get_parser(subparser=None, subcommand_name="timeseries") -> argparse.Argumen help=( "Reference point (row, col) used if performing a time series inversion. " "If not provided, a point will be selected from a consistent connected " - "component with low amplitude dispersion or high temporal coherence." + "component with low amplitude dispersion or high temporal coherence" ), ) parser.add_argument( - "--correlation-threshold", - type=range_limited_float_type, - default=0.2, - metavar="[0-1]", - help="Pixels with correlation below this value will be masked out.", + "--wavelength", + type=float, + help=( + "The wavelength of the radar signal, in meters. If provided, the output " + "rasters are in meters and meters/year for displacement and velocity. " + "If not provided, outputs are in radians" + ), + ) + parser.add_argument( + "--add-overviews", + action="store_true", + default=True, + help="If True, creates overviews of the new velocity raster", + ) + parser.add_argument( + "--extra-reference-date", + type=lambda s: datetime.strptime(s, "%Y%m%d"), + help=( + "If provided, makes another set of interferograms referenced to this " + "for all dates later than it. Format: YYYYMMDD" + ), ) parser.set_defaults(run_func=_run_timeseries) diff --git a/src/dolphin/timeseries.py b/src/dolphin/timeseries.py index e1c90109..f61928f2 100644 --- a/src/dolphin/timeseries.py +++ b/src/dolphin/timeseries.py @@ -41,16 +41,16 @@ class ReferencePointError(ValueError): def run( unwrapped_paths: Sequence[PathOrStr], - conncomp_paths: Sequence[PathOrStr], + conncomp_paths: Sequence[PathOrStr] | None, condition_file: PathOrStr, condition: CallFunc, output_dir: PathOrStr, - method: InversionMethod = InversionMethod.L2, + method: InversionMethod = InversionMethod.L1, run_velocity: bool = False, corr_paths: Sequence[PathOrStr] | None = None, weight_velocity_by_corr: bool = False, velocity_file: Optional[PathOrStr] = None, - correlation_threshold: float = 0.2, + correlation_threshold: float = 0.0, block_shape: tuple[int, int] = (256, 256), num_threads: int = 4, reference_point: tuple[int, int] = (-1, -1), From 48bb9b866161d686b0a7765da1c0170168bd7c28 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Thu, 24 Oct 2024 15:47:35 -0400 Subject: [PATCH 2/2] Add a `layover_shadow_mask_files` to mask pixels in layover/shadow during `wrapped_phase` --- src/dolphin/workflows/_utils.py | 4 ++ src/dolphin/workflows/config/_displacement.py | 8 ++++ src/dolphin/workflows/config/_ps.py | 8 ++++ src/dolphin/workflows/displacement.py | 7 ++++ src/dolphin/workflows/ps.py | 39 ++++++++++++------- src/dolphin/workflows/wrapped_phase.py | 38 ++++++++++-------- 6 files changed, 72 insertions(+), 32 deletions(-) diff --git a/src/dolphin/workflows/_utils.py b/src/dolphin/workflows/_utils.py index a2ab9b33..805f8591 100644 --- a/src/dolphin/workflows/_utils.py +++ b/src/dolphin/workflows/_utils.py @@ -22,6 +22,7 @@ def _create_burst_cfg( grouped_slc_files: dict[str, list[Path]], grouped_amp_mean_files: dict[str, list[Path]], grouped_amp_dispersion_files: dict[str, list[Path]], + grouped_layover_shadow_mask_files: dict[str, list[Path]], ) -> DisplacementWorkflow: cfg_temp_dict = cfg.model_dump(exclude={"cslc_file_list"}) @@ -31,6 +32,9 @@ def _create_burst_cfg( cfg_temp_dict["cslc_file_list"] = grouped_slc_files[burst_id] cfg_temp_dict["amplitude_mean_files"] = grouped_amp_mean_files[burst_id] cfg_temp_dict["amplitude_dispersion_files"] = grouped_amp_dispersion_files[burst_id] + cfg_temp_dict["layover_shadow_mask_files"] = grouped_layover_shadow_mask_files[ + burst_id + ] return DisplacementWorkflow(**cfg_temp_dict) diff --git a/src/dolphin/workflows/config/_displacement.py b/src/dolphin/workflows/config/_displacement.py index 9c86330f..4c4077f2 100644 --- a/src/dolphin/workflows/config/_displacement.py +++ b/src/dolphin/workflows/config/_displacement.py @@ -131,6 +131,14 @@ class DisplacementWorkflow(WorkflowBase): " calculation. If none provided, computed using the input SLC stack." ), ) + layover_shadow_mask_files: list[Path] = Field( + default_factory=list, + description=( + "Paths to layover/shadow binary masks, where 0 indicates a pixel in" + " layover/shadow, 1 is a good pixel. If none provided, no masking is" + " performed for layover/shadow." + ), + ) phase_linking: PhaseLinkingOptions = Field(default_factory=PhaseLinkingOptions) interferogram_network: InterferogramNetwork = Field( diff --git a/src/dolphin/workflows/config/_ps.py b/src/dolphin/workflows/config/_ps.py index c4fb9173..418a1004 100644 --- a/src/dolphin/workflows/config/_ps.py +++ b/src/dolphin/workflows/config/_ps.py @@ -37,6 +37,14 @@ class PsWorkflow(WorkflowBase): # Options for each step in the workflow ps_options: PsOptions = Field(default_factory=PsOptions) output_options: OutputOptions = Field(default_factory=OutputOptions) + layover_shadow_mask_files: list[Path] = Field( + default_factory=list, + description=( + "Paths to layover/shadow binary masks, where 0 indicates a pixel in" + " layover/shadow, 1 is a good pixel. If none provided, no masking is" + " performed for layover/shadow." + ), + ) # internal helpers model_config = ConfigDict( diff --git a/src/dolphin/workflows/displacement.py b/src/dolphin/workflows/displacement.py index 7950107b..eebc1247 100755 --- a/src/dolphin/workflows/displacement.py +++ b/src/dolphin/workflows/displacement.py @@ -90,6 +90,12 @@ def run( grouped_amp_mean_files = group_by_burst(cfg.amplitude_mean_files) else: grouped_amp_mean_files = defaultdict(list) + if cfg.layover_shadow_mask_files: + grouped_layover_shadow_mask_files = group_by_burst( + cfg.layover_shadow_mask_files + ) + else: + grouped_layover_shadow_mask_files = defaultdict(list) grouped_iono_files = parse_ionosphere_files( cfg.correction_options.ionosphere_files, cfg.correction_options._iono_date_fmt @@ -109,6 +115,7 @@ def run( grouped_slc_files, grouped_amp_mean_files, grouped_amp_dispersion_files, + grouped_layover_shadow_mask_files, ), ) for burst in grouped_slc_files diff --git a/src/dolphin/workflows/ps.py b/src/dolphin/workflows/ps.py index 91320771..59c34278 100644 --- a/src/dolphin/workflows/ps.py +++ b/src/dolphin/workflows/ps.py @@ -7,19 +7,20 @@ import opera_utils import dolphin.ps -from dolphin import __version__ +from dolphin import __version__, masking from dolphin._log import log_runtime, setup_logging from dolphin.io import VRTStack from dolphin.utils import get_max_memory_usage +from dolphin.workflows.wrapped_phase import _get_mask -from .config import PsWorkflow +from .config import DisplacementWorkflow, PsWorkflow logger = logging.getLogger(__name__) @log_runtime def run( - cfg: PsWorkflow, + cfg: PsWorkflow | DisplacementWorkflow, compute_looked: bool = False, debug: bool = False, ) -> list[Path]: @@ -70,25 +71,32 @@ def run( msg = "No input files found" raise ValueError(msg) - # ############################################# - # Make a VRT pointing to the input SLC files - # ############################################# subdataset = cfg.input_options.subdataset vrt_stack = VRTStack( input_file_list, subdataset=subdataset, outfile=cfg.work_directory / "slc_stack.vrt", ) + # Mark any files beginning with "compressed" as compressed + is_compressed = ["compressed" in str(f).lower() for f in input_file_list] - # Make the nodata mask from the polygons, if we're using OPERA CSLCs - try: - nodata_mask_file = cfg.work_directory / "nodata_mask.tif" - opera_utils.make_nodata_mask( - vrt_stack.file_list, out_file=nodata_mask_file, buffer_pixels=200 - ) - except Exception as e: - logger.warning(f"Could not make nodata mask: {e}") - nodata_mask_file = None + non_compressed_slcs = [ + f for f, is_comp in zip(input_file_list, is_compressed) if not is_comp + ] + + layover_shadow_mask = ( + cfg.layover_shadow_mask_files[0] if cfg.layover_shadow_mask_files else None + ) + mask_filename = _get_mask( + output_dir=cfg.work_directory, + output_bounds=cfg.output_options.bounds, + output_bounds_wkt=cfg.output_options.bounds_wkt, + output_bounds_epsg=cfg.output_options.bounds_epsg, + like_filename=vrt_stack.outfile, + layover_shadow_mask=layover_shadow_mask, + cslc_file_list=non_compressed_slcs, + ) + nodata_mask = masking.load_mask_as_numpy(mask_filename) if mask_filename else None logger.info(f"Creating persistent scatterer file {ps_output}") dolphin.ps.create_ps( @@ -98,6 +106,7 @@ def run( output_amp_dispersion_file=output_file_list[2], like_filename=vrt_stack.outfile, amp_dispersion_threshold=cfg.ps_options.amp_dispersion_threshold, + nodata_mask=nodata_mask, block_shape=cfg.worker_settings.block_shape, ) # Save a looked version of the PS mask too diff --git a/src/dolphin/workflows/wrapped_phase.py b/src/dolphin/workflows/wrapped_phase.py index 0e1ec082..31008cf5 100644 --- a/src/dolphin/workflows/wrapped_phase.py +++ b/src/dolphin/workflows/wrapped_phase.py @@ -80,7 +80,9 @@ def run( non_compressed_slcs = [ f for f, is_comp in zip(input_file_list, is_compressed) if not is_comp ] - + layover_shadow_mask = ( + cfg.layover_shadow_mask_files[0] if cfg.layover_shadow_mask_files else None + ) # Create a mask file from input bounding polygons and/or specified output bounds mask_filename = _get_mask( output_dir=cfg.work_directory, @@ -88,6 +90,7 @@ def run( output_bounds_wkt=cfg.output_options.bounds_wkt, output_bounds_epsg=cfg.output_options.bounds_epsg, like_filename=vrt_stack.outfile, + layover_shadow_mask=layover_shadow_mask, cslc_file_list=non_compressed_slcs, ) @@ -461,9 +464,11 @@ def _get_mask( output_bounds_wkt: str | None, output_bounds_epsg: int, like_filename: Filename, + layover_shadow_mask: Filename | None, cslc_file_list: Sequence[Filename], ) -> Path | None: # Make the nodata mask from the polygons, if we're using OPERA CSLCs + mask_files: list[Path] = [] try: nodata_mask_file = output_dir / "nodata_mask.tif" @@ -472,11 +477,11 @@ def _get_mask( out_file=nodata_mask_file, buffer_pixels=200, ) + mask_files.append(nodata_mask_file) except Exception as e: logger.warning(f"Could not make nodata mask: {e}") nodata_mask_file = None - mask_filename: Path | None = None # Also mask outside the area of interest if we've specified a small bounds if output_bounds is not None or output_bounds_wkt is not None: # Make a mask just from the bounds @@ -488,23 +493,22 @@ def _get_mask( output_filename=bounds_mask_filename, like_filename=like_filename, ) + mask_files.append(bounds_mask_filename) - # Then combine with the nodata mask - if nodata_mask_file is not None: - combined_mask_filename = output_dir / "combined_mask.tif" - if not combined_mask_filename.exists(): - masking.combine_mask_files( - mask_files=[bounds_mask_filename, nodata_mask_file], - output_file=combined_mask_filename, - output_convention=masking.MaskConvention.ZERO_IS_NODATA, - ) - mask_filename = combined_mask_filename - else: - mask_filename = bounds_mask_filename - else: - mask_filename = nodata_mask_file + if layover_shadow_mask is not None: + mask_files.append(Path(layover_shadow_mask)) + + if not mask_files: + return None - return mask_filename + combined_mask_filename = output_dir / "combined_mask.tif" + if not combined_mask_filename.exists(): + masking.combine_mask_files( + mask_files=mask_files, + output_file=combined_mask_filename, + output_convention=masking.MaskConvention.ZERO_IS_NODATA, + ) + return combined_mask_filename def _is_single_reference_network(