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

Issue:
[Bug]: Compressed SLCs are getting created incorrectly for the single update workflow
Here's a sample correlation matrix using 1 Fringe-created CompSlc and 5 normal SLCs:
![2023-02-17-1676665212891](https://user-images.githubusercontent.com/8291800/219786307-14918d44-d1ab-4289-9519-de4fbcdedb57.jpeg)

Here's the same correlation matrix, testing with dolphin functions (top row)
![2023-02-17-1676664648266](https://user-images.githubusercontent.com/8291800/219786312-eb76011a-28be-4d85-900d-1b5c9c4d97f8.jpeg)
(edit: this was not the result of a workflow, but was from testing each step without referencing to the first date)

The fact that it's lower correlation than all others is a clear bug, though it was basically getting ignored in the full sequential workflow (why it's been undetected until now).



Discussion / hints:

Fix patch:
diff --git a/scripts/create_vrt.py b/scripts/create_vrt.py
index 183449e2..c672e830 100644
--- a/scripts/create_vrt.py
+++ b/scripts/create_vrt.py
@@ -4,7 +4,7 @@
 from dolphin._log import get_log
 from dolphin.stack import VRTStack
 
-logger = get_log()
+logger = get_log(__name__)
 
 
 def get_cli_args():
diff --git a/src/dolphin/_background.py b/src/dolphin/_background.py
index c8c16516..4193e11a 100644
--- a/src/dolphin/_background.py
+++ b/src/dolphin/_background.py
@@ -5,9 +5,9 @@
 
 from dolphin._log import get_log
 
-logger = get_log()
+logger = get_log(__name__)
 
-_DEFAULT_TIMEOUT = 1.0
+_DEFAULT_TIMEOUT = 0.5
 
 
 def is_main_thread_active() -> bool:
diff --git a/src/dolphin/_log.py b/src/dolphin/_log.py
index 4a6ec8b7..509bb198 100644
--- a/src/dolphin/_log.py
+++ b/src/dolphin/_log.py
@@ -6,7 +6,7 @@
 Usage:
 
     from ._log import get_log
-    logger = get_log()
+    logger = get_log(__name__)
 
     logger.info("Something happened")
     logger.warning("Something concerning happened")
@@ -100,7 +100,7 @@ def log_runtime(f: Callable) -> Callable:
     def test_func():
         return 2 + 4
     """
-    logger = get_log()
+    logger = get_log(__name__)
 
     @wraps(f)
     def wrapper(*args, **kwargs):
diff --git a/src/dolphin/interferogram.py b/src/dolphin/interferogram.py
index f517bc18..58c042ba 100644
--- a/src/dolphin/interferogram.py
+++ b/src/dolphin/interferogram.py
@@ -14,7 +14,7 @@
 
 gdal.UseExceptions()
 
-logger = get_log()
+logger = get_log(__name__)
 
 
 class VRTInterferogram(BaseModel):
diff --git a/src/dolphin/io.py b/src/dolphin/io.py
index 2ac0f170..a5fbac6b 100644
--- a/src/dolphin/io.py
+++ b/src/dolphin/io.py
@@ -48,7 +48,7 @@
 )
 DEFAULT_DATETIME_FORMAT = "%Y%m%d"
 
-logger = get_log()
+logger = get_log(__name__)
 
 
 def load_gdal(
diff --git a/src/dolphin/phase_link/__init__.py b/src/dolphin/phase_link/__init__.py
index fd349497..564cfed9 100644
--- a/src/dolphin/phase_link/__init__.py
+++ b/src/dolphin/phase_link/__init__.py
@@ -3,4 +3,5 @@
 Currently implements the eigenvalue-based maximum likelihood (EMI)
 algorithm from (Ansari, 2018).
 """
+from ._compress import compress  # noqa: F401
 from .mle import PhaseLinkRuntimeError, run_mle  # noqa: F401
diff --git a/src/dolphin/phase_link/_compress.py b/src/dolphin/phase_link/_compress.py
new file mode 100644
index 00000000..e6976c5a
--- /dev/null
+++ b/src/dolphin/phase_link/_compress.py
@@ -0,0 +1,40 @@
+import warnings
+
+import numpy as np
+
+from dolphin.utils import upsample_nearest
+
+
+def compress(
+    slc_stack: np.ndarray,
+    mle_estimate: np.ndarray,
+):
+    """Compress the stack of SLC data using the estimated phase.
+
+    Parameters
+    ----------
+    slc_stack : np.array
+        The stack of complex SLC data, shape (nslc, rows, cols)
+    mle_estimate : np.array
+        The estimated phase from [`run_mle`][dolphin.phase_link.mle.run_mle],
+        shape (nslc, rows // strides['y'], cols // strides['x'])
+
+    Returns
+    -------
+    np.array
+        The compressed SLC data, shape (rows, cols)
+    """
+    # If the output is downsampled, we need to make `mle_estimate` the same shape
+    # as the output
+    mle_estimate_upsampled = upsample_nearest(mle_estimate, slc_stack.shape[1:])
+    # For each pixel, project the SLCs onto the (normalized) estimated phase
+    # by performing a pixel-wise complex dot product
+    mle_norm = np.linalg.norm(mle_estimate_upsampled, axis=0)
+    # Avoid divide by zero (there may be 0s at the upsampled boundary)
+    mle_norm[mle_norm == 0] = np.nan
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore", message="invalid value encountered")
+        return (
+            np.nansum(slc_stack * np.conjugate(mle_estimate_upsampled), axis=0)
+            / mle_norm
+        )
diff --git a/src/dolphin/phase_link/mle.py b/src/dolphin/phase_link/mle.py
index 967a593e..81ed5498 100644
--- a/src/dolphin/phase_link/mle.py
+++ b/src/dolphin/phase_link/mle.py
@@ -220,7 +220,7 @@ def mle_stack(C_arrays, beta: float = 0.01, reference_idx: float = 0):
     # Return the phase (still as a GPU array)
     phase_stack = xp.angle(evd_estimate)
     # Move the SLC dimension to the front (to match the SLC stack shape)
-    return np.moveaxis(phase_stack, -1, 0)
+    return xp.moveaxis(phase_stack, -1, 0)
 
 
 def _get_eigvecs(C):
diff --git a/src/dolphin/ps.py b/src/dolphin/ps.py
index 0b2953d3..1ffcd453 100644
--- a/src/dolphin/ps.py
+++ b/src/dolphin/ps.py
@@ -15,7 +15,7 @@
 
 gdal.UseExceptions()
 
-logger = get_log()
+logger = get_log(__name__)
 
 
 def create_ps(
diff --git a/src/dolphin/sequential.py b/src/dolphin/sequential.py
deleted file mode 100644
index eb41ae08..00000000
--- a/src/dolphin/sequential.py
+++ /dev/null
@@ -1,423 +0,0 @@
-"""Estimate wrapped phase using batches of ministacks.
-
-References
-----------
-    [1] Ansari, H., De Zan, F., & Bamler, R. (2017). Sequential estimator: Toward
-    efficient InSAR time series analysis. IEEE Transactions on Geoscience and
-    Remote Sensing, 55(10), 5637-5652.
-"""
-import warnings
-from collections import defaultdict
-from math import nan
-from os import fspath
-from pathlib import Path
-from typing import Dict, List, Optional
-
-import numpy as np
-from osgeo_utils import gdal_calc
-
-from dolphin import io
-from dolphin._log import get_log
-from dolphin._types import Filename
-from dolphin.interferogram import VRTInterferogram
-from dolphin.phase_link import PhaseLinkRuntimeError, run_mle
-from dolphin.stack import VRTStack
-from dolphin.utils import upsample_nearest
-
-logger = get_log()
-
-__all__ = ["run_evd_sequential", "compress"]
-
-
-def run_evd_sequential(
-    *,
-    slc_vrt_file: Filename,
-    # weight_file: Filename,
-    output_folder: Filename,
-    half_window: dict,
-    strides: dict = {"x": 1, "y": 1},
-    ministack_size: int = 10,
-    mask_file: Optional[Filename] = None,
-    ps_mask_file: Optional[Filename] = None,
-    beta: float = 0.01,
-    max_bytes: float = 32e6,
-    n_workers: int = 1,
-    gpu_enabled: bool = True,
-):
-    """Estimate wrapped phase using batches of ministacks."""
-    output_folder = Path(output_folder)
-    v_all = VRTStack.from_vrt_file(slc_vrt_file)
-    file_list_all = v_all.file_list
-    date_list_all = v_all.dates
-
-    logger.info(f"{v_all}: from {v_all.file_list[0]} to {v_all.file_list[-1]}")
-
-    # Map of {ministack_index: [output_slc_files]}
-    output_slc_files: Dict[int, List] = defaultdict(list)
-    comp_slc_files: List[Path] = []
-    tcorr_files: List[Path] = []
-
-    nrows, ncols = v_all.shape[-2:]
-    if mask_file is not None:
-        nodata_mask = io.load_gdal(mask_file).astype(bool)
-    else:
-        nodata_mask = np.zeros((nrows, ncols), dtype=bool)
-
-    if ps_mask_file is not None:
-        ps_mask = io.load_gdal(ps_mask_file, masked=True)
-        # Fill the nodata values with false
-        ps_mask = ps_mask.astype(bool).filled(False)
-    else:
-        ps_mask = np.zeros_like(nodata_mask)
-
-    xhalf, yhalf = half_window["x"], half_window["y"]
-    xs, ys = strides["x"], strides["y"]
-
-    # If we were passed any compressed SLCs in `file_list_all`,
-    # then we want the first `mini_idx` start at the first non-compressed index
-    first_non_comp_idx = 0
-    for filename in file_list_all:
-        if not Path(filename).name.startswith("compressed"):
-            break
-        first_non_comp_idx += 1
-
-    # Solve each ministack using the current chunk (and the previous compressed SLCs)
-    ministack_starts = range(0, len(file_list_all), ministack_size)
-    for mini_idx, full_stack_idx in enumerate(
-        ministack_starts, start=first_non_comp_idx
-    ):
-        cur_slice = slice(full_stack_idx, full_stack_idx + ministack_size)
-        cur_files = file_list_all[cur_slice].copy()
-        cur_dates = date_list_all[cur_slice].copy()
-
-        # Make the current ministack output folder using the start/end dates
-        d0 = cur_dates[first_non_comp_idx][0]
-        d1 = cur_dates[-1][0]
-        start_end = io._format_date_pair(d0, d1)
-        cur_output_folder = output_folder / start_end
-        cur_output_folder.mkdir(parents=True, exist_ok=True)
-
-        msg = (
-            f"Processing {len(cur_files) - first_non_comp_idx} SLCs +"
-            f" {len(comp_slc_files) + first_non_comp_idx} compressed SLCs. "
-        )
-        msg += f"Output folder: {cur_output_folder}"
-        logger.info(msg)
-        # Add the existing compressed SLC files to the start
-        cur_files = comp_slc_files + cur_files
-        cur_vrt = VRTStack(
-            cur_files,
-            outfile=cur_output_folder / f"{start_end}.vrt",
-            sort_files=False,
-            subdataset=v_all.subdataset,
-        )
-        # Create the background writer for this ministack
-        writer = io.Writer()
-
-        # mini_idx is first non-compressed SLC
-        logger.info(
-            f"{cur_vrt}: from {Path(cur_vrt.file_list[mini_idx]).name} to"
-            f" {Path(cur_vrt.file_list[-1]).name}"
-        )
-        # Set up the output folder with empty files to write into
-        cur_output_files = _setup_output_folder(
-            cur_vrt, driver="GTiff", start_idx=mini_idx, strides=strides
-        )
-        # Save these for the final adjustment later
-        # Keep the list of compressed SLCs to prepend to next VRTStack.file_list
-        output_slc_files[mini_idx] = cur_output_files
-
-        # Create the empty compressed SLC file
-        cur_comp_slc_file = cur_output_folder / f"compressed_{start_end}.tif"
-        io.write_arr(
-            arr=None,
-            like_filename=cur_vrt.outfile,
-            output_name=cur_comp_slc_file,
-            nbands=1,
-            # Note that the compressed SLC is the same size as the original SLC
-        )
-        comp_slc_files.append(cur_comp_slc_file)
-
-        # Create the empty compressed temporal coherence file
-        tcorr_file = cur_output_folder / f"tcorr_{start_end}.tif"
-        io.write_arr(
-            arr=None,
-            like_filename=cur_vrt.outfile,
-            output_name=tcorr_file,
-            nbands=1,
-            dtype=np.float32,
-            strides=strides,
-        )
-        tcorr_files.append(tcorr_file)
-
-        # Iterate over the ministack in blocks
-        # Note the overlap to redo the edge effects
-        # TODO: adjust the writing to avoid the overlap
-
-        # Note: dividing by len(stack) since cov is shape (rows, cols, nslc, nslc)
-        # so we need to load less to not overflow memory
-        stack_max_bytes = max_bytes / len(cur_vrt)
-        overlaps = (yhalf, xhalf)
-        block_gen = cur_vrt.iter_blocks(
-            overlaps=overlaps,
-            max_bytes=stack_max_bytes,
-            skip_empty=True,
-            # TODO: get the nodata value from the vrt stack
-            # this involves verifying that COMPASS correctly sets the nodata value
-        )
-        for cur_data, (rows, cols) in block_gen:
-            if np.all(cur_data == 0):
-                continue
-            cur_data = cur_data.astype(np.complex64)
-
-            # Run the phase linking process on the current ministack
-            try:
-                cur_mle_stack, tcorr = run_mle(
-                    cur_data,
-                    half_window=half_window,
-                    strides=strides,
-                    beta=beta,
-                    reference_idx=mini_idx,
-                    nodata_mask=nodata_mask[rows, cols],
-                    ps_mask=ps_mask[rows, cols],
-                    n_workers=n_workers,
-                    gpu_enabled=gpu_enabled,
-                )
-            except PhaseLinkRuntimeError as e:
-                # note: this is a warning instead of info, since it should
-                # get caught at the "skip_empty" step
-                logger.warning(f"Exception at ({rows}, {cols}): {e}")
-                continue
-
-            # Save each of the MLE estimates (ignoring the compressed SLCs)
-            assert len(cur_mle_stack[mini_idx:]) == len(cur_output_files)
-            # Get the location within the output file, shrinking down the slices
-            out_row_start = rows.start // ys
-            out_col_start = cols.start // xs
-            for img, f in zip(cur_mle_stack[mini_idx:], cur_output_files):
-                writer.queue_write(img, f, out_row_start, out_col_start)
-
-            # Save the temporal coherence blocks
-            writer.queue_write(tcorr, tcorr_file, out_row_start, out_col_start)
-
-            # Compress the ministack using only the non-compressed SLCs
-            cur_comp_slc = compress(
-                cur_data[mini_idx:],
-                cur_mle_stack[mini_idx:],
-            )
-            # Save the compressed SLC block
-            writer.queue_write(
-                cur_comp_slc, cur_comp_slc_file, out_row_start, out_col_start
-            )
-            # logger.debug(f"Saved compressed block SLC to {cur_comp_slc_file}")
-            # tqdm.write(" Finished block, loading next block.")
-
-        # Block until all the writers for this ministack have finished
-        logger.info(f"Waiting to write {writer.num_queued} blocks of data.")
-        writer.notify_finished()
-        logger.info(f"Finished ministack {mini_idx} of size {cur_vrt.shape}.")
-
-    ##############################################
-    # Set up the output folder with empty files to write into
-    # final_output_folder = output_folder / "final"
-    # final_output_folder.mkdir(parents=True, exist_ok=True)
-
-    # Average the temporal coherence files in each ministack
-    # TODO: do we want to include the date span in this filename?
-    output_tcorr_file = output_folder / "tcorr_average.tif"
-    # Find the offsets between stacks by doing a phase linking only compressed SLCs
-    # (But only if we have >1 ministacks. If only one, just rename the outputs)
-    if len(comp_slc_files) == 1:
-        # There was only one ministack, so we can skip this step
-        logger.info("Only one ministack, skipping offset calculation.")
-        assert len(output_slc_files) == 1
-        assert len(tcorr_files) == 1
-        for slc_fname in output_slc_files[first_non_comp_idx]:
-            slc_fname.rename(output_folder / slc_fname.name)
-
-        tcorr_files[0].rename(output_tcorr_file)
-
-        output_comp_slc_file = output_folder / comp_slc_files[0].name
-        comp_slc_files[0].rename(output_comp_slc_file)
-        return
-
-    # Compute the adjustments by running EVD on the compressed SLCs
-    comp_output_folder = output_folder / "adjustments"
-    comp_output_folder.mkdir(parents=True, exist_ok=True)
-    adjustment_vrt_stack = VRTStack(
-        comp_slc_files, outfile=comp_output_folder / "compressed_stack.vrt"
-    )
-
-    logger.info(f"Running EVD on compressed files: {adjustment_vrt_stack}")
-    adjusted_comp_slc_files = _setup_output_folder(
-        adjustment_vrt_stack, driver="GTiff", strides=strides
-    )
-
-    writer = io.Writer()
-    # Iterate over the ministack in blocks
-    # Note the overlap to redo the edge effects
-    block_gen = adjustment_vrt_stack.iter_blocks(
-        overlaps=(yhalf, xhalf),
-        # Note: dividing by len of stack because cov is shape (rows, cols, nslc, nslc)
-        max_bytes=max_bytes / len(adjustment_vrt_stack),
-        skip_empty=True,
-    )
-    for cur_data, (rows, cols) in block_gen:
-        msg = f"Processing block {rows.start}:{rows.stop}, {cols.start}:{cols.stop}"
-        logger.debug(msg)
-
-        # Run the phase linking process on the current adjustment stack
-        cur_mle_stack, tcorr = run_mle(
-            cur_data,
-            half_window=half_window,
-            strides=strides,
-            beta=beta,
-            reference_idx=0,
-            nodata_mask=nodata_mask[rows, cols],
-            ps_mask=None,  # PS mask doesn't matter for the adjustments
-            use_slc_amp=False,  # Make adjustments unit-amplitude
-            n_workers=n_workers,
-            gpu_enabled=gpu_enabled,
-        )
-
-        # Get the location within the output file, shrinking down the slices
-        out_row_start = rows.start // ys
-        out_col_start = cols.start // xs
-        # Save each of the MLE estimates (ignoring the compressed SLCs)
-        for img, f in zip(cur_mle_stack, adjusted_comp_slc_files):
-            writer.queue_write(img, f, out_row_start, out_col_start)
-        # Don't think I care about the temporal coherence here
-
-    writer.notify_finished()
-    # Compensate for the offsets between ministacks (aka "datum adjustments")
-    for mini_idx, slc_files in output_slc_files.items():
-        adjustment_fname = adjusted_comp_slc_files[mini_idx]
-
-        for slc_fname in slc_files:
-            logger.info(f"Compensating {slc_fname} with {adjustment_fname}")
-            outfile = output_folder / f"{slc_fname.name}"
-            VRTInterferogram(
-                ref_slc=slc_fname,
-                sec_slc=adjustment_fname,
-                path=outfile,
-                pixel_function="mul",
-            )
-
-    # Can pass the list of files to gdal_calc, which interprets it
-    # as a multi-band file
-    logger.info(f"Averaging temporal coherence files into: {output_tcorr_file}")
-    gdal_calc.Calc(
-        NoDataValue=nan,
-        format="GTiff",
-        outfile=fspath(output_tcorr_file),
-        type="Float32",
-        quiet=True,
-        overwrite=True,
-        creation_options=io.DEFAULT_TIFF_OPTIONS,
-        A=tcorr_files,
-        calc="numpy.nanmean(A, axis=0)",
-    )
-
-
-def compress(
-    slc_stack: np.ndarray,
-    mle_estimate: np.ndarray,
-):
-    """Compress the stack of SLC data using the estimated phase.
-
-    Parameters
-    ----------
-    slc_stack : np.array
-        The stack of complex SLC data, shape (nslc, rows, cols)
-    mle_estimate : np.array
-        The estimated phase from [`run_mle`][dolphin.phase_link.mle.run_mle],
-        shape (nslc, rows // strides['y'], cols // strides['x'])
-
-    Returns
-    -------
-    np.array
-        The compressed SLC data, shape (rows, cols)
-    """
-    # If the output is downsampled, we need to make `mle_estimate` the same shape
-    # as the output
-    mle_estimate_upsampled = upsample_nearest(mle_estimate, slc_stack.shape[1:])
-    # For each pixel, project the SLCs onto the (normalized) estimated phase
-    # by performing a pixel-wise complex dot product
-    mle_norm = np.linalg.norm(mle_estimate_upsampled, axis=0)
-    # Avoid divide by zero (there may be 0s at the upsampled boundary)
-    mle_norm[mle_norm == 0] = np.nan
-    with warnings.catch_warnings():
-        warnings.filterwarnings("ignore", message="invalid value encountered")
-        return (
-            np.nansum(slc_stack * np.conjugate(mle_estimate_upsampled), axis=0)
-            / mle_norm
-        )
-
-
-def _setup_output_folder(
-    vrt_stack,
-    driver: str = "GTiff",
-    dtype="complex64",
-    start_idx: int = 0,
-    strides: Dict[str, int] = {"y": 1, "x": 1},
-    creation_options: Optional[List] = None,
-) -> List[Path]:
-    """Create empty output files for each band after `start_idx` in `vrt_stack`.
-
-    Also creates an empty file for the compressed SLC.
-    Used to prepare output for block processing.
-
-    Parameters
-    ----------
-    vrt_stack : VRTStack
-        object containing the current stack of SLCs
-    driver : str, optional
-        Name of GDAL driver, by default "GTiff"
-    dtype : str, optional
-        Numpy datatype of output files, by default "complex64"
-    start_idx : int, optional
-        Index of vrt_stack to begin making output files.
-        This should match the ministack index to avoid re-creating the
-        past compressed SLCs.
-    strides : Dict[str, int], optional
-        Strides to use when creating the empty files, by default {"y": 1, "x": 1}
-        Larger strides will create smaller output files, computed using
-        [dolphin.io.compute_out_shape][]
-    creation_options : list, optional
-        List of options to pass to the GDAL driver, by default None
-
-    Returns
-    -------
-    List[Path]
-        List of saved empty files.
-    """
-    output_folder = vrt_stack.outfile.parent
-
-    date_strs: List[str] = []
-    for d in vrt_stack.dates[start_idx:]:
-        if len(d) == 1:
-            # normal SLC files will have a single date
-            s = d[0].strftime(io.DEFAULT_DATETIME_FORMAT)
-        else:
-            # Compressed SLCs will have 2 dates in the name marking the start and end
-            s = io._format_date_pair(d[0], d[1])
-        date_strs.append(s)
-
-    output_files = []
-    for filename in date_strs:
-        slc_name = Path(filename).stem
-        output_path = output_folder / f"{slc_name}.slc.tif"
-
-        io.write_arr(
-            arr=None,
-            like_filename=vrt_stack.outfile,
-            output_name=output_path,
-            driver=driver,
-            nbands=1,
-            dtype=dtype,
-            strides=strides,
-            options=creation_options,
-        )
-
-        output_files.append(output_path)
-    return output_files
diff --git a/src/dolphin/stack.py b/src/dolphin/stack.py
index c317855b..f86f3c64 100755
--- a/src/dolphin/stack.py
+++ b/src/dolphin/stack.py
@@ -3,7 +3,7 @@
 from os import fspath
 from pathlib import Path
 from pprint import pformat
-from typing import Generator, Optional, Sequence, Tuple
+from typing import Generator, List, Optional, Sequence, Tuple
 
 import numpy as np
 from osgeo import gdal
@@ -13,7 +13,7 @@
 from dolphin._types import Filename
 
 gdal.UseExceptions()
-logger = get_log()
+logger = get_log(__name__)
 
 
 DEFAULT_BLOCK_BYTES = 32e6
@@ -80,10 +80,10 @@ def __init__(
                 "To create from an existing VRT, use the `from_vrt_file` method."
             )
 
-        files = [Path(f) for f in file_list]
+        files: List[Filename] = [Path(f) for f in file_list]
         self._use_abs_path = use_abs_path
         if use_abs_path:
-            files = [p.resolve() for p in files]
+            files = [utils._resolve_gdal_path(p) for p in files]
         # Extract the date/datetimes from the filenames
         dates = [utils.get_dates(f, fmt=file_date_fmt) for f in files]
         if sort_files:
diff --git a/src/dolphin/stitching.py b/src/dolphin/stitching.py
index 1372a0a2..0998619e 100644
--- a/src/dolphin/stitching.py
+++ b/src/dolphin/stitching.py
@@ -14,7 +14,7 @@
 from dolphin._log import get_log
 from dolphin._types import Filename
 
-logger = get_log()
+logger = get_log(__name__)
 
 
 def merge_by_date(
diff --git a/src/dolphin/unwrap.py b/src/dolphin/unwrap.py
index f480f8bf..40f7ef19 100644
--- a/src/dolphin/unwrap.py
+++ b/src/dolphin/unwrap.py
@@ -14,7 +14,7 @@
 from dolphin.io import get_raster_xysize
 from dolphin.utils import numpy_to_gdal_type
 
-logger = get_log()
+logger = get_log(__name__)
 
 gdal.UseExceptions()
 
diff --git a/src/dolphin/utils.py b/src/dolphin/utils.py
index 68bc746d..a14ddd54 100644
--- a/src/dolphin/utils.py
+++ b/src/dolphin/utils.py
@@ -13,7 +13,7 @@
 from dolphin._types import Filename
 
 gdal.UseExceptions()
-logger = get_log()
+logger = get_log(__name__)
 
 
 def numpy_to_gdal_type(np_dtype: DTypeLike) -> int:
@@ -103,20 +103,26 @@ def _get_path_from_gdal_str(name: Filename) -> Path:
     return Path(p)
 
 
-def _resolve_gdal_path(gdal_str):
+def _resolve_gdal_path(gdal_str: Filename) -> Filename:
     """Resolve the file portion of a gdal-openable string to an absolute path."""
     s = str(gdal_str)
     if s.startswith("DERIVED_SUBDATASET"):
         # like DERIVED_SUBDATASET:AMPLITUDE:slc_filepath.tif
         file_part = s.split(":")[-1]
+        is_gdal_str = True
     elif ":" in s and (s.startswith("NETCDF") or s.startswith("HDF")):
         # like NETCDF:"slc_filepath.nc":slc_var
         file_part = s.split(":")[1]
+        is_gdal_str = True
+    else:
+        file_part = s
+        is_gdal_str = False
 
     # strip quotes to add back in after
     file_part = file_part.strip('"').strip("'")
     file_part_resolved = Path(file_part).resolve()
-    return gdal_str.replace(file_part, str(file_part_resolved))
+    resolved = s.replace(file_part, str(file_part_resolved))
+    return Path(resolved) if not is_gdal_str else resolved
 
 
 def _date_format_to_regex(date_format):
diff --git a/src/dolphin/workflows/__init__.py b/src/dolphin/workflows/__init__.py
index e7b4de65..a63bcdbc 100644
--- a/src/dolphin/workflows/__init__.py
+++ b/src/dolphin/workflows/__init__.py
@@ -1,4 +1,4 @@
 """Tools for creating workflows and running displacement workflows."""
 from ._enums import *  # noqa
-from ._grouping import *  # noqa
+from ._utils import *  # noqa
 from .config import *  # noqa
diff --git a/src/dolphin/workflows/_config_cli.py b/src/dolphin/workflows/_config_cli.py
index 91c27ac2..fefaabe4 100644
--- a/src/dolphin/workflows/_config_cli.py
+++ b/src/dolphin/workflows/_config_cli.py
@@ -4,7 +4,7 @@
 from pathlib import Path
 from typing import List, Optional, Tuple, Union
 
-from .config import OPERA_DATASET_NAME, InterferogramNetworkType, Workflow
+from .config import OPERA_DATASET_NAME, InterferogramNetworkType, Workflow, WorkflowName
 
 
 def create_config(
@@ -28,12 +28,15 @@ def create_config(
             network_type=InterferogramNetworkType.MANUAL_INDEX,
             indexes=[(0, -1)],
         )
+        workflow_name = WorkflowName.SINGLE
         # Override the ministack size so that only one phase linking is run
         ministack_size = 1000
     else:
         interferogram_network = {}  # Use default
+        workflow_name = WorkflowName.STACK
 
     cfg = Workflow(
+        workflow_name=workflow_name,
         inputs=dict(
             cslc_file_list=slc_files,
             mask_files=mask_files,
diff --git a/src/dolphin/workflows/_grouping.py b/src/dolphin/workflows/_grouping.py
deleted file mode 100644
index ea1a2800..00000000
--- a/src/dolphin/workflows/_grouping.py
+++ /dev/null
@@ -1,82 +0,0 @@
-import itertools
-import re
-from pathlib import Path
-from typing import Dict, List, Pattern, Sequence, Union
-
-from dolphin._log import get_log
-from dolphin._types import Filename
-
-from .config import OPERA_BURST_RE
-
-logger = get_log(__name__)
-
-__all__ = ["group_by_burst"]
-
-
-def group_by_burst(
-    file_list: Sequence[Filename],
-    burst_id_fmt: Union[str, Pattern[str]] = OPERA_BURST_RE,
-    minimum_slcs: int = 2,
-) -> Dict[str, List[Path]]:
-    """Group Sentinel CSLC files by burst.
-
-    Parameters
-    ----------
-    file_list: List[Filename]
-        path to folder containing CSLC files
-    burst_id_fmt: str
-        format of the burst id in the filename.
-        Default is [OPERA_BURST_RE][]
-    minimum_slcs: int
-        Minimum number of SLCs needed to run the workflow for each burst.
-        If there are fewer SLCs in a burst, it will be skipped and
-        a warning will be logged.
-
-    Returns
-    -------
-    dict
-        key is the burst id of the SLC acquisition
-        Value is a list of Paths on that burst:
-        {
-            't087_185678_iw2': [Path(...), Path(...),],
-            't087_185678_iw3': [Path(...),... ],
-        }
-    """
-
-    def get_burst_id(filename):
-        m = re.search(burst_id_fmt, str(filename))
-        if not m:
-            raise ValueError(f"Could not parse burst id from {filename}")
-        return m.group()
-
-    def sort_by_burst_id(file_list):
-        """Sort files by burst id."""
-        burst_ids = [get_burst_id(f) for f in file_list]
-        file_burst_tups = sorted(
-            [(Path(f), d) for f, d in zip(file_list, burst_ids)],
-            # use the date or dates as the key
-            key=lambda f_d_tuple: f_d_tuple[1],  # type: ignore
-        )
-        # Unpack the sorted pairs with new sorted values
-        file_list, burst_ids = zip(*file_burst_tups)  # type: ignore
-        return file_list
-
-    sorted_file_list = sort_by_burst_id(file_list)
-    # Now collapse into groups, sorted by the burst_id
-    grouped_images = {
-        burst_id: list(g)
-        for burst_id, g in itertools.groupby(
-            sorted_file_list, key=lambda x: get_burst_id(x)
-        )
-    }
-    # Make sure that each burst has at least the minimum number of SLCs
-    out = {}
-    for burst_id, slc_list in grouped_images.items():
-        if len(slc_list) < minimum_slcs:
-            logger.warning(
-                f"Skipping burst {burst_id} because it has only {len(slc_list)} SLCs."
-                f"Minimum number of SLCs is {minimum_slcs}"
-            )
-        else:
-            out[burst_id] = slc_list
-    return out
diff --git a/src/dolphin/workflows/_run_cli.py b/src/dolphin/workflows/_run_cli.py
index 0663706d..f790b597 100644
--- a/src/dolphin/workflows/_run_cli.py
+++ b/src/dolphin/workflows/_run_cli.py
@@ -4,7 +4,7 @@
 
 from dolphin._log import log_runtime
 
-from ._enums import WorkflowName
+from . import s1_disp
 from .config import Workflow
 
 
@@ -23,17 +23,8 @@ def run(config_file: str, debug: bool = False, log_file: Optional[str] = None):
     """
     cfg = Workflow.from_yaml(config_file)
     cfg.create_dir_tree(debug=debug)
-    if cfg.workflow_name == "stack":
-        from dolphin.workflows import s1_disp
 
-        s1_disp.run(cfg, debug=debug, log_file=log_file)
-    elif cfg.workflow_name == "single":
-        raise NotImplementedError("Single interferogram workflow not yet implemented")
-    else:
-        choices = WorkflowName.__members__.values()
-        raise ValueError(
-            f"Unknown workflow name: {cfg.workflow_name}. Must be one of {choices}"
-        )
+    s1_disp.run(cfg, debug=debug, log_file=log_file)
 
 
 def get_parser(subparser=None, subcommand_name="run"):
diff --git a/src/dolphin/workflows/_utils.py b/src/dolphin/workflows/_utils.py
new file mode 100644
index 00000000..45036c49
--- /dev/null
+++ b/src/dolphin/workflows/_utils.py
@@ -0,0 +1,157 @@
+import itertools
+import re
+from pathlib import Path
+from typing import Dict, List, Optional, Pattern, Sequence, Union
+
+from dolphin import io
+from dolphin._log import get_log
+from dolphin._types import Filename
+
+from .config import OPERA_BURST_RE
+
+logger = get_log(__name__)
+
+__all__ = ["group_by_burst", "setup_output_folder"]
+
+
+def group_by_burst(
+    file_list: Sequence[Filename],
+    burst_id_fmt: Union[str, Pattern[str]] = OPERA_BURST_RE,
+    minimum_slcs: int = 2,
+) -> Dict[str, List[Path]]:
+    """Group Sentinel CSLC files by burst.
+
+    Parameters
+    ----------
+    file_list: List[Filename]
+        path to folder containing CSLC files
+    burst_id_fmt: str
+        format of the burst id in the filename.
+        Default is [OPERA_BURST_RE][]
+    minimum_slcs: int
+        Minimum number of SLCs needed to run the workflow for each burst.
+        If there are fewer SLCs in a burst, it will be skipped and
+        a warning will be logged.
+
+    Returns
+    -------
+    dict
+        key is the burst id of the SLC acquisition
+        Value is a list of Paths on that burst:
+        {
+            't087_185678_iw2': [Path(...), Path(...),],
+            't087_185678_iw3': [Path(...),... ],
+        }
+    """
+
+    def get_burst_id(filename):
+        m = re.search(burst_id_fmt, str(filename))
+        if not m:
+            raise ValueError(f"Could not parse burst id from {filename}")
+        return m.group()
+
+    def sort_by_burst_id(file_list):
+        """Sort files by burst id."""
+        burst_ids = [get_burst_id(f) for f in file_list]
+        file_burst_tups = sorted(
+            [(Path(f), d) for f, d in zip(file_list, burst_ids)],
+            # use the date or dates as the key
+            key=lambda f_d_tuple: f_d_tuple[1],  # type: ignore
+        )
+        # Unpack the sorted pairs with new sorted values
+        file_list, burst_ids = zip(*file_burst_tups)  # type: ignore
+        return file_list
+
+    sorted_file_list = sort_by_burst_id(file_list)
+    # Now collapse into groups, sorted by the burst_id
+    grouped_images = {
+        burst_id: list(g)
+        for burst_id, g in itertools.groupby(
+            sorted_file_list, key=lambda x: get_burst_id(x)
+        )
+    }
+    # Make sure that each burst has at least the minimum number of SLCs
+    out = {}
+    for burst_id, slc_list in grouped_images.items():
+        if len(slc_list) < minimum_slcs:
+            logger.warning(
+                f"Skipping burst {burst_id} because it has only {len(slc_list)} SLCs."
+                f"Minimum number of SLCs is {minimum_slcs}"
+            )
+        else:
+            out[burst_id] = slc_list
+    return out
+
+
+def setup_output_folder(
+    vrt_stack,
+    driver: str = "GTiff",
+    dtype="complex64",
+    start_idx: int = 0,
+    strides: Dict[str, int] = {"y": 1, "x": 1},
+    creation_options: Optional[List] = None,
+    output_folder: Optional[Path] = None,
+) -> List[Path]:
+    """Create empty output files for each band after `start_idx` in `vrt_stack`.
+
+    Also creates an empty file for the compressed SLC.
+    Used to prepare output for block processing.
+
+    Parameters
+    ----------
+    vrt_stack : VRTStack
+        object containing the current stack of SLCs
+    driver : str, optional
+        Name of GDAL driver, by default "GTiff"
+    dtype : str, optional
+        Numpy datatype of output files, by default "complex64"
+    start_idx : int, optional
+        Index of vrt_stack to begin making output files.
+        This should match the ministack index to avoid re-creating the
+        past compressed SLCs.
+    strides : Dict[str, int], optional
+        Strides to use when creating the empty files, by default {"y": 1, "x": 1}
+        Larger strides will create smaller output files, computed using
+        [dolphin.io.compute_out_shape][]
+    creation_options : list, optional
+        List of options to pass to the GDAL driver, by default None
+    output_folder : Path, optional
+        Path to output folder, by default None
+        If None, will use the same folder as the first SLC in `vrt_stack`
+
+    Returns
+    -------
+    List[Path]
+        List of saved empty files.
+    """
+    if output_folder is None:
+        output_folder = vrt_stack.outfile.parent
+
+    date_strs: List[str] = []
+    for d in vrt_stack.dates[start_idx:]:
+        if len(d) == 1:
+            # normal SLC files will have a single date
+            s = d[0].strftime(io.DEFAULT_DATETIME_FORMAT)
+        else:
+            # Compressed SLCs will have 2 dates in the name marking the start and end
+            s = io._format_date_pair(d[0], d[1])
+        date_strs.append(s)
+
+    output_files = []
+    for filename in date_strs:
+        slc_name = Path(filename).stem
+        output_path = output_folder / f"{slc_name}.slc.tif"
+
+        io.write_arr(
+            arr=None,
+            like_filename=vrt_stack.outfile,
+            output_name=output_path,
+            driver=driver,
+            nbands=1,
+            dtype=dtype,
+            strides=strides,
+            options=creation_options,
+        )
+
+        output_files.append(output_path)
+    return output_files
diff --git a/src/dolphin/workflows/config.py b/src/dolphin/workflows/config.py
index 4b7706ff..ca568e10 100644
--- a/src/dolphin/workflows/config.py
+++ b/src/dolphin/workflows/config.py
@@ -34,7 +34,7 @@
     "Workflow",
 ]
 
-logger = get_log()
+logger = get_log(__name__)
 
 # Specific to OPERA CSLC products:
 OPERA_DATASET_NAME = "science/SENTINEL1/CSLC/grids/VV"
diff --git a/src/dolphin/workflows/s1_disp.py b/src/dolphin/workflows/s1_disp.py
index 155659ec..0445ed1d 100755
--- a/src/dolphin/workflows/s1_disp.py
+++ b/src/dolphin/workflows/s1_disp.py
@@ -7,7 +7,7 @@
 from dolphin.interferogram import VRTInterferogram
 
 from . import stitch_and_unwrap, wrapped_phase
-from ._grouping import group_by_burst
+from ._utils import group_by_burst
 from .config import Workflow
 
 
diff --git a/src/dolphin/workflows/sequential.py b/src/dolphin/workflows/sequential.py
new file mode 100644
index 00000000..88fa2eb0
--- /dev/null
+++ b/src/dolphin/workflows/sequential.py
@@ -0,0 +1,217 @@
+"""Estimate wrapped phase using batches of ministacks.
+
+References
+----------
+    [1] Ansari, H., De Zan, F., & Bamler, R. (2017). Sequential estimator: Toward
+    efficient InSAR time series analysis. IEEE Transactions on Geoscience and
+    Remote Sensing, 55(10), 5637-5652.
+"""
+from collections import defaultdict
+from math import nan
+from os import fspath
+from pathlib import Path
+from typing import Dict, List, Optional
+
+import numpy as np
+from osgeo_utils import gdal_calc
+
+from dolphin import io
+from dolphin._log import get_log
+from dolphin._types import Filename
+from dolphin.interferogram import VRTInterferogram
+from dolphin.phase_link import run_mle
+from dolphin.stack import VRTStack
+
+from ._utils import setup_output_folder
+from .single import run_evd_single
+
+logger = get_log(__name__)
+
+__all__ = ["run_evd_sequential"]
+
+
+def run_evd_sequential(
+    *,
+    slc_vrt_file: Filename,
+    # weight_file: Filename,
+    output_folder: Filename,
+    half_window: dict,
+    strides: dict = {"x": 1, "y": 1},
+    ministack_size: int = 10,
+    mask_file: Optional[Filename] = None,
+    ps_mask_file: Optional[Filename] = None,
+    beta: float = 0.01,
+    max_bytes: float = 32e6,
+    n_workers: int = 1,
+    gpu_enabled: bool = True,
+):
+    """Estimate wrapped phase using batches of ministacks."""
+    output_folder = Path(output_folder)
+    v_all = VRTStack.from_vrt_file(slc_vrt_file)
+    file_list_all = v_all.file_list
+    date_list_all = v_all.dates
+
+    logger.info(f"{v_all}: from {v_all.file_list[0]} to {v_all.file_list[-1]}")
+
+    # Map of {ministack_index: [output_slc_files]}
+    output_slc_files: Dict[int, List] = defaultdict(list)
+    comp_slc_files: List[Path] = []
+    tcorr_files: List[Path] = []
+
+    nrows, ncols = v_all.shape[-2:]
+    if mask_file is not None:
+        nodata_mask = io.load_gdal(mask_file).astype(bool)
+    else:
+        nodata_mask = np.zeros((nrows, ncols), dtype=bool)
+
+    if ps_mask_file is not None:
+        ps_mask = io.load_gdal(ps_mask_file, masked=True)
+        # Fill the nodata values with false
+        ps_mask = ps_mask.astype(bool).filled(False)
+    else:
+        ps_mask = np.zeros_like(nodata_mask)
+
+    xhalf, yhalf = half_window["x"], half_window["y"]
+    xs, ys = strides["x"], strides["y"]
+
+    # Solve each ministack using the current chunk (and the previous compressed SLCs)
+    ministack_starts = range(0, len(file_list_all), ministack_size)
+    for mini_idx, full_stack_idx in enumerate(ministack_starts):
+        cur_slice = slice(full_stack_idx, full_stack_idx + ministack_size)
+        cur_files = file_list_all[cur_slice].copy()
+        cur_dates = date_list_all[cur_slice].copy()
+
+        # 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)
+        cur_output_folder = output_folder / start_end
+        cur_output_folder.mkdir(parents=True, exist_ok=True)
+
+        msg = f"Processing {len(cur_files)} SLCs."
+        msg += f"Output folder: {cur_output_folder}"
+        logger.info(msg)
+        # Add the existing compressed SLC files to the start
+        cur_files = comp_slc_files + cur_files
+        cur_vrt = VRTStack(
+            cur_files,
+            outfile=cur_output_folder / f"{start_end}.vrt",
+            sort_files=False,
+            subdataset=v_all.subdataset,
+        )
+        cur_output_files, cur_comp_slc_file, tcorr_file = run_evd_single(
+            slc_vrt_file=cur_vrt,
+            output_folder=cur_output_folder,
+            half_window=half_window,
+            strides=strides,
+            reference_idx=mini_idx,
+            mask_file=mask_file,
+            ps_mask_file=ps_mask_file,
+            beta=beta,
+            max_bytes=max_bytes,
+            n_workers=n_workers,
+            gpu_enabled=gpu_enabled,
+        )
+
+        output_slc_files[mini_idx] = cur_output_files
+        comp_slc_files.append(cur_comp_slc_file)
+        tcorr_files.append(tcorr_file)
+
+    ##############################################
+    # Set up the output folder with empty files to write into
+
+    # Average the temporal coherence files in each ministack
+    # TODO: do we want to include the date span in this filename?
+    output_tcorr_file = output_folder / "tcorr_average.tif"
+    # Find the offsets between stacks by doing a phase linking only compressed SLCs
+    # (But only if we have >1 ministacks. If only one, just rename the outputs)
+    if len(comp_slc_files) == 1:
+        # There was only one ministack, so we can skip this step
+        logger.info("Only one ministack, skipping offset calculation.")
+        assert len(output_slc_files) == 1
+        assert len(tcorr_files) == 1
+        for slc_fname in output_slc_files[0]:
+            slc_fname.rename(output_folder / slc_fname.name)
+
+        tcorr_files[0].rename(output_tcorr_file)
+
+        output_comp_slc_file = output_folder / comp_slc_files[0].name
+        comp_slc_files[0].rename(output_comp_slc_file)
+        return
+
+    # Compute the adjustments by running EVD on the compressed SLCs
+    comp_output_folder = output_folder / "adjustments"
+    comp_output_folder.mkdir(parents=True, exist_ok=True)
+    adjustment_vrt_stack = VRTStack(
+        comp_slc_files, outfile=comp_output_folder / "compressed_stack.vrt"
+    )
+
+    logger.info(f"Running EVD on compressed files: {adjustment_vrt_stack}")
+    adjusted_comp_slc_files = setup_output_folder(
+        adjustment_vrt_stack, driver="GTiff", strides=strides
+    )
+
+    writer = io.Writer()
+    # Iterate over the ministack in blocks
+    # Note the overlap to redo the edge effects
+    block_gen = adjustment_vrt_stack.iter_blocks(
+        overlaps=(yhalf, xhalf),
+        # Note: dividing by len of stack because cov is shape (rows, cols, nslc, nslc)
+        max_bytes=max_bytes / len(adjustment_vrt_stack),
+        skip_empty=True,
+    )
+    for cur_data, (rows, cols) in block_gen:
+        msg = f"Processing block {rows.start}:{rows.stop}, {cols.start}:{cols.stop}"
+        logger.debug(msg)
+
+        # Run the phase linking process on the current adjustment stack
+        cur_mle_stack, tcorr = run_mle(
+            cur_data,
+            half_window=half_window,
+            strides=strides,
+            beta=beta,
+            reference_idx=0,
+            nodata_mask=nodata_mask[rows, cols],
+            ps_mask=None,  # PS mask doesn't matter for the adjustments
+            use_slc_amp=False,  # Make adjustments unit-amplitude
+            n_workers=n_workers,
+            gpu_enabled=gpu_enabled,
+        )
+
+        # Get the location within the output file, shrinking down the slices
+        out_row_start = rows.start // ys
+        out_col_start = cols.start // xs
+        # Save each of the MLE estimates (ignoring the compressed SLCs)
+        for img, f in zip(cur_mle_stack, adjusted_comp_slc_files):
+            writer.queue_write(img, f, out_row_start, out_col_start)
+        # Don't think I care about the temporal coherence here
+
+    writer.notify_finished()
+    # Compensate for the offsets between ministacks (aka "datum adjustments")
+    for mini_idx, slc_files in output_slc_files.items():
+        adjustment_fname = adjusted_comp_slc_files[mini_idx]
+
+        for slc_fname in slc_files:
+            logger.info(f"Compensating {slc_fname} with {adjustment_fname}")
+            outfile = output_folder / f"{slc_fname.name}"
+            VRTInterferogram(
+                ref_slc=slc_fname,
+                sec_slc=adjustment_fname,
+                path=outfile,
+                pixel_function="mul",
+            )
+
+    # Can pass the list of files to gdal_calc, which interprets it
+    # as a multi-band file
+    logger.info(f"Averaging temporal coherence files into: {output_tcorr_file}")
+    gdal_calc.Calc(
+        NoDataValue=nan,
+        format="GTiff",
+        outfile=fspath(output_tcorr_file),
+        type="Float32",
+        quiet=True,
+        overwrite=True,
+        creation_options=io.DEFAULT_TIFF_OPTIONS,
+        A=tcorr_files,
+        calc="numpy.nanmean(A, axis=0)",
+    )
diff --git a/src/dolphin/workflows/single.py b/src/dolphin/workflows/single.py
new file mode 100644
index 00000000..4856560d
--- /dev/null
+++ b/src/dolphin/workflows/single.py
@@ -0,0 +1,192 @@
+"""Estimate wrapped phase for one updated SLC using the online algorithm.
+
+References
+----------
+    [1] Mirzaee, Sara, Falk Amelung, and Heresh Fattahi. "Non-linear phase
+    linking using joined distributed and persistent scatterers." Computers &
+    Geosciences (2022): 105291.
+"""
+from pathlib import Path
+from typing import Optional
+
+import numpy as np
+
+from dolphin import io
+from dolphin._log import get_log
+from dolphin._types import Filename
+from dolphin.phase_link import PhaseLinkRuntimeError, compress, run_mle
+from dolphin.stack import VRTStack
+
+from ._utils import setup_output_folder
+
+logger = get_log(__name__)
+
+__all__ = ["run_evd_single"]
+
+
+def run_evd_single(
+    *,
+    slc_vrt_file: Filename,
+    # weight_file: Filename,
+    output_folder: Filename,
+    half_window: dict,
+    strides: dict = {"x": 1, "y": 1},
+    reference_idx: int = 0,
+    mask_file: Optional[Filename] = None,
+    ps_mask_file: Optional[Filename] = None,
+    beta: float = 0.01,
+    max_bytes: float = 32e6,
+    n_workers: int = 1,
+    gpu_enabled: bool = True,
+):
+    """Estimate wrapped phase for one ministack."""
+    # TODO: extract common stuff between here and sequential
+    output_folder = Path(output_folder)
+    vrt = VRTStack.from_vrt_file(slc_vrt_file)
+    file_list_all = vrt.file_list
+    date_list_all = vrt.dates
+
+    logger.info(f"{vrt}: from {vrt.file_list[0]} to {vrt.file_list[-1]}")
+
+    nrows, ncols = vrt.shape[-2:]
+    if mask_file is not None:
+        nodata_mask = io.load_gdal(mask_file).astype(bool)
+    else:
+        nodata_mask = np.zeros((nrows, ncols), dtype=bool)
+
+    if ps_mask_file is not None:
+        ps_mask = io.load_gdal(ps_mask_file, masked=True)
+        # Fill the nodata values with false
+        ps_mask = ps_mask.astype(bool).filled(False)
+    else:
+        ps_mask = np.zeros_like(nodata_mask)
+
+    xhalf, yhalf = half_window["x"], half_window["y"]
+    xs, ys = strides["x"], strides["y"]
+
+    # If we were passed any compressed SLCs in `file_list_all`,
+    # then we want that index for when we create new compressed SLCs.
+    # We skip the old compressed SLCs to create new ones
+    first_non_comp_idx = 0
+    for filename in file_list_all:
+        if not Path(filename).name.startswith("compressed"):
+            break
+        first_non_comp_idx += 1
+
+    # Make the output folder using the start/end dates
+    d0 = date_list_all[first_non_comp_idx][0]
+    d1 = date_list_all[-1][0]
+    start_end = io._format_date_pair(d0, d1)
+
+    msg = (
+        f"Processing {len(file_list_all) - first_non_comp_idx} SLCs +"
+        f" {first_non_comp_idx} compressed SLCs. "
+    )
+    logger.info(msg)
+
+    # Create the background writer for this ministack
+    writer = io.Writer()
+
+    logger.info(
+        f"{vrt}: from {Path(vrt.file_list[first_non_comp_idx]).name} to"
+        f" {Path(vrt.file_list[-1]).name}"
+    )
+    # Set up the output folder with empty files to write into
+    output_slc_files = setup_output_folder(
+        vrt,
+        driver="GTiff",
+        start_idx=first_non_comp_idx,
+        strides=strides,
+        output_folder=output_folder,
+    )
+
+    # Create the empty compressed SLC file
+    comp_slc_file = output_folder / f"compressed_{start_end}.tif"
+    io.write_arr(
+        arr=None,
+        like_filename=vrt.outfile,
+        output_name=comp_slc_file,
+        nbands=1,
+        # Note that the compressed SLC is the same size as the original SLC
+    )
+
+    # Create the empty compressed temporal coherence file
+    tcorr_file = output_folder / f"tcorr_{start_end}.tif"
+    io.write_arr(
+        arr=None,
+        like_filename=vrt.outfile,
+        output_name=tcorr_file,
+        nbands=1,
+        dtype=np.float32,
+        strides=strides,
+    )
+
+    # Iterate over the stack in blocks
+    # Note the overlap to redo the edge effects
+    # TODO: adjust the writing to avoid the overlap
+
+    # Note: dividing by len(stack) since cov is shape (rows, cols, nslc, nslc)
+    # so we need to load less to not overflow memory
+    # TODO: more robust way to do this... or change this to a "loading_size",
+    # then calculate the max_bytes somehow based on testing
+    stack_max_bytes = max_bytes / len(vrt)
+    overlaps = (yhalf, xhalf)
+    block_gen = vrt.iter_blocks(
+        overlaps=overlaps,
+        max_bytes=stack_max_bytes,
+        skip_empty=True,
+        # TODO: get the nodata value from the vrt stack
+        # this involves verifying that COMPASS correctly sets the nodata value
+    )
+    for cur_data, (rows, cols) in block_gen:
+        if np.all(cur_data == 0):
+            continue
+        cur_data = cur_data.astype(np.complex64)
+
+        # Run the phase linking process on the current ministack
+        try:
+            cur_mle_stack, tcorr = run_mle(
+                cur_data,
+                half_window=half_window,
+                strides=strides,
+                beta=beta,
+                reference_idx=reference_idx,
+                nodata_mask=nodata_mask[rows, cols],
+                ps_mask=ps_mask[rows, cols],
+                n_workers=n_workers,
+                gpu_enabled=gpu_enabled,
+            )
+        except PhaseLinkRuntimeError as e:
+            # note: this is a warning instead of info, since it should
+            # get caught at the "skip_empty" step
+            logger.warning(f"Exception at ({rows}, {cols}): {e}")
+            continue
+
+        # Save each of the MLE estimates (ignoring the compressed SLCs)
+        assert len(cur_mle_stack[first_non_comp_idx:]) == len(output_slc_files)
+        # Get the location within the output file, shrinking down the slices
+        out_row_start = rows.start // ys
+        out_col_start = cols.start // xs
+        for img, f in zip(cur_mle_stack[first_non_comp_idx:], output_slc_files):
+            writer.queue_write(img, f, out_row_start, out_col_start)
+
+        # Save the temporal coherence blocks
+        writer.queue_write(tcorr, tcorr_file, out_row_start, out_col_start)
+
+        # Compress the ministack using only the non-compressed SLCs
+        cur_comp_slc = compress(
+            cur_data[first_non_comp_idx:],
+            cur_mle_stack[first_non_comp_idx:],
+        )
+        # Save the compressed SLC block
+        # TODO: make a flag? We don't always need to save the compressed SLCs
+        writer.queue_write(cur_comp_slc, comp_slc_file, out_row_start, out_col_start)
+        # logger.debug(f"Saved compressed block SLC to {cur_comp_slc_file}")
+        # tqdm.write(" Finished block, loading next block.")
+
+    # Block until all the writers for this ministack have finished
+    logger.info(f"Waiting to write {writer.num_queued} blocks of data.")
+    writer.notify_finished()
+    logger.info(f"Finished ministack of size {vrt.shape}.")
+
+    return output_slc_files, comp_slc_file, tcorr_file
diff --git a/src/dolphin/workflows/wrapped_phase.py b/src/dolphin/workflows/wrapped_phase.py
index 85db4660..e5808641 100644
--- a/src/dolphin/workflows/wrapped_phase.py
+++ b/src/dolphin/workflows/wrapped_phase.py
@@ -1,10 +1,11 @@
 from pathlib import Path
 from typing import List
 
-from dolphin import ps, sequential, stack
+from dolphin import ps, stack
 from dolphin._log import get_log, log_runtime
 from dolphin.interferogram import Network, VRTInterferogram
 
+from . import sequential, single
 from .config import Workflow
 
 
@@ -67,19 +68,34 @@ def run(cfg: Workflow, debug: bool = False) -> List[VRTInterferogram]:
         logger.info(f"Skipping EVD step, {len(existing_files)} files already exist")
     else:
         logger.info(f"Running sequential EMI step in {pl_path}")
-        sequential.run_evd_sequential(
-            slc_vrt_file=vrt_stack.outfile,
-            output_folder=pl_path,
-            half_window=cfg.phase_linking.half_window.dict(),
-            strides=cfg.outputs.strides,
-            ministack_size=cfg.phase_linking.ministack_size,
-            # mask_file=cfg.inputs.mask_file,
-            ps_mask_file=cfg.ps_options.output_file,
-            max_bytes=cfg.worker_settings.max_ram_gb * 1e9,
-            n_workers=cfg.worker_settings.n_workers,
-            gpu_enabled=cfg.worker_settings.gpu_enabled,
-            beta=cfg.phase_linking.beta,
-        )
+        if cfg.workflow_name == "single":
+            single.run_evd_single(
+                slc_vrt_file=vrt_stack.outfile,
+                output_folder=pl_path,
+                half_window=cfg.phase_linking.half_window.dict(),
+                strides=cfg.outputs.strides,
+                reference_idx=0,
+                # mask_file=cfg.inputs.mask_file,
+                ps_mask_file=cfg.ps_options.output_file,
+                max_bytes=cfg.worker_settings.max_ram_gb * 1e9,
+                n_workers=cfg.worker_settings.n_workers,
+                gpu_enabled=cfg.worker_settings.gpu_enabled,
+                beta=cfg.phase_linking.beta,
+            )
+        else:
+            sequential.run_evd_sequential(
+                slc_vrt_file=vrt_stack.outfile,
+                output_folder=pl_path,
+                half_window=cfg.phase_linking.half_window.dict(),
+                strides=cfg.outputs.strides,
+                ministack_size=cfg.phase_linking.ministack_size,
+                # mask_file=cfg.inputs.mask_file,
+                ps_mask_file=cfg.ps_options.output_file,
+                max_bytes=cfg.worker_settings.max_ram_gb * 1e9,
+                n_workers=cfg.worker_settings.n_workers,
+                gpu_enabled=cfg.worker_settings.gpu_enabled,
+                beta=cfg.phase_linking.beta,
+            )
 
     # ###################################################
     # 4. Form interferograms from estimated wrapped phase


Test patch:
diff --git a/tests/test_workflows_s1_disp.py b/tests/test_workflows_s1_disp.py
index 1ec9fd51..e69de29b 100644
--- a/tests/test_workflows_s1_disp.py
+++ b/tests/test_workflows_s1_disp.py
@@ -1,47 +0,0 @@
-import random
-from itertools import chain
-from pathlib import Path
-
-import pytest
-
-from dolphin.workflows import group_by_burst
-
-
-def test_group_by_burst():
-    expected = {
-        "t087_185678_iw2": [
-            Path("t087_185678_iw2/20180210/t087_185678_iw2_20180210_VV.h5"),
-            Path("t087_185678_iw2/20180318/t087_185678_iw2_20180318_VV.h5"),
-            Path("t087_185678_iw2/20180423/t087_185678_iw2_20180423_VV.h5"),
-        ],
-        "t087_185678_iw3": [
-            Path("t087_185678_iw3/20180210/t087_185678_iw3_20180210_VV.h5"),
-            Path("t087_185678_iw3/20180318/t087_185678_iw3_20180318_VV.h5"),
-            Path("t087_185678_iw3/20180517/t087_185678_iw3_20180517_VV.h5"),
-        ],
-        "t087_185679_iw1": [
-            Path("t087_185679_iw1/20180210/t087_185679_iw1_20180210_VV.h5"),
-            Path("t087_185679_iw1/20180318/t087_185679_iw1_20180318_VV.h5"),
-        ],
-    }
-    in_files = list(chain.from_iterable(expected.values()))
-
-    assert group_by_burst(in_files) == expected
-
-    # Any order should work
-    random.shuffle(in_files)
-    # but the order of the lists of each key may be different
-    for burst, file_list in group_by_burst(in_files).items():
-        assert sorted(file_list) == sorted(expected[burst])
-
-
-def test_group_by_burst_non_opera():
-    with pytest.raises(ValueError, match="Could not parse burst id"):
-        group_by_burst(["20200101.slc", "20200202.slc"])
-        # A combination should still error
-        group_by_burst(
-            [
-                "20200101.slc",
-                Path("t087_185679_iw1/20180210/t087_185679_iw1_20180210_VV.h5"),
-            ]
-        )
diff --git a/tests/test_sequential.py b/tests/test_workflows_sequential.py
similarity index 65%
rename from tests/test_sequential.py
rename to tests/test_workflows_sequential.py
index 7f6b31b0..a5222986 100644
--- a/tests/test_sequential.py
+++ b/tests/test_workflows_sequential.py
@@ -1,11 +1,10 @@
-import numpy as np
 import numpy.testing as npt
 import pytest
-from osgeo import gdal
 
-from dolphin import io, sequential, stack
+from dolphin import io, stack
 from dolphin.phase_link import mle, simulate
 from dolphin.utils import gpu_is_available
+from dolphin.workflows import sequential
 
 GPU_AVAILABLE = gpu_is_available()
 simulate._seed(1234)
@@ -16,54 +15,6 @@
 )
 
 
-def test_setup_output_folder(tmpdir, tiled_file_list):
-    vrt_stack = stack.VRTStack(tiled_file_list, outfile=tmpdir / "stack.vrt")
-    out_file_list = sequential._setup_output_folder(
-        vrt_stack, driver="GTiff", dtype=np.complex64
-    )
-    for out_file in out_file_list:
-        assert out_file.exists()
-        assert out_file.suffix == ".tif"
-        assert out_file.parent == tmpdir
-        ds = gdal.Open(str(out_file))
-        assert ds.GetRasterBand(1).DataType == gdal.GDT_CFloat32
-        ds = None
-
-    out_file_list = sequential._setup_output_folder(
-        vrt_stack,
-        driver="GTiff",
-        dtype="float32",
-        start_idx=1,
-    )
-    assert len(out_file_list) == len(vrt_stack) - 1
-    for out_file in out_file_list:
-        assert out_file.exists()
-        ds = gdal.Open(str(out_file))
-        assert ds.GetRasterBand(1).DataType == gdal.GDT_Float32
-
-
-@pytest.mark.parametrize(
-    "strides", [{"x": 1, "y": 1}, {"x": 1, "y": 2}, {"x": 2, "y": 3}, {"x": 4, "y": 2}]
-)
-def test_setup_output_folder_strided(tmpdir, tiled_file_list, strides):
-    vrt_stack = stack.VRTStack(tiled_file_list, outfile=tmpdir / "stack.vrt")
-
-    out_file_list = sequential._setup_output_folder(
-        vrt_stack, driver="GTiff", dtype=np.complex64, strides=strides
-    )
-    rows, cols = vrt_stack.shape[-2:]
-    for out_file in out_file_list:
-        assert out_file.exists()
-        assert out_file.suffix == ".tif"
-        assert out_file.parent == tmpdir
-
-        ds = gdal.Open(str(out_file))
-        assert ds.GetRasterBand(1).DataType == gdal.GDT_CFloat32
-        assert ds.RasterXSize == cols // strides["x"]
-        assert ds.RasterYSize == rows // strides["y"]
-        ds = None
-
-
 @pytest.mark.parametrize("gpu_enabled", [True, False])
 def test_sequential_gtiff(tmp_path, slc_file_list, gpu_enabled):
     """Run through the sequential estimation with a GeoTIFF stack."""
diff --git a/tests/test_workflows_utils.py b/tests/test_workflows_utils.py
new file mode 100644
index 00000000..133982ee
--- /dev/null
+++ b/tests/test_workflows_utils.py
@@ -0,0 +1,97 @@
+import random
+from itertools import chain
+from pathlib import Path
+
+import numpy as np
+import pytest
+from osgeo import gdal
+
+from dolphin import stack
+from dolphin.workflows import group_by_burst
+from dolphin.workflows._utils import setup_output_folder
+
+
+def test_group_by_burst():
+    expected = {
+        "t087_185678_iw2": [
+            Path("t087_185678_iw2/20180210/t087_185678_iw2_20180210_VV.h5"),
+            Path("t087_185678_iw2/20180318/t087_185678_iw2_20180318_VV.h5"),
+            Path("t087_185678_iw2/20180423/t087_185678_iw2_20180423_VV.h5"),
+        ],
+        "t087_185678_iw3": [
+            Path("t087_185678_iw3/20180210/t087_185678_iw3_20180210_VV.h5"),
+            Path("t087_185678_iw3/20180318/t087_185678_iw3_20180318_VV.h5"),
+            Path("t087_185678_iw3/20180517/t087_185678_iw3_20180517_VV.h5"),
+        ],
+        "t087_185679_iw1": [
+            Path("t087_185679_iw1/20180210/t087_185679_iw1_20180210_VV.h5"),
+            Path("t087_185679_iw1/20180318/t087_185679_iw1_20180318_VV.h5"),
+        ],
+    }
+    in_files = list(chain.from_iterable(expected.values()))
+
+    assert group_by_burst(in_files) == expected
+
+    # Any order should work
+    random.shuffle(in_files)
+    # but the order of the lists of each key may be different
+    for burst, file_list in group_by_burst(in_files).items():
+        assert sorted(file_list) == sorted(expected[burst])
+
+
+def test_group_by_burst_non_opera():
+    with pytest.raises(ValueError, match="Could not parse burst id"):
+        group_by_burst(["20200101.slc", "20200202.slc"])
+        # A combination should still error
+        group_by_burst(
+            [
+                "20200101.slc",
+                Path("t087_185679_iw1/20180210/t087_185679_iw1_20180210_VV.h5"),
+            ]
+        )
+
+
+def test_setup_output_folder(tmpdir, tiled_file_list):
+    vrt_stack = stack.VRTStack(tiled_file_list, outfile=tmpdir / "stack.vrt")
+    out_file_list = setup_output_folder(vrt_stack, driver="GTiff", dtype=np.complex64)
+    for out_file in out_file_list:
+        assert out_file.exists()
+        assert out_file.suffix == ".tif"
+        assert out_file.parent == tmpdir
+        ds = gdal.Open(str(out_file))
+        assert ds.GetRasterBand(1).DataType == gdal.GDT_CFloat32
+        ds = None
+
+    out_file_list = setup_output_folder(
+        vrt_stack,
+        driver="GTiff",
+        dtype="float32",
+        start_idx=1,
+    )
+    assert len(out_file_list) == len(vrt_stack) - 1
+    for out_file in out_file_list:
+        assert out_file.exists()
+        ds = gdal.Open(str(out_file))
+        assert ds.GetRasterBand(1).DataType == gdal.GDT_Float32
+
+
+@pytest.mark.parametrize(
+    "strides", [{"x": 1, "y": 1}, {"x": 1, "y": 2}, {"x": 2, "y": 3}, {"x": 4, "y": 2}]
+)
+def test_setup_output_folder_strided(tmpdir, tiled_file_list, strides):
+    vrt_stack = stack.VRTStack(tiled_file_list, outfile=tmpdir / "stack.vrt")
+
+    out_file_list = setup_output_folder(
+        vrt_stack, driver="GTiff", dtype=np.complex64, strides=strides
+    )
+    rows, cols = vrt_stack.shape[-2:]
+    for out_file in out_file_list:
+        assert out_file.exists()
+        assert out_file.suffix == ".tif"
+        assert out_file.parent == tmpdir
+
+        ds = gdal.Open(str(out_file))
+        assert ds.GetRasterBand(1).DataType == gdal.GDT_CFloat32
+        assert ds.RasterXSize == cols // strides["x"]
+        assert ds.RasterYSize == rows // strides["y"]
+        ds = None