Add phase_similarity raster, add water masking to recommended_mask

src/disp_s1/main.py

@@ -9,7 +9,7 @@
 from pathlib import Path
 from typing import NamedTuple
 
-from dolphin import interferogram
+from dolphin import interferogram, stitching
 from dolphin._log import log_runtime, setup_logging
 from dolphin.io import load_gdal
 from dolphin.unwrap import grow_conncomp_snaphu
@@ -47,12 +47,14 @@ def run(
 
     """
     setup_logging(logger_name="disp_s1", debug=debug, filename=cfg.log_file)
+    cfg.work_directory.mkdir(exist_ok=True, parents=True)
     # Setup the binary mask as dolphin expects
     if pge_runconfig.dynamic_ancillary_file_group.mask_file:
         water_binary_mask = cfg.work_directory / "water_binary_mask.tif"
         create_mask_from_distance(
             water_distance_file=pge_runconfig.dynamic_ancillary_file_group.mask_file,
             output_file=water_binary_mask,
+            # Set a little conservative for the general processing
             land_buffer=2,
             ocean_buffer=2,
         )
@@ -127,6 +129,29 @@ def run(
             disp_date_keys, out_paths.ionospheric_corrections, "ionosphere"
         )
 
+    if pge_runconfig.dynamic_ancillary_file_group.mask_file:
+        aggressive_water_binary_mask = (
+            cfg.work_directory / "water_binary_mask_nobuffer.tif"
+        )
+        tmp_outfile = aggressive_water_binary_mask.with_suffix(".temp.tif")
+        create_mask_from_distance(
+            water_distance_file=pge_runconfig.dynamic_ancillary_file_group.mask_file,
+            # Make the file in lat/lon
+            output_file=tmp_outfile,
+            # Still don't trust the land water 100%
+            land_buffer=1,
+            # Trust the ocean buffer
+            ocean_buffer=0,
+        )
+        # Then need to warp to match the output UTM files
+        stitching.warp_to_match(
+            input_file=tmp_outfile,
+            match_file=out_paths.timeseries_paths[0],
+            output_file=aggressive_water_binary_mask,
+        )
+    else:
+        aggressive_water_binary_mask = None
+
     logger.info(f"Creating {len(out_paths.timeseries_paths)} outputs in {out_dir}")
     create_displacement_products(
         out_paths,
@@ -138,6 +163,7 @@ def run(
         reference_point=ref_point,
         los_east_file=los_east_file,
         los_north_file=los_north_file,
+        water_mask=aggressive_water_binary_mask,
     )
     logger.info("Finished creating output products.")
 
@@ -197,6 +223,8 @@ class ProductFiles(NamedTuple):
     troposphere: Path | None
     ionosphere: Path | None
     unwrapper_mask: Path | None
+    similarity: Path
+    water_mask: Path | None
 
 
 def process_product(
@@ -277,6 +305,8 @@ def process_product(
         ps_mask_filename=files.ps_mask,
         shp_count_filename=files.shp_counts,
         unwrapper_mask_filename=files.unwrapper_mask,
+        similarity_filename=files.similarity,
+        water_mask_filename=files.water_mask,
         los_east_file=los_east_file,
         los_north_file=los_north_file,
         pge_runconfig=pge_runconfig,
@@ -297,10 +327,11 @@ def create_displacement_products(
     date_to_cslc_files: Mapping[tuple[datetime], list[Path]],
     pge_runconfig: RunConfig,
     dolphin_config: DisplacementWorkflow,
-    wavelength_cutoff: float = 50_000.0,
+    wavelength_cutoff: float = 25_000.0,
     reference_point: ReferencePoint | None = None,
     los_east_file: Path | None = None,
     los_north_file: Path | None = None,
+    water_mask: Path | None = None,
     max_workers: int = 3,
 ) -> None:
     """Run parallel processing for all interferograms.
@@ -322,14 +353,17 @@ def create_displacement_products(
         If None, will record empty in the dataset's attributes
     wavelength_cutoff : float
         Wavelength cutoff (in meters) for filtering long wavelengths.
-        Default is 50_000.
+        Default is 25_000.
     reference_point : ReferencePoint, optional
         Reference point recorded from dolphin after unwrapping.
         If none, leaves product attributes empty.
     los_east_file : Path, optional
         Path to the east component of line of sight unit vector
     los_north_file : Path, optional
         Path to the north component of line of sight unit vector
+    water_mask : Path, optional
+        Binary water mask to use for output product.
+        If provided, is used in the `recommended_mask`.
     max_workers : int
         Number of parallel products to process.
         Default is 3.
@@ -360,6 +394,8 @@ def create_displacement_products(
             troposphere=tropo,
             ionosphere=iono,
             unwrapper_mask=mask_f,
+            similarity=out_paths.stitched_similarity_file,
+            water_mask=water_mask,
         )
         for unw, cc, cor, tropo, iono, mask_f in zip(
             out_paths.timeseries_paths,

src/disp_s1/pge_runconfig.py

@@ -191,7 +191,7 @@ class AlgorithmParameters(YamlModel):
         description="Name of the subdataset to use in the input NetCDF files.",
     )
     spatial_wavelength_cutoff: float = Field(
-        50_000,
+        25_000,
         description=(
             "Spatial wavelength cutoff (in meters) for the spatial filter. Used to"
             " create the short wavelength displacement layer"

src/disp_s1/product.py

@@ -85,6 +85,8 @@ def create_output_product(
     ps_mask_filename: Filename,
     shp_count_filename: Filename,
     unwrapper_mask_filename: Filename | None,
+    similarity_filename: Filename,
+    water_mask_filename: Filename | None,
     pge_runconfig: RunConfig,
     dolphin_config: DisplacementWorkflow,
     reference_cslc_files: list[Filename],
@@ -93,7 +95,7 @@ def create_output_product(
     los_north_file: Filename | None = None,
     reference_point: ReferencePoint | None = None,
     corrections: Optional[dict[str, ArrayLike]] = None,
-    wavelength_cutoff: float = 50_000.0,
+    wavelength_cutoff: float = 25_000.0,
 ):
     """Create the OPERA output product in NetCDF format.
 
@@ -115,6 +117,10 @@ def create_output_product(
         The path to statistically homogeneous pixels (SHP) counts.
     unwrapper_mask_filename : Filename, optional
         The path to the boolean mask used during unwrapping to ignore pixels.
+    similarity_filename : Filename
+        The path to the cosine similarity image.
+    water_mask_filename : Filename, optional
+        Path to the binary water mask to use in creating a recommended mask.
     pge_runconfig : Optional[RunConfig], optional
         The PGE run configuration, by default None
         Used to add extra metadata to the output file.
@@ -137,7 +143,7 @@ def create_output_product(
         A dictionary of corrections to write to the output file, by default None
     wavelength_cutoff : float, optional
         The wavelength cutoff for filtering long wavelengths.
-        Default is 50_000.0
+        Default is 25_000.0
 
 
     """
@@ -225,12 +231,30 @@ def _get_start_end_cslcs(files):
         "Creating short wavelength displacement product with %s meter cutoff",
         wavelength_cutoff,
     )
+
+    # Create the commended mask:
+    temporal_coherence = io.load_gdal(temp_coh_filename, masked=True).mean()
+    bad_temporal_coherence = temporal_coherence < 0.6
+    # Get summary statistics on the layers for CMR filtering/searching purposes
+    average_temporal_coherence = temporal_coherence.mean()
+
+    if water_mask_filename:
+        water_mask_data = io.load_gdal(water_mask_filename, masked=True).filled(0)
+        is_water = water_mask_data == 0
+    else:
+        # Not provided: Don't indicate anything is water in this mask.
+        is_water = np.zeros(temporal_coherence.shape, dtype=bool)
+
     conncomps = io.load_gdal(conncomp_filename, masked=True).filled(0)
     bad_conncomp = conncomps == 0
 
-    temporal_coherence = io.load_gdal(temp_coh_filename, masked=True).mean()
-    bad_temporal_coherence = temporal_coherence < 0.5
-    bad_pixel_mask = bad_temporal_coherence | bad_conncomp
+    similarity = io.load_gdal(similarity_filename, masked=True).mean()
+    bad_similarity = similarity < 0.5
+    bad_pixel_mask = is_water | bad_conncomp | (bad_temporal_coherence & bad_similarity)
+    # Note: An alternate way to view this:
+    # good_conncomp & is_no_water & (good_temporal_coherence | good_similarity)
+    recommended_mask = ~bad_pixel_mask
+    del temporal_coherence, conncomps, similarity
 
     filtered_disp_arr = filtering.filter_long_wavelength(
         unwrapped_phase=disp_arr,
@@ -245,8 +269,6 @@ def _get_start_end_cslcs(files):
     disp_arr[mask] = np.nan
     filtered_disp_arr[mask] = np.nan
 
-    recommended_mask = ~bad_pixel_mask
-
     product_infos: list[ProductInfo] = list(DISPLACEMENT_PRODUCTS)
 
     with h5netcdf.File(output_name, "w", **FILE_OPTS) as f:
@@ -296,16 +318,16 @@ def _get_start_end_cslcs(files):
 
         # For the others, load and save each individually
         data_files = [
-            conncomp_filename,
             conncomp_filename,
             temp_coh_filename,
             ifg_corr_filename,
             ps_mask_filename,
             shp_count_filename,
             unwrapper_mask_filename,
+            similarity_filename,
         ]
 
-        for info, filename in zip(product_infos[3:], data_files):
+        for info, filename in zip(product_infos[3:], data_files, strict=True):
             if filename is not None and Path(filename).exists():
                 data = io.load_gdal(filename).astype(info.dtype)
             else:
@@ -323,7 +345,6 @@ def _get_start_end_cslcs(files):
                 fillvalue=info.fillvalue,
                 attrs=info.attrs,
             )
-            del data  # Free up memory
 
     if los_east_file is not None and los_north_file is not None:
         logger.info("Calculating solid earth tide")
@@ -354,8 +375,6 @@ def _get_start_end_cslcs(files):
         reference_point=reference_point,
     )
 
-    # Get summary statistics on the layers for CMR filtering/searching purposes
-    average_temporal_coherence = temporal_coherence.mean()
     _create_identification_group(
         output_name=output_name,
         pge_runconfig=pge_runconfig,

src/disp_s1/product_info.py

@@ -139,6 +139,21 @@ class DisplacementProducts:
             dtype=np.uint8,
         )
     )
+    phase_similarity: ProductInfo = field(
+        default_factory=lambda: ProductInfo(
+            name="phase_similarity",
+            long_name="Phase Similarity",
+            description=(
+                "Median cosine similarity of wrapped phase in each pixel's"
+                " neighborhood."
+            ),
+            fillvalue=np.nan,
+            attrs={"units": "unitless"},
+            # 8 bits (between 0 and 1) is around .001 precision
+            keep_bits=8,
+            dtype=np.float32,
+        )
+    )
 
     def __iter__(self):
         """Return all displacement dataset info as an iterable."""