Repository: isce-framework/dolphin (language: python)

Issue:
Create the single-update version of the s1_disp workflow
Currently we have a version of the ministack-based workflow (`dolphin/workflows/s1_disp_stack.py`) to process a historical archive of data. We need the ability to ingest one new SLC, pull some number of old SLCs + saved compressed SLCs to produce a new wrapped phase estimate.

The structure of the core modules (`ps.py` and the 'phase_link` subpackage) should make it straighforward to add the single-update version; it will be mostly bookkeeping on where/how to pull the required data.

The general outline of this new workflow should be

- Locate the SLCs  and the Compressed SLCs we will use
- Create a VRTStack https://github.com/opera-adt/dolphin/blob/eda6c25057568cf421c8937f0e01a69c774ba514/src/dolphin/workflows/wrapped_phase.py#L37-L41 which should sort the compressed SLCs before the normal ones
- (possibly deferred until later): update the pre-computed PS map: https://github.com/opera-adt/dolphin/blob/eda6c25057568cf421c8937f0e01a69c774ba514/src/dolphin/ps.py#L146-L151 . Re-save to some location
- Pass the PS pixels and VRT stack through the `phase_link.run_mle` : https://github.com/opera-adt/dolphin/blob/eda6c25057568cf421c8937f0e01a69c774ba514/src/dolphin/sequential.py#L147-L178
    - Note that we won't need to use the `sequential.py` module; we're only passing one stack through `run_mle`, so we don't need to do all the mini-stack juggling
    - We also only need to save the last layer that will come out of `run_mle`, since it a number of images equal to the number of input SLCs
- Reference the newly estimated SLC phase to something (TODO: make a separate issue for that research topic)

This will probably be divided into several issues for easier testing/tracking.


Discussion / hints:
The structure for this was created in #47, but we still need to add some kind of testing script to repeatedly run this and ensure we get fine results for a timeseries.

Fix patch:
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index a5864928..3b50c06d 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -6,6 +6,13 @@ version: 2
 mkdocs:
   configuration: mkdocs.yml
 
+# Conda is causing "excessive memory usage"
+# https://docs.readthedocs.io/en/stable/guides/conda.html#making-builds-faster-with-mamba
+build:
+  os: "ubuntu-20.04"
+  tools:
+    python: "mambaforge-4.10"
+
 conda:
   environment: conda-env.yml
 
diff --git a/CHANGELOG.md b/CHANGELOG.md
index b4974486..14b66a3f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,11 +13,13 @@
 - Also included a script to view the window in an interactive matplotlib figure (matplotlib must be installed separately)
 - Added a simple method to check for adjacent-pixel unwrapping errors in `unwrap.compute_phase_diffs`
 - Adds a method `utils.get_cpu_count` which returns either `os.cpu_count`, or (if running in a Docker container) the number of CPUs allocated by Docker
+- Created a script to test the incremental/near-real-time version of phase linking
 
 **Changes**
 
 - Passing an existing file to `VRTStack` will no longer error unless `fail_on_overwrite=True`. The default just prints out the overwrite is happening. This prevents multiple runs in the same folder from errorings just for creating a reference to the SLC files.
 - The environment variable `NUMBA_NUM_THREADS` is set using the passed in config to prevent numba from using all CPUs during `prange` calls
+- The `sequential.py` module uses a different implementation of the sequential estimator to avoid the need for a datum adjustment.
 
 **Dependencies**
 
@@ -93,7 +95,6 @@ Added requirements:
 
 Added requirements:
 
-- pyproj>=3.2
 - rich>=12.0
 - threadpoolctl>=3.0
 - isce3>=0.8.0
diff --git a/conda-env.yml b/conda-env.yml
index 527425cc..96a51cbd 100644
--- a/conda-env.yml
+++ b/conda-env.yml
@@ -15,7 +15,7 @@ dependencies:
   - pillow>=7.0
   - pydantic>=1.10
   - pymp-pypi>=0.4.5
-  - pyproj>=3.2
+  - pyproj>=3.3
   - rich>=12.0
   - ruamel.yaml>=0.15
   - scipy>=1.5 # "scipy 0.16+ is required for linear algebra", numba. 1.5 is the oldest version that supports Python 3.7
diff --git a/requirements.txt b/requirements.txt
index 95503275..5c081cef 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -6,7 +6,7 @@ numpy>=1.20
 pillow>=7.0
 pydantic>=1.10
 pymp-pypi>=0.4.5
-pyproj>=3.2
+pyproj>=3.3
 rich>=12.0
 ruamel_yaml>=0.15
 scipy>=1.5
diff --git a/scripts/run_repeated_nrt.py b/scripts/run_repeated_nrt.py
new file mode 100755
index 00000000..6323960c
--- /dev/null
+++ b/scripts/run_repeated_nrt.py
@@ -0,0 +1,421 @@
+#!/usr/bin/env python
+from __future__ import annotations
+
+import argparse
+import os
+from concurrent.futures import ThreadPoolExecutor  # , as_completed
+from itertools import chain
+from pathlib import Path
+from typing import Sequence
+
+from dolphin import io, ps, stack, utils
+from dolphin._log import get_log, log_runtime
+from dolphin._types import Filename
+from dolphin.workflows import s1_disp
+from dolphin.workflows._utils import group_by_burst, make_nodata_mask
+from dolphin.workflows.config import (
+    OPERA_DATASET_NAME,
+    InterferogramNetworkType,
+    ShpMethod,
+    Workflow,
+    WorkflowName,
+)
+
+logger = get_log("dolphin.run_repeated_nrt")
+
+
+def _create_cfg(
+    *,
+    slc_files: Sequence[Filename],
+    half_window_size: tuple[int, int] = (11, 5),
+    first_ministack: bool = False,
+    run_unwrap: bool = False,
+    shp_method: ShpMethod = ShpMethod.GLRT,
+    amplitude_mean_files=[],
+    amplitude_dispersion_files=[],
+):
+    strides = {"x": 6, "y": 3}
+    if first_ministack:
+        interferogram_network = dict(
+            network_type=InterferogramNetworkType.SINGLE_REFERENCE
+        )
+        workflow_name = WorkflowName.STACK
+    else:
+        interferogram_network = dict(
+            network_type=InterferogramNetworkType.MANUAL_INDEX,
+            indexes=[(0, -1)],
+        )
+        workflow_name = WorkflowName.SINGLE
+
+    cfg = Workflow(
+        # Things that change with each workflow run
+        cslc_file_list=slc_files,
+        interferogram_network=interferogram_network,
+        amplitude_mean_files=amplitude_mean_files,
+        amplitude_dispersion_files=amplitude_dispersion_files,
+        # Configurable from CLI inputs:
+        output_options=dict(
+            strides=strides,
+        ),
+        phase_linking=dict(
+            ministack_size=1000,  # for single update, process in one ministack
+            half_window={"x": half_window_size[0], "y": half_window_size[1]},
+            shp_method=shp_method,
+        ),
+        # worker_settings=dict(
+        #     block_size_gb=block_size_gb,
+        #     n_workers=n_workers,
+        #     threads_per_worker=threads_per_worker,
+        # ),
+        #     ps_options=dict(
+        #         amp_dispersion_threshold=amp_dispersion_threshold,
+        #     ),
+        #     log_file=log_file,
+        # )
+        # Definite hard coded things
+        unwrap_options=dict(
+            unwrap_method="snaphu",
+            run_unwrap=run_unwrap,
+            # CHANGEME: or else run in backgroun somehow?
+        ),
+        save_compressed_slc=True,  # always save, and only sometimes will we grab it
+        workflow_name=workflow_name,
+    )
+    return cfg
+
+
+def get_cli_args() -> argparse.Namespace:
+    """Set up the command line interface."""
+    parser = argparse.ArgumentParser(
+        description="Repeatedly run the dolphin single update mode.",
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+        # https://docs.python.org/3/library/argparse.html#fromfile-prefix-chars
+        fromfile_prefix_chars="@",
+    )
+    parser.add_argument(
+        "--slc-files",
+        nargs=argparse.ONE_OR_MORE,
+        help="List the paths of all SLC files to include.",
+        required=True,
+    )
+    parser.add_argument(
+        "-ms",
+        "--ministack-size",
+        type=int,
+        default=15,
+        help="Strides/decimation factor (x, y) (in pixels) to use when determining",
+    )
+    parser.add_argument(
+        "-hw",
+        "--half-window-size",
+        type=int,
+        nargs=2,
+        default=(11, 5),
+        metavar=("X", "Y"),
+        help="Half window size for the phase linking algorithm",
+    )
+    parser.add_argument(
+        "--shp-method",
+        type=ShpMethod,
+        choices=[s.value for s in ShpMethod],
+        default=ShpMethod.GLRT,
+        help="Method used to calculate the SHP.",
+    )
+    parser.add_argument(
+        "--run-unwrap",
+        action="store_true",
+        help="Run the unwrapping stack after phase linking.",
+    )
+    parser.add_argument(
+        "--pre-compute",
+        action="store_true",
+        help=(
+            "Run the amplitude mean/dispersion pre-compute step (not the main"
+            " workflow)."
+        ),
+    )
+    return parser.parse_args()
+
+
+@log_runtime
+def compute_ps_files(
+    burst_grouped_slc_files: dict[str, list[Filename]],
+    burst_to_nodata_mask: dict[str, Filename],
+    # max_workers: int = 3,
+    ps_stack_size: int = 60,
+    output_folder: Path = Path("precomputed_ps"),
+):
+    """Compute the mean/DA/PS files for each burst group."""
+    all_amp_files, all_disp_files, all_ps_files = [], [], []
+    # future_burst_dict = {}
+    # with ThreadPoolExecutor(max_workers=max_workers) as exc:
+    for burst, file_list in burst_grouped_slc_files.items():
+        nodata_mask_file = burst_to_nodata_mask[burst]
+        # fut = exc.submit(
+        # _compute_burst_ps_files,
+        amp_files, disp_files, ps_files = _compute_burst_ps_files(
+            burst,
+            file_list,
+            nodata_mask_file=nodata_mask_file,
+            ps_stack_size=ps_stack_size,
+            output_folder=output_folder,
+            # show_progress=False,
+        )
+        # future_burst_dict[fut] = burst
+
+        # for future in as_completed(future_burst_dict.keys()):
+        # burst = future_burst_dict[future]
+        # amp_files, disp_files, ps_files = future.result()
+
+        all_amp_files.extend(amp_files)
+        all_disp_files.extend(disp_files)
+        all_ps_files.extend(ps_files)
+
+        logger.info(f"Done with {burst}")
+
+    return all_amp_files, all_disp_files, all_ps_files
+
+
+@log_runtime
+def _compute_burst_ps_files(
+    burst: str,
+    file_list_all: list[Filename],
+    nodata_mask_file: Filename,
+    ps_stack_size: int = 60,
+    output_folder: Path = Path("precomputed_ps"),
+) -> tuple[list[Path], list[Path], list[Path]]:
+    """Pre-compute the PS files (mean / amp. dispersion) for one burst."""
+    logger.info(f"Computing PS files for {burst} into {output_folder}")
+    vrt_all = stack.VRTStack(
+        file_list_all, subdataset=OPERA_DATASET_NAME, write_file=False
+    )
+    # logger.info("Created total vrt")
+    date_list_all = vrt_all.dates
+
+    nodata_mask = io.load_gdal(nodata_mask_file, masked=True).astype(bool).filled(False)
+    # invert the mask so 1s are the missing data pixels
+    nodata_mask = ~nodata_mask
+
+    # TODO: fixed number of PS files? fixed time window?
+    amp_files, disp_files, ps_files = [], [], []
+    for full_stack_idx in range(0, len(file_list_all), ps_stack_size):
+        cur_slice = slice(full_stack_idx, full_stack_idx + ps_stack_size)
+        cur_files = file_list_all[cur_slice]
+        cur_dates = date_list_all[cur_slice]
+
+        # Make the current ministack output folder using the start/end dates
+        d0 = cur_dates[0][0]
+        d1 = cur_dates[-1][0]
+        start_end = io._format_date_pair(d0, d1)
+        basename = f"{burst}_{start_end}"
+
+        # output_folder = output_folder / start_end
+        output_folder.mkdir(parents=True, exist_ok=True)
+        cur_vrt = stack.VRTStack(
+            cur_files,
+            outfile=output_folder / f"{basename}.vrt",
+            subdataset=OPERA_DATASET_NAME,
+        )
+        cur_ps_file = (output_folder / f"{basename}_ps_pixels.tif").resolve()
+        cur_amp_mean = (output_folder / f"{basename}_amp_mean.tif").resolve()
+        cur_amp_dispersion = (
+            output_folder / f"{basename}_amp_dispersion.tif"
+        ).resolve()
+        if not all(f.exists() for f in [cur_ps_file, cur_amp_mean, cur_amp_dispersion]):
+            ps.create_ps(
+                slc_vrt_file=cur_vrt,
+                output_amp_mean_file=cur_amp_mean,
+                output_amp_dispersion_file=cur_amp_dispersion,
+                output_file=cur_ps_file,
+                nodata_mask=nodata_mask,
+                block_size_gb=0.2,
+            )
+        else:
+            logger.info(f"Skipping existing {basename} files in {output_folder}")
+        amp_files.append(cur_amp_mean)
+        disp_files.append(cur_amp_dispersion)
+        ps_files.append(cur_ps_file)
+        logger.info(f"Finished with PS processing for {burst}")
+    return amp_files, disp_files, ps_files
+
+
+def create_nodata_masks(
+    # date_grouped_slc_files: dict[tuple[datetime.date], list[Filename]],
+    burst_grouped_slc_files: dict[str, list[Filename]],
+    buffer_pixels: int = 30,
+    max_workers: int = 3,
+    output_folder: Path = Path("nodata_masks"),
+):
+    """Create the nodata binary masks for each burst."""
+    output_folder.mkdir(exist_ok=True, parents=True)
+    futures = []
+    out_burst_to_file = {}
+    with ThreadPoolExecutor(max_workers=max_workers) as exc:
+        for burst, file_list in burst_grouped_slc_files.items():
+            outfile = output_folder / f"{burst}.tif"
+            fut = exc.submit(
+                make_nodata_mask,
+                file_list,
+                outfile,
+                buffer_pixels=buffer_pixels,
+            )
+            futures.append(fut)
+            out_burst_to_file[burst] = outfile
+        for fut in futures:
+            fut.result()
+    return out_burst_to_file
+
+
+def _form_burst_vrt_stacks(
+    burst_grouped_slc_files: dict[str, list[Filename]]
+) -> dict[str, stack.VRTStack]:
+    logger.info("For each burst, creating a VRTStack...")
+    # Each burst needs to be the same size
+    burst_to_vrt_stack = {}
+    for b, file_list in burst_grouped_slc_files.items():
+        logger.info(f"Checking {len(file_list)} files for {b}")
+        outfile = Path(f"slc_stack_{b}.vrt")
+        if not outfile.exists():
+            vrt = stack.VRTStack(
+                file_list, subdataset=OPERA_DATASET_NAME, outfile=outfile
+            )
+        else:
+            vrt = stack.VRTStack.from_vrt_file(outfile)
+        burst_to_vrt_stack[b] = vrt
+    logger.info("Done.")
+    return burst_to_vrt_stack
+
+
+@log_runtime
+def precompute_ps_files(arg_dict) -> None:
+    """Run the pre-compute step to get means/amp. dispersion for each burst."""
+    all_slc_files = arg_dict.pop("slc_files")
+
+    burst_grouped_slc_files = group_by_burst(all_slc_files)
+    #  {'t173_370312_iw2': [PosixPath('t173_370312_iw2_20170203.h5'),... ] }
+    date_grouped_slc_files = utils.group_by_date(all_slc_files)
+    #  { (datetime.date(2017, 5, 22),) : [PosixPath('t173_370311_iw1_20170522.h5'), ] }
+    logger.info(f"Found {len(all_slc_files)} total SLC files")
+    logger.info(f"  {len(date_grouped_slc_files)} unique dates,")
+    logger.info(f"  {len(burst_grouped_slc_files)} unique bursts.")
+
+    burst_to_nodata_mask = create_nodata_masks(burst_grouped_slc_files)
+
+    all_amp_files, all_disp_files, all_ps_files = compute_ps_files(
+        burst_grouped_slc_files, burst_to_nodata_mask
+    )
+
+
+def _get_all_slc_files(burst_to_vrt_stack, start_idx, end_idx):
+    return list(
+        chain.from_iterable(
+            [v.file_list[start_idx:end_idx] for v in burst_to_vrt_stack.values()]
+        )
+    )
+
+
+@log_runtime
+def main(arg_dict: dict) -> None:
+    """Get the command line arguments and run the workflow."""
+    arg_dict = vars(args)
+    ministack_size = arg_dict.pop("ministack_size")
+    # TODO: verify this is fine to sort them by date?
+    all_slc_files = sorted(arg_dict.pop("slc_files"))
+    logger.info(f"Found {len(all_slc_files)} total SLC files")
+
+    # format of `group_by_burst`:
+    #   {'t173_370312_iw2': [PosixPath('t173_370312_iw2_20170203.h5'),... ] }
+    burst_grouped_slc_files = group_by_burst(all_slc_files)
+    num_bursts = len(burst_grouped_slc_files)
+    logger.info(f"  {num_bursts} unique bursts.")
+    # format of `group_by_date`:
+    #  { (datetime.date(2017, 5, 22),) : [PosixPath('t173_370311_iw1_20170522.h5'), ] }
+    date_grouped_slc_files = utils.group_by_date(all_slc_files)
+    num_dates = len(date_grouped_slc_files)
+    logger.info(f"  {num_dates} unique dates,")
+
+    burst_to_vrt_stack = _form_burst_vrt_stacks(
+        burst_grouped_slc_files=burst_grouped_slc_files
+    )
+
+    # Get the pre-compted PS files (assuming --pre-compute has been run)
+    all_amp_files = sorted(Path("precomputed_ps/").resolve().glob("*_amp_mean.tif"))
+    all_disp_files = sorted(
+        Path("precomputed_ps/").resolve().glob("*_amp_dispersion.tif")
+    )
+    # TODO: how to make it shift when the year changes for PS files
+
+    slc_idx_start = 0
+    slc_idx_end = ministack_size
+    max_stack_size = 2 * ministack_size - 1  # Size at which we archive/shrink
+    cur_path = Path(f"stack_{slc_idx_start}_{slc_idx_end}")
+    cur_path.mkdir(exist_ok=True)
+    # TODO: use the context-manager change-then-back
+    os.chdir(cur_path)
+
+    # TODO: do i wanna somehow provide the burst nodata file?
+
+    # Make the first ministack
+    cur_slc_files = _get_all_slc_files(burst_to_vrt_stack, slc_idx_start, slc_idx_end)
+    cfg = _create_cfg(
+        slc_files=cur_slc_files,
+        first_ministack=True,
+        amplitude_mean_files=all_amp_files,
+        amplitude_dispersion_files=all_disp_files,
+        **arg_dict,
+    )
+
+    cfg.to_yaml("dolphin_config.yaml")
+    s1_disp.run(cfg)
+    os.chdir("..")
+
+    # Rest of mini stacks in incremental-mode
+    comp_slc_files: list[Path] = []
+    slc_idx_end = ministack_size + 1
+    while slc_idx_end <= num_dates:
+        cur_path = Path(f"stack_{slc_idx_start}_{slc_idx_end}")
+        cur_path.mkdir(exist_ok=True)
+        os.chdir(cur_path)
+
+        logger.info(f"***** START: {cur_path} *****")
+        # Get the nearest amplitude mean/dispersion files
+        cur_slc_files = _get_all_slc_files(
+            burst_to_vrt_stack, slc_idx_start, slc_idx_end
+        )
+        cfg = _create_cfg(
+            slc_files=comp_slc_files + cur_slc_files,
+            amplitude_mean_files=all_amp_files,
+            amplitude_dispersion_files=all_disp_files,
+            **arg_dict,
+        )
+        cfg.to_yaml("dolphin_config.yaml")
+        s1_disp.run(cfg)
+
+        # On the step before we hit double `ministack_size`,
+        # archive, shrink, and pull another compressed SLC to replace.
+        stack_size = slc_idx_end - slc_idx_start
+        if stack_size == max_stack_size:
+            # time to shrink!
+            # Get the compressed SLC that was output
+            comp_slc_path = Path("output/compressed_slcs/").resolve()
+            new_comp_slcs = list(comp_slc_path.glob("*.h5"))
+            # assert len(new_comp_slcs) == 1
+            # comp_slc_files.append(new_comp_slcs[0].resolve())
+            comp_slc_files.extend(new_comp_slcs)
+
+            # Move the front forward one ministack
+            slc_idx_start += ministack_size
+
+        slc_idx_end += 1
+        logger.info(f"***** END: {cur_path} *****")
+
+        os.chdir("..")
+
+
+if __name__ == "__main__":
+    args = get_cli_args()
+    arg_dict = vars(args)
+    if arg_dict.pop("pre_compute"):
+        precompute_ps_files(arg_dict)
+    else:
+        main(arg_dict)
diff --git a/src/dolphin/_types.py b/src/dolphin/_types.py
index c9156cef..1d231d4d 100644
--- a/src/dolphin/_types.py
+++ b/src/dolphin/_types.py
@@ -1,5 +1,7 @@
+from __future__ import annotations
+
 from os import PathLike
-from typing import TYPE_CHECKING, Union
+from typing import TYPE_CHECKING, Tuple, Union
 
 # Some classes are declared as generic in stubs, but not at runtime.
 # In Python 3.9 and earlier, os.PathLike is not subscriptable, results in a runtime error
@@ -10,3 +12,6 @@
     PathLikeStr = PathLike
 
 Filename = Union[str, PathLikeStr]
+
+# left, bottom, right, top
+Bbox = Tuple[float, float, float, float]
diff --git a/src/dolphin/io.py b/src/dolphin/io.py
index ca40e1a8..081fe251 100644
--- a/src/dolphin/io.py
+++ b/src/dolphin/io.py
@@ -18,7 +18,7 @@
 
 from dolphin._background import _DEFAULT_TIMEOUT, BackgroundReader, BackgroundWriter
 from dolphin._log import get_log
-from dolphin._types import Filename
+from dolphin._types import Bbox, Filename
 from dolphin.utils import gdal_to_numpy_type, numpy_to_gdal_type, progress
 
 gdal.UseExceptions()
@@ -56,7 +56,7 @@ def load_gdal(
     filename: Filename,
     *,
     band: Optional[int] = None,
-    subsample_factor: int = 1,
+    subsample_factor: Union[int, tuple[int, int]] = 1,
     rows: Optional[slice] = None,
     cols: Optional[slice] = None,
     masked: bool = False,
@@ -69,7 +69,7 @@ def load_gdal(
         Path to the file to load.
     band : int, optional
         Band to load. If None, load all bands as 3D array.
-    subsample_factor : int, optional
+    subsample_factor : int or tuple[int, int], optional
         Subsample the data by this factor. Default is 1 (no subsampling).
         Uses nearest neighbor resampling.
     rows : slice, optional
@@ -91,6 +91,9 @@ def load_gdal(
     # Make an output object of the right size
     dt = gdal_to_numpy_type(ds.GetRasterBand(1).DataType)
 
+    if isinstance(subsample_factor, int):
+        subsample_factor = (subsample_factor, subsample_factor)
+
     if rows is not None and cols is not None:
         xoff, yoff = cols.start, rows.start
         row_stop = min(rows.stop, nrows)
@@ -101,11 +104,17 @@ def load_gdal(
                 f"Invalid row/col slices: {rows}, {cols} for file {filename} of size"
                 f" {nrows}x{ncols}"
             )
-        nrows_out, ncols_out = ysize // subsample_factor, xsize // subsample_factor
+        nrows_out, ncols_out = (
+            ysize // subsample_factor[0],
+            xsize // subsample_factor[1],
+        )
     else:
         xoff, yoff = 0, 0
         xsize, ysize = ncols, nrows
-        nrows_out, ncols_out = nrows // subsample_factor, ncols // subsample_factor
+        nrows_out, ncols_out = (
+            nrows // subsample_factor[0],
+            ncols // subsample_factor[1],
+        )
     # Read the data, and decimate if specified
     resamp = gdal.GRA_NearestNeighbour
     if band is None:
@@ -246,8 +255,8 @@ def get_raster_gt(filename: Filename) -> list[float]:
 
     Returns
     -------
-    tuple[float, float, float, float, float, float]
-        Geotransform.
+    List[float]
+        6 floats representing a GDAL Geotransform.
     """
     ds = gdal.Open(fspath(filename))
     gt = ds.GetGeoTransform()
@@ -292,7 +301,7 @@ def get_raster_driver(filename: Filename) -> str:
 
 def get_raster_bounds(
     filename: Optional[Filename] = None, ds: Optional[gdal.Dataset] = None
-) -> tuple[float, float, float, float]:
+) -> Bbox:
     """Get the (left, bottom, right, top) bounds of the image."""
     if ds is None:
         if filename is None:
@@ -638,6 +647,7 @@ def __init__(
         nodata_mask: Optional[ArrayLike] = None,
         queue_size: int = 1,
         timeout: float = _DEFAULT_TIMEOUT,
+        show_progress: bool = True,
     ):
         super().__init__(nq=queue_size, timeout=timeout, name="EagerLoader")
         self.filename = filename
@@ -656,6 +666,7 @@ def __init__(
         self._nodata_mask = nodata_mask
         self._block_shape = block_shape
         self._nodata = get_raster_nodata(filename)
+        self._show_progress = show_progress
         if self._nodata is None:
             self._nodata = np.nan
 
@@ -681,7 +692,7 @@ def iter_blocks(
 
         s_iter = range(len(queued_slices))
         desc = f"Processing {self._block_shape} sized blocks..."
-        with progress() as p:
+        with progress(dummy=not self._show_progress) as p:
             for _ in p.track(s_iter, description=desc):
                 cur_block, (rows, cols) = self.get_data()
                 logger.debug(f"got data for {rows, cols}: {cur_block.shape}")
diff --git a/src/dolphin/phase_link/mle.py b/src/dolphin/phase_link/mle.py
index ecd2eb14..048e8498 100644
--- a/src/dolphin/phase_link/mle.py
+++ b/src/dolphin/phase_link/mle.py
@@ -150,7 +150,9 @@ def run_mle(
 
     # Fill in the PS pixels from the original SLC stack, if it was given
     if np.any(ps_mask):
-        _fill_ps_pixels(mle_est, temp_coh, slc_stack, ps_mask, strides, avg_mag)
+        _fill_ps_pixels(
+            mle_est, temp_coh, slc_stack, ps_mask, strides, avg_mag, reference_idx
+        )
 
     return mle_est, temp_coh
 
@@ -286,7 +288,14 @@ def _check_all_nans(slc_stack):
 
 
 def _fill_ps_pixels(
-    mle_est, temp_coh, slc_stack, ps_mask, strides, avg_mag, use_max_ps: bool = False
+    mle_est: np.ndarray,
+    temp_coh: np.ndarray,
+    slc_stack: np.ndarray,
+    ps_mask: np.ndarray,
+    strides: dict[str, int],
+    avg_mag: np.ndarray,
+    reference_idx: int = 0,
+    use_max_ps: bool = False,
 ):
     """Fill in the PS locations in the MLE estimate with the original SLC data.
 
@@ -307,6 +316,9 @@ def _fill_ps_pixels(
     avg_mag : np.ndarray, optional
         The average magnitude of the SLC stack, used to to find the brightest
         PS pixels to fill within each look window.
+    reference_idx : int, default = 0
+        SLC to use as reference for PS pixels. All pixel values are multiplied
+        by the conjugate of
     use_max_ps : bool, optional
         If True, use the brightest PS pixel in each look window to fill in the
         MLE estimate. If False, use the average of all PS pixels in each look window.
@@ -334,7 +346,7 @@ def _fill_ps_pixels(
         # Get the indices of the brightest pixels within each look window
         slc_r_idxs, slc_c_idxs = _get_max_idxs(mag, strides["y"], strides["x"])
         # we're only filling where there are PS pixels
-        ref = np.conj(slc_stack[0][slc_r_idxs, slc_c_idxs])
+        ref = np.exp(-1j * np.angle(slc_stack[reference_idx][slc_r_idxs, slc_c_idxs]))
         for i in range(len(slc_stack)):
             mle_est[i][ps_mask_looked] = slc_stack[i][slc_r_idxs, slc_c_idxs] * ref
     else:
@@ -356,7 +368,9 @@ def _get_avg_ps(
     slc_stack_nanned = slc_stack.copy()
     slc_stack_nanned[:, ~ps_mask] = np.nan
     # Reference all ps pixels in the SLC stack to the first SLC
-    slc_stack_nanned[:, ps_mask] *= np.conj(slc_stack_nanned[0, ps_mask])[None]
+    slc_stack_nanned[:, ps_mask] *= np.exp(
+        -1j * np.angle(slc_stack_nanned[0, ps_mask])
+    )[None]
     # Then, take the average of all PS pixels within each look window
     return take_looks(
         slc_stack_nanned,
diff --git a/src/dolphin/phase_link/simulate.py b/src/dolphin/phase_link/simulate.py
index 7b754a9d..7d078043 100644
--- a/src/dolphin/phase_link/simulate.py
+++ b/src/dolphin/phase_link/simulate.py
@@ -105,7 +105,6 @@ def simulate_C(
 
     C = _sim_coherence_mat(t, gamma0, gamma_inf, Tau0, signal_phase)
     return C, truth
-    # return C
 
 
 @njit(cache=True)
@@ -142,11 +141,11 @@ def _sim_signal(
     # we divided std by 2 since we're subtracting the first value
     signal_phase = signal_phase - signal_phase[0]
 
-    # wrap the phase to -pi to p
-    signal_phase = np.angle(np.exp(1j * signal_phase))
-    truth = np.angle(np.exp(1j * (truth - truth[0])))
+    # wrap the signal_phase to -pi to pi
+    signal_phase = np.mod(signal_phase + np.pi, 2 * np.pi) - np.pi
+    truth = np.mod((truth - truth[0]) + np.pi, 2 * np.pi) - np.pi
 
-    return signal_phase, truth
+    return signal_phase.astype(np.float64), truth.astype(np.float64)
 
 
 def rmse(x, y):
diff --git a/src/dolphin/ps.py b/src/dolphin/ps.py
index 691e548e..ef4d0c77 100644
--- a/src/dolphin/ps.py
+++ b/src/dolphin/ps.py
@@ -31,8 +31,10 @@ def create_ps(
     amp_dispersion_threshold: float = 0.25,
     existing_amp_mean_file: Optional[Filename] = None,
     existing_amp_dispersion_file: Optional[Filename] = None,
+    nodata_mask: Optional[np.ndarray] = None,
     update_existing: bool = False,
     block_size_gb: float = 1.0,
+    show_progress: bool = True,
 ):
     """Create the amplitude dispersion, mean, and PS files.
 
@@ -52,6 +54,9 @@ def create_ps(
         An existing amplitude mean file to use, by default None.
     existing_amp_dispersion_file : Optional[Filename], optional
         An existing amplitude dispersion file to use, by default None.
+    nodata_mask : Optional[np.ndarray]
+        If provided, skips computing PS over areas where the mask is False
+        Otherwise, loads input data from everywhere and calculates.
     update_existing : bool, optional
         If providing existing amp mean/dispersion files, combine them with the
         data from the current SLC stack.
@@ -60,6 +65,8 @@ def create_ps(
     block_size_gb : int, optional
         The maximum amount of data to read at a time (in GB).
         Default is 1.0 GB.
+    show_progress : bool, default=True
+        If true, displays a `rich` ProgressBar.
     """
     if existing_amp_dispersion_file and existing_amp_mean_file and not update_existing:
         logger.info("Using existing amplitude dispersion file, skipping calculation.")
@@ -94,18 +101,27 @@ def create_ps(
     # Initialize the intermediate arrays for the calculation
     magnitude = np.zeros((len(vrt_stack), *block_shape), dtype=np.float32)
 
+    skip_empty = nodata_mask is None
+
     writer = io.Writer()
     # Make the generator for the blocks
     block_gen = vrt_stack.iter_blocks(
         max_bytes=max_bytes,
-        skip_empty=False,
+        skip_empty=skip_empty,
+        nodata_mask=nodata_mask,
+        show_progress=show_progress,
     )
     for cur_data, (rows, cols) in block_gen:
         cur_rows, cur_cols = cur_data.shape[-2:]
 
         if not (np.all(cur_data == 0) or np.all(np.isnan(cur_data))):
             magnitude_cur = np.abs(cur_data, out=magnitude[:, :cur_rows, :cur_cols])
-            mean, amp_disp, ps = calc_ps_block(magnitude_cur, amp_dispersion_threshold)
+            mean, amp_disp, ps = calc_ps_block(
+                # use min_count == size of stack so that ALL need to be not Nan
+                magnitude_cur,
+                amp_dispersion_threshold,
+                min_count=len(magnitude_cur),
+            )
 
             # Use the UInt8 type for the PS to save.
             # For invalid pixels, set to max Byte value
diff --git a/src/dolphin/stack.py b/src/dolphin/stack.py
index 9e775835..e59849bb 100755
--- a/src/dolphin/stack.py
+++ b/src/dolphin/stack.py
@@ -11,7 +11,7 @@
 
 from dolphin import io, utils
 from dolphin._log import get_log
-from dolphin._types import Filename
+from dolphin._types import Bbox, Filename
 
 gdal.UseExceptions()
 logger = get_log(__name__)
@@ -57,6 +57,9 @@ class VRTStack:
     file_date_fmt : str, optional (default = "%Y%m%d")
         Format string for parsing the dates from the filenames.
         Passed to [dolphin.utils.get_dates][].
+    num_read_threads : int, optional
+        If provided and `num_read_threads > 1`, uses `threading` to read stack layers
+        in parallel. Default = 1, no threading.
     """
 
     def __init__(
@@ -67,7 +70,7 @@ def __init__(
         subdataset: Optional[str] = None,
         pixel_bbox: Optional[tuple[int, int, int, int]] = None,
         target_extent: Optional[tuple[float, float, float, float]] = None,
-        latlon_bbox: Optional[tuple[float, float, float, float]] = None,
+        latlon_bbox: Optional[Bbox] = None,
         sort_files: bool = True,
         nodata_value: Optional[float] = None,
         file_date_fmt: str = "%Y%m%d",
@@ -383,6 +386,7 @@ def iter_blocks(
         max_bytes: Optional[float] = DEFAULT_BLOCK_BYTES,
         skip_empty: bool = True,
         nodata_mask: Optional[np.ndarray] = None,
+        show_progress: bool = True,
     ) -> Generator[tuple[np.ndarray, tuple[slice, slice]], None, None]:
         """Iterate over blocks of the stack.
 
@@ -406,6 +410,8 @@ def iter_blocks(
             Optional mask indicating nodata values. If provided, will skip
             blocks that are entirely nodata.
             1s are the nodata values, 0s are valid data.
+        show_progress : bool, default=True
+            If true, displays a `rich` ProgressBar.
 
         Yields
         ------
@@ -422,6 +428,7 @@ def iter_blocks(
             overlaps=overlaps,
             nodata_mask=nodata_mask,
             skip_empty=skip_empty,
+            show_progress=show_progress,
         )
         yield from self._loader.iter_blocks()
 
@@ -465,6 +472,8 @@ def __getitem__(self, index):
                 index = len(self) + index
             return self.read_stack(band=index + 1)
 
+        # TODO: raise an error if they try to skip like [::2, ::2]
+        # or pass it to read_stack... but I dont think I need to support it.
         n, rows, cols = index
         if isinstance(rows, int):
             rows = slice(rows, rows + 1)
@@ -476,12 +485,14 @@ def __getitem__(self, index):
             return self.read_stack(band=n + 1, rows=rows, cols=cols)
 
         bands = list(range(1, 1 + len(self)))[n]
-        data = [self.read_stack(band=i, rows=rows, cols=cols) for i in bands]
-        # TODO: should multithread this read?
-        # data = Parallel(n_jobs=min(len(bands), 10), prefer="threads")(
-        #     delayed(self.read_stack)(band=i, rows=rows, cols=cols) for i in bands
-        # )
-        return np.stack(data, axis=0).squeeze()
+        if len(bands) == len(self):
+            # This will use gdal's ds.ReadAsRaster, no iteration needed
+            data = self.read_stack(band=None, rows=rows, cols=cols)
+        else:
+            data = np.stack(
+                [self.read_stack(band=i, rows=rows, cols=cols) for i in bands], axis=0
+            )
+        return data.squeeze()
 
     @property
     def dtype(self):
diff --git a/src/dolphin/stitching.py b/src/dolphin/stitching.py
index a13b0c3c..86f8fe87 100644
--- a/src/dolphin/stitching.py
+++ b/src/dolphin/stitching.py
@@ -1,7 +1,6 @@
 """stitching.py: utilities for combining interferograms into larger images."""
 from __future__ import annotations
 
-import itertools
 import math
 import tempfile
 from datetime import date
@@ -12,10 +11,11 @@
 import numpy as np
 from numpy.typing import DTypeLike
 from osgeo import gdal, osr
+from pyproj import Transformer
 
 from dolphin import io, utils
 from dolphin._log import get_log
-from dolphin._types import Filename
+from dolphin._types import Bbox, Filename
 
 logger = get_log(__name__)
 
@@ -56,7 +56,7 @@ def merge_by_date(
     This function is intended to be used with filenames that contain date pairs
     (from interferograms).
     """
-    grouped_images = _group_by_date(image_file_list, file_date_fmt=file_date_fmt)
+    grouped_images = utils.group_by_date(image_file_list, file_date_fmt=file_date_fmt)
     stitched_acq_times = {}
 
     for dates, cur_images in grouped_images.items():
@@ -86,6 +86,8 @@ def merge_images(
     file_list: Sequence[Filename],
     outfile: Filename,
     target_aligned_pixels: bool = True,
+    out_bounds: Optional[Bbox] = None,
+    out_bounds_epsg: Optional[int] = None,
     driver: str = "ENVI",
     out_nodata: Optional[float] = 0,
     out_dtype: Optional[DTypeLike] = None,
@@ -106,6 +108,13 @@ def merge_images(
         are integer multiples of pixel size, matching the ``-tap``
         options of GDAL utilities.
         Default is True.
+    out_bounds: Optional[tuple[float]]
+        if provided, forces the output image bounds to
+            (left, bottom, right, top).
+        Otherwise, computes from the outside of all input images.
+    out_bounds_epsg: Optional[int]
+        EPSG code for the `out_bounds`.
+        If not provided, assumed to match the projections of `file_list`.
     driver : str
         GDAL driver to use for output file. Default is ENVI.
     out_nodata : Optional[float]
@@ -146,7 +155,7 @@ def merge_images(
     temp_dir = tempfile.TemporaryDirectory()
     # temp_dir = Path("test_ifgs")
 
-    warped_file_list = _warp_to_projection(
+    warped_file_list = warp_to_projection(
         file_list,
         # temp_dir,
         Path(temp_dir.name),
@@ -155,8 +164,11 @@ def merge_images(
     )
     # Compute output array shape. We guarantee it will cover the output
     # bounds completely
-    bounds, gt = _get_combined_bounds_gt(  # type: ignore
-        *warped_file_list, target_aligned_pixels=target_aligned_pixels
+    bounds, gt = get_combined_bounds_gt(  # type: ignore
+        *warped_file_list,
+        target_aligned_pixels=target_aligned_pixels,
+        out_bounds=out_bounds,
+        out_bounds_epsg=out_bounds_epsg,
     )
 
     out_shape = _get_output_shape(bounds, res)
@@ -236,47 +248,6 @@ def merge_images(
     temp_dir.cleanup()
 
 
-def _group_by_date(
-    file_list: Iterable[Filename], file_date_fmt: str = io.DEFAULT_DATETIME_FORMAT
-):
-    """Combine Sentinel objects by date.
-
-    Parameters
-    ----------
-    file_list: Iterable[Filename]
-        path to folder containing CSLC files
-    file_date_fmt: str
-        format of the date in the filename.
-        Default is [dolphin.io.DEFAULT_DATETIME_FORMAT][]
-
-    Returns
-    -------
-    dict
-        key is the date of the SLC acquisition
-        Value is a list of Paths on that date:
-        [(datetime.date(2017, 10, 13),
-          [Path(...)
-            Path(...),
-            ...]),
-         (datetime.date(2017, 10, 25),
-          [Path(...)
-            Path(...),
-            ...]),
-    """
-    sorted_file_list, _ = utils.sort_files_by_date(
-        file_list, file_date_fmt=file_date_fmt
-    )
-
-    # Now collapse into groups, sorted by the date
-    grouped_images = {
-        dates: list(g)
-        for dates, g in itertools.groupby(
-            sorted_file_list, key=lambda x: tuple(utils.get_dates(x))
-        )
-    }
-    return grouped_images
-
-
 def _blend_new_arr(
     cur_arr: np.ndarray, new_arr: np.ndarray, nodata_vals: Iterable[Optional[float]]
 ):
@@ -314,8 +285,8 @@ def _blend_new_arr(
     return cur_arr
 
 
-def _warp_to_projection(
-    filenames: Iterable[Filename],
+def warp_to_projection(
+    filenames: Sequence[Filename],
     dirname: Filename,
     projection: str,
     res: tuple[float, float],
@@ -328,7 +299,7 @@ def _warp_to_projection(
 
     Parameters
     ----------
-    filenames : Iterable[Filename]
+    filenames : Sequence[Filename]
         list of filenames to warp.
     dirname : Filename
         The directory to write the warped files to.
@@ -342,6 +313,11 @@ def _warp_to_projection(
     list[Filename]
         The warped filenames.
     """
+    if projection is None:
+        projection = _get_mode_projection(filenames)
+    if res is None:
+        res = _get_resolution(filenames)
+
     warped_files = []
     for idx, fn in enumerate(filenames):
         fn = Path(fn)
@@ -387,10 +363,12 @@ def _get_resolution(filenames: Iterable[Filename]) -> tuple[float, float]:
     return res[0]
 
 
-def _get_combined_bounds_gt(
+def get_combined_bounds_gt(
     *filenames: Filename,
     target_aligned_pixels: bool = False,
-) -> tuple[tuple[float, float, float, float], list]:
+    out_bounds: Optional[Bbox] = None,
+    out_bounds_epsg: Optional[int] = None,
+) -> tuple[Bbox, list[float]]:
     """Get the bounds and geotransform of the combined image.
 
     Parameters
@@ -399,11 +377,18 @@ def _get_combined_bounds_gt(
         list of filenames to combine
     target_aligned_pixels : bool
         if True, adjust output image bounds so that pixel coordinates
-        are integer multiples of pixel size, matching the ``-tap``.
+        are integer multiples of pixel size, matching the `-tap` GDAL option.
+    out_bounds: Optional[Bbox]
+        if provided, forces the output image bounds to
+            (left, bottom, right, top).
+        Otherwise, computes from the outside of all input images.
+    out_bounds_epsg: Optional[int]
+        The EPSG of `out_bounds`. If not provided, assumed to be the same
+        as the EPSG of all `filenames`.
 
     Returns
     -------
-    bounds : tuple[float]
+    bounds : Bbox
         (min_x, min_y, max_x, max_y)
     gt : list[float]
         geotransform of the combined image.
@@ -413,6 +398,8 @@ def _get_combined_bounds_gt(
     ys = []
     resolutions = set()
     projs = set()
+
+    # Check all files match in resolution/projection
     for fn in filenames:
         ds = gdal.Open(fspath(fn))
         left, bottom, right, top = io.get_raster_bounds(fn)
@@ -429,9 +416,17 @@ def _get_combined_bounds_gt(
         raise ValueError(f"The input files have different resolutions: {resolutions}. ")
     if len(projs) > 1:
         raise ValueError(f"The input files have different projections: {projs}. ")
-
     res = (abs(dx), abs(dy))
-    bounds = min(xs), min(ys), max(xs), max(ys)
+
+    if out_bounds is not None:
+        if out_bounds_epsg is not None:
+            dst_epsg = io.get_raster_crs(filenames[0]).to_epsg()
+            bounds = _reproject_bounds(out_bounds, out_bounds_epsg, dst_epsg)
+        else:
+            bounds = out_bounds  # type: ignore
+    else:
+        bounds = min(xs), min(ys), max(xs), max(ys)
+
     if target_aligned_pixels:
         bounds = _align_bounds(bounds, res)
 
@@ -457,6 +452,13 @@ def _align_bounds(bounds: Iterable[float], res: tuple[float, float]):
     return (left, bottom, right, top)
 
 
+def _reproject_bounds(bounds: Bbox, src_epsg: int, dst_epsg: int) -> Bbox:
+    t = Transformer.from_crs(src_epsg, dst_epsg, always_xy=True)
+    left, bottom, right, top = bounds
+    bbox: Bbox = (*t.transform(left, bottom), *t.transform(right, top))  # type: ignore
+    return bbox
+
+
 def _nodata_to_zero(
     infile: Filename,
     outfile: Optional[Filename] = None,
diff --git a/src/dolphin/utils.py b/src/dolphin/utils.py
index 1bee4fa0..1e7062dc 100644
--- a/src/dolphin/utils.py
+++ b/src/dolphin/utils.py
@@ -1,11 +1,13 @@
 from __future__ import annotations
 
 import datetime
+import itertools
 import math
 import re
 import resource
 import sys
 import warnings
+from contextlib import nullcontext
 from multiprocessing import cpu_count
 from pathlib import Path
 from typing import Iterable, Optional, Union
@@ -23,9 +25,14 @@
 logger = get_log(__name__)
 
 
-def progress():
+def progress(dummy=False):
     """Create a Progress bar context manager.
 
+    Parameters
+    ----------
+    dummy : bool, default = False
+        If True, skips showing and calls `contextlib.nullcontext`
+
     Usage
     -----
     >>> with progress() as p:
@@ -33,6 +40,9 @@ def progress():
     ...         pass
     10/10 Working... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00
     """
+    if dummy:
+        return nullcontext()
+
     return Progress(
         SpinnerColumn(),
         MofNCompleteColumn(),
@@ -230,6 +240,51 @@ def sort_key(file_date_tuple):
     return list(file_list), list(dates)
 
 
+def group_by_date(
+    file_list: Iterable[Filename], file_date_fmt: Optional[str] = None
+) -> dict[tuple[datetime.date, ...], list[Filename]]:
+    """Combine files by date into a dict.
+
+    Parameters
+    ----------
+    file_list: Iterable[Filename]
+        Path to folder containing files with dates in the filename.
+    file_date_fmt: str
+        Format of the date in the filename.
+        Default is [dolphin.io.DEFAULT_DATETIME_FORMAT][]
+
+    Returns
+    -------
+    dict
+        key is a list of dates in the filenames.
+        Value is a list of Paths on that date.
+        E.g.:
+        {(datetime.date(2017, 10, 13),
+          [Path(...)
+            Path(...),
+            ...]),
+         (datetime.date(2017, 10, 25),
+          [Path(...)
+            Path(...),
+            ...]),
+        }
+    """
+    if file_date_fmt is None:
+        import dolphin.io
+
+        file_date_fmt = dolphin.io.DEFAULT_DATETIME_FORMAT
+    sorted_file_list, _ = sort_files_by_date(file_list, file_date_fmt=file_date_fmt)
+
+    # Now collapse into groups, sorted by the date
+    grouped_images = {
+        dates: list(g)
+        for dates, g in itertools.groupby(
+            sorted_file_list, key=lambda x: tuple(get_dates(x))
+        )
+    }
+    return grouped_images
+
+
 @njit
 def _get_slices(half_r: int, half_c: int, r: int, c: int, rows: int, cols: int):
     """Get the slices for the given pixel and half window size."""
diff --git a/src/dolphin/workflows/_product.py b/src/dolphin/workflows/_product.py
index df1785e0..73cae207 100644
--- a/src/dolphin/workflows/_product.py
+++ b/src/dolphin/workflows/_product.py
@@ -395,11 +395,15 @@ def form_name(filename: Path, burst: str):
     group_name = "/".join(parts)
 
     for burst, comp_slc_file in comp_slc_dict.items():
+        outname = Path(output_dir) / form_name(comp_slc_file, burst)
+        if outname.exists():
+            logger.info(f"Skipping existing {outname}")
+            continue
+
         crs = io.get_raster_crs(comp_slc_file)
         gt = io.get_raster_gt(comp_slc_file)
         data = _zero_mantissa(io.load_gdal(comp_slc_file))
 
-        outname = Path(output_dir) / form_name(comp_slc_file, burst)
         logger.info(f"Writing {outname}")
         with h5py.File(outname, "w") as hf:
             # add type to root for GDAL recognition of complex datasets in NetCDF
diff --git a/src/dolphin/workflows/_utils.py b/src/dolphin/workflows/_utils.py
index 838a910d..9b8497eb 100644
--- a/src/dolphin/workflows/_utils.py
+++ b/src/dolphin/workflows/_utils.py
@@ -71,6 +71,9 @@ def sort_by_burst_id(file_list):
         file_list, burst_ids = zip(*file_burst_tups)  # type: ignore
         return file_list
 
+    if not file_list:
+        return {}
+
     sorted_file_list = sort_by_burst_id(file_list)
     # Now collapse into groups, sorted by the burst_id
     grouped_images = {
@@ -229,6 +232,13 @@ def make_nodata_mask(
             logger.info(f"Overwriting {out_file} since overwrite=True.")
             Path(out_file).unlink()
 
+    # Check these are the right format to get nodata polygons
+    try:
+        test_f = f"NETCDF:{opera_file_list[0]}:{OPERA_DATASET_NAME}"
+        io.get_raster_xysize(test_f)
+    except RuntimeError:
+        raise ValueError(f"Unable to open {test_f}")
+
     # convert pixels to degrees lat/lon
     # TODO: more robust way to get the pixel size... this is a hack
     # maybe just use pyproj to warp lat/lon to meters and back?
diff --git a/src/dolphin/workflows/s1_disp.py b/src/dolphin/workflows/s1_disp.py
index b5f860c2..9eae2175 100755
--- a/src/dolphin/workflows/s1_disp.py
+++ b/src/dolphin/workflows/s1_disp.py
@@ -84,7 +84,9 @@ def run(
         _remove_dir_if_empty(cfg.ps_options._directory)
 
     else:
-        wrapped_phase_cfgs = [("", cfg)]
+        # grab the only key (either a burst, or "") and use that
+        b = list(grouped_slc_files.keys())[0]
+        wrapped_phase_cfgs = [(b, cfg)]
     # ###########################
     # 1. Wrapped phase estimation
     # ###########################
diff --git a/src/dolphin/workflows/sequential.py b/src/dolphin/workflows/sequential.py
index 421db82d..29532c59 100644
--- a/src/dolphin/workflows/sequential.py
+++ b/src/dolphin/workflows/sequential.py
@@ -102,7 +102,8 @@ def run_wrapped_phase_sequential(
             output_folder=cur_output_folder,
             half_window=half_window,
             strides=strides,
-            reference_idx=mini_idx,
+            # TODO: what situations do we need to set reference-idx != 0
+            reference_idx=0,
             beta=beta,
             mask_file=mask_file,
             ps_mask_file=ps_mask_file,
diff --git a/src/dolphin/workflows/single.py b/src/dolphin/workflows/single.py
index 9fd53477..e0970308 100644
--- a/src/dolphin/workflows/single.py
+++ b/src/dolphin/workflows/single.py
@@ -78,8 +78,9 @@ def run_wrapped_phase_single(
         ps_mask = np.zeros_like(nodata_mask)
 
     if amp_mean_file is not None and amp_dispersion_file is not None:
-        amp_mean = io.load_gdal(amp_mean_file, masked=True)
-        amp_dispersion = io.load_gdal(amp_dispersion_file, masked=True)
+        # Note: have to fill, since numba (as of 0.57) can't do masked arrays
+        amp_mean = io.load_gdal(amp_mean_file, masked=True).filled(np.nan)
+        amp_dispersion = io.load_gdal(amp_dispersion_file, masked=True).filled(np.nan)
         # convert back to variance from dispersion: amp_disp = std_dev / mean
         amp_variance = (amp_dispersion * amp_mean) ** 2
     else:
@@ -197,9 +198,6 @@ def run_wrapped_phase_single(
             amp_stack=amp_stack,
             method=shp_method,
         )
-        # Save the count for each pixel
-        shp_counts = np.sum(neighbor_arrays, axis=(-2, -1))
-
         # Run the phase linking process on the current ministack
         try:
             cur_mle_stack, tcorr = run_mle(
@@ -243,7 +241,8 @@ def run_wrapped_phase_single(
         # Save the temporal coherence blocks
         writer.queue_write(tcorr, tcorr_file, out_row_start, out_col_start)
 
-        # Save the SHP counts (if not using Rect window)
+        # Save the SHP counts for each pixel (if not using Rect window)
+        shp_counts = np.sum(neighbor_arrays[rows, cols], axis=(-2, -1))
         writer.queue_write(shp_counts, shp_counts_file, out_row_start, out_col_start)
 
         # Compress the ministack using only the non-compressed SLCs
diff --git a/src/dolphin/workflows/wrapped_phase.py b/src/dolphin/workflows/wrapped_phase.py
index a42d1f46..ec57f4a3 100644
--- a/src/dolphin/workflows/wrapped_phase.py
+++ b/src/dolphin/workflows/wrapped_phase.py
@@ -45,15 +45,11 @@ def run(cfg: Workflow, debug: bool = False) -> tuple[list[Path], Path, Path]:
     # Make a VRT pointing to the input SLC files
     # #############################################
     subdataset = cfg.input_options.subdataset
-    vrt_path = cfg.scratch_directory / "slc_stack.vrt"
-    if vrt_path.exists():
-        vrt_stack = stack.VRTStack.from_vrt_file(vrt_path)
-    else:
-        vrt_stack = stack.VRTStack(
-            input_file_list,
-            subdataset=subdataset,
-            outfile=cfg.scratch_directory / "slc_stack.vrt",
-        )
+    vrt_stack = stack.VRTStack(
+        input_file_list,
+        subdataset=subdataset,
+        outfile=cfg.scratch_directory / "slc_stack.vrt",
+    )
 
     # Make the nodata mask from the polygons, if we're using OPERA CSLCs
     try:
@@ -187,11 +183,12 @@ def run(cfg: Workflow, debug: bool = False) -> tuple[list[Path], Path, Path]:
             ref_idx, sec_idx = idxs[0]
             file1, file2 = vrt_stack.file_list[ref_idx], vrt_stack.file_list[sec_idx]
             date1, date2 = utils.get_dates(file1)[0], utils.get_dates(file2)[0]
-            suffix = utils.full_suffix(file2)
+            # We're just copying, so get the extension of the file to copy
+            to_copy = phase_linked_slcs[sec_idx]
+            suffix = utils.full_suffix(to_copy)
             ifg_name = ifg_dir / (io._format_date_pair(date1, date2) + suffix)
-            # We're just copying
-            shutil.copyfile(phase_linked_slcs[sec_idx], ifg_name)
-            ifg_file_list = [ifg_name]
+            shutil.copyfile(to_copy, ifg_name)
+            ifg_file_list = [ifg_name]  # return just the one as the "network"
         else:
             network = Network(
                 slc_list=phase_linked_slcs,


Test patch:
diff --git a/tests/test_stitching.py b/tests/test_stitching.py
index 09ec5c17..020ee501 100644
--- a/tests/test_stitching.py
+++ b/tests/test_stitching.py
@@ -38,7 +38,7 @@ def test_get_combined_bounds_gt(shifted_slc_files):
     assert bnds == expected_bnds
 
     # should be same as test #1 above
-    bnds, gt = stitching._get_combined_bounds_gt(
+    bnds, gt = stitching.get_combined_bounds_gt(
         shifted_slc_files[0],
     )
     expected_bnds = (-5.5, -4.5, 4.5, 5.5)
@@ -47,7 +47,7 @@ def test_get_combined_bounds_gt(shifted_slc_files):
     assert gt == expected_gt
 
     # check same file twice
-    bnds, gt = stitching._get_combined_bounds_gt(
+    bnds, gt = stitching.get_combined_bounds_gt(
         shifted_slc_files[0], shifted_slc_files[0]
     )
     assert bnds == expected_bnds
@@ -55,7 +55,7 @@ def test_get_combined_bounds_gt(shifted_slc_files):
 
     # Now combined one: should have the mins as the first two values,
     # and the last two values should be the -1 file
-    bnds, gt = stitching._get_combined_bounds_gt(*shifted_slc_files)
+    bnds, gt = stitching.get_combined_bounds_gt(*shifted_slc_files)
     expected_bnds = (-5.5, -4.5, 8.5, 9.5)
     # only the top left corner should change
     expected_gt[0], expected_gt[3] = -5.5, 9.5
@@ -66,14 +66,12 @@ def test_get_combined_bounds_gt(shifted_slc_files):
 def test_get_combined_bounds_gt_different_proj(
     slc_file_list_nc, slc_file_list_nc_wgs84
 ):
-    bnds, gt = stitching._get_combined_bounds_gt(*slc_file_list_nc)
+    bnds, gt = stitching.get_combined_bounds_gt(*slc_file_list_nc)
     assert bnds == (-5.5, -2.0, 4.5, 3.0)
     assert gt == [-5.5, 1.0, 0, 3.0, 0, -1.0]
 
     with pytest.raises(ValueError):
-        stitching._get_combined_bounds_gt(
-            slc_file_list_nc_wgs84[0], slc_file_list_nc[0]
-        )
+        stitching.get_combined_bounds_gt(slc_file_list_nc_wgs84[0], slc_file_list_nc[0])
 
 
 def test_get_mode_projection(slc_file_list_nc, slc_file_list_nc_wgs84):