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(