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

Issue:
Add automatic reference point selection based on highest quality PS
The default reference point at the end can be a PS with
- the lowest value of $D_A$
- the highest amplitude
[Bug]: multiscale_unwrap duplicates unwrapping when tiles=(1,1)

Unwrapping with tiles=(1,1) causes output like this:

```
iter: 0  supplys_left: 9512
iter: 0  supplys_left: 9512
iter: 1  supplys_left: 1750
iter: 1  supplys_left: 1750
iter: 2  supplys_left: 341
iter: 2  supplys_left: 341
```

The duplicated output is due to unwrapping being run twice for the identical "course" and "tiled".




Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index f414abe0..c268812a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,13 +4,16 @@
 - Added `dolphin.timeseries` module with basic functionality:
   - Invert a stack of unwrapped interferograms to a timeseries (using correlation weighting optionally)
   - Estimate a (weighted) linear velocity from a timeseries
+- Added inversion and velocity estimation as options to `DisplacementWorkflow`
 - Create `DatasetStackWriter` protocol, with `BackgroundStackWriter` implementation
 
 **Changed**
 - Rename `GdalWriter` to `BackgroundBlockWriter`
+- Displacement workflow now also creates/returns a stitched, multi-looked version of the amplitude dispersion
 
 **Fixed**
 - `BackgroundRasterWriter`  was not creating the files necessary before writing
+- Allow user to specify more than one type of interferogram in `Network` configuration
 
 # [0.15.3](https://github.com/isce-framework/dolphin/compare/v0.15.2...0.15.3) - 2024-02-27
 
diff --git a/src/dolphin/_types.py b/src/dolphin/_types.py
index 51b0c259..6cbe2596 100644
--- a/src/dolphin/_types.py
+++ b/src/dolphin/_types.py
@@ -75,6 +75,11 @@ class HalfWindow(NamedTuple):
 DateOrDatetime = Union[datetime.datetime, datetime.date]
 
 
+class ReferencePoint(NamedTuple):
+    row: int
+    col: int
+
+
 class TropoModel(str, Enum):
     """Enumeration representing different tropospheric models."""
 
diff --git a/src/dolphin/interferogram.py b/src/dolphin/interferogram.py
index 3575d727..905daad2 100644
--- a/src/dolphin/interferogram.py
+++ b/src/dolphin/interferogram.py
@@ -375,33 +375,44 @@ def __post_init__(
     def _make_ifg_pairs(self) -> list[IfgPairT]:
         """Form interferogram pairs from a list of SLC files sorted by date."""
         assert self.dates is not None
+        ifgs: list[IfgPairT] = []
+        # For each type of possible network, add the ifgs to the list
         if self.indexes is not None:
             # Give the option to select exactly which interferograms to create
-            ifgs = [
-                (self.slc_list[ref_idx], self.slc_list[sec_idx])
-                for ref_idx, sec_idx in self.indexes
-            ]
-        elif self.max_bandwidth is not None:
-            ifgs = Network._limit_by_bandwidth(self.slc_list, self.max_bandwidth)
-        elif self.max_temporal_baseline is not None:
-            ifgs = Network._limit_by_temporal_baseline(
-                self.slc_list,
-                dates=self.dates,
-                max_temporal_baseline=self.max_temporal_baseline,
+            ifgs.extend(
+                [
+                    (self.slc_list[ref_idx], self.slc_list[sec_idx])
+                    for ref_idx, sec_idx in self.indexes
+                ]
             )
-        elif self.reference_idx is not None:
-            ifgs = Network._single_reference_network(self.slc_list, self.reference_idx)
-        else:
+        if self.max_bandwidth is not None:
+            ifgs.extend(Network._limit_by_bandwidth(self.slc_list, self.max_bandwidth))
+        if self.max_temporal_baseline is not None:
+            ifgs.extend(
+                Network._limit_by_temporal_baseline(
+                    self.slc_list,
+                    dates=self.dates,
+                    max_temporal_baseline=self.max_temporal_baseline,
+                )
+            )
+        if self.reference_idx is not None:
+            ifgs.extend(
+                Network._single_reference_network(self.slc_list, self.reference_idx)
+            )
+
+        if self.include_annual:
+            # Add in the annual pairs, then re-sort
+            annual_ifgs = Network._find_annuals(
+                self.slc_list, self.dates, buffer_days=self.annual_buffer_days
+            )
+            ifgs.extend(annual_ifgs)
+
+        if not ifgs:
             msg = "No valid ifg list generation method specified"
             raise ValueError(msg)
 
-        if not self.include_annual:
-            return ifgs
-        # Add in the annual pairs, then re-sort
-        annual_ifgs = Network._find_annuals(
-            self.slc_list, self.dates, buffer_days=self.annual_buffer_days
-        )
-        return sorted(ifgs + annual_ifgs)
+        # Sort and dedupe them
+        return sorted(set(ifgs))
 
     def _create_vrt_ifgs(self) -> list[VRTInterferogram]:
         """Write out a VRTInterferogram for each ifg."""
diff --git a/src/dolphin/io/_readers.py b/src/dolphin/io/_readers.py
index ff9b7d6a..0caa172d 100644
--- a/src/dolphin/io/_readers.py
+++ b/src/dolphin/io/_readers.py
@@ -682,9 +682,6 @@ class VRTStack(StackReader):
     sort_files : bool, optional (default = True)
         Sort the files in `file_list`. Assumes that the naming convention
         will sort the files in increasing time order.
-    nodata_mask_file : pathlib.Path, optional
-        Path to file containing a mask of pixels containing with nodata
-        in every images. Used for skipping the loading of these pixels.
     file_date_fmt : str, optional (default = "%Y%m%d")
         Format string for parsing the dates from the filenames.
         Passed to [opera_utils.get_dates][].
@@ -703,6 +700,7 @@ def __init__(
         fail_on_overwrite: bool = False,
         skip_size_check: bool = False,
         num_threads: int = 1,
+        read_masked: bool = False,
     ):
         if Path(outfile).exists() and write_file:
             if fail_on_overwrite:
@@ -730,6 +728,7 @@ def __init__(
         self.file_list = files
         self.dates = dates
         self.num_threads = num_threads
+        self._read_masked = read_masked
 
         self.outfile = Path(outfile).resolve()
         # Assumes that all files use the same subdataset (if NetCDF)
@@ -740,15 +739,21 @@ def __init__(
 
         # Use the first file in the stack to get size, transform info
         ds = gdal.Open(fspath(self._gdal_file_strings[0]))
+        bnd1 = ds.GetRasterBand(1)
         self.xsize = ds.RasterXSize
         self.ysize = ds.RasterYSize
+        self.nodatavals = []
+        for i in range(1, ds.RasterCount + 1):
+            bnd = ds.GetRasterBand(i)
+            self.nodatavals.append(bnd.GetNoDataValue())
+        self.nodata = self.nodatavals[0]
         # Should be CFloat32
-        self.gdal_dtype = gdal.GetDataTypeName(ds.GetRasterBand(1).DataType)
+        self.gdal_dtype = gdal.GetDataTypeName(bnd1.DataType)
         # Save these for setting at the end
         self.gt = ds.GetGeoTransform()
         self.proj = ds.GetProjection()
         self.srs = ds.GetSpatialRef()
-        ds = None
+        ds = bnd1 = None
         # Save the subset info
 
         self.xoff, self.yoff = 0, 0
@@ -794,6 +799,12 @@ def _write(self):
         ds.SetGeoTransform(self.gt)
         ds.SetProjection(self.proj)
         ds.SetSpatialRef(self.srs)
+        if self.nodata is not None:
+            for i in range(ds.RasterCount):
+                # ds.GetRasterBand(i + 1).SetNoDataValue(self.nodatavals[i])
+                # Force to be the same nodataval for all bands
+                ds.GetRasterBand(i + 1).SetNoDataValue(self.nodata)
+
         ds = None
 
     @property
@@ -899,9 +910,11 @@ def read_stack(
         subsample_factor: int = 1,
         rows: Optional[slice] = None,
         cols: Optional[slice] = None,
-        masked: bool = False,
+        masked: bool | None = None,
     ):
         """Read in the SLC stack."""
+        if masked is None:
+            masked = self._read_masked
         return io.load_gdal(
             self.outfile,
             band=band,
diff --git a/src/dolphin/ps.py b/src/dolphin/ps.py
index cbc2a794..18c5f3df 100644
--- a/src/dolphin/ps.py
+++ b/src/dolphin/ps.py
@@ -261,12 +261,12 @@ def _use_existing_files(
     shutil.copy(existing_amp_mean_file, output_amp_mean_file)
 
 
-def multilook_ps_mask(
+def multilook_ps_files(
     strides: dict[str, int],
     ps_mask_file: Filename,
-    output_file: Optional[Filename] = None,
-) -> Path:
-    """Create a multilooked version of the full-res PS mask.
+    amp_dispersion_file: Filename,
+) -> tuple[Path, Path]:
+    """Create a multilooked version of the full-res PS mask/ampltiude dispersion.
 
     Parameters
     ----------
@@ -274,43 +274,71 @@ def multilook_ps_mask(
         Decimation factor for 'x', 'y'
     ps_mask_file : Filename
         Name of input full-res uint8 PS mask file
-    output_file : Optional[Filename], optional
-        Name of file to save result to.
-        Defaults to same as `ps_mask_file`, but with "_looked" added before suffix.
+    amp_dispersion_file : Filename
+        Name of input full-res float32 amplitude dispersion file
 
     Returns
     -------
-    output_file : Path
+    output_ps_file : Path
+        Multilooked PS mask file
+        Will be same as `ps_mask_file`, but with "_looked" added before suffix.
+    output_amp_disp_file : Path
+        Multilooked amplitude dispersion file
+        Similar naming scheme to `output_ps_file`
 
     """
     if strides == {"x": 1, "y": 1}:
         logger.info("No striding request, skipping multilook.")
-        return Path(ps_mask_file)
-    if output_file is None:
-        ps_suffix = Path(ps_mask_file).suffix
-        out_path = Path(str(ps_mask_file).replace(ps_suffix, f"_looked{ps_suffix}"))
-        logger.info(f"Saving a looked PS mask to {out_path}")
-    else:
-        out_path = Path(output_file)
+        return Path(ps_mask_file), Path(amp_dispersion_file)
+    full_cols, full_rows = io.get_raster_xysize(ps_mask_file)
+    out_rows, out_cols = full_rows // strides["y"], full_cols // strides["x"]
 
-    if Path(out_path).exists():
-        logger.info(f"{out_path} exists, skipping.")
-        return out_path
+    ps_suffix = Path(ps_mask_file).suffix
+    ps_out_path = Path(str(ps_mask_file).replace(ps_suffix, f"_looked{ps_suffix}"))
+    logger.info(f"Saving a looked PS mask to {ps_out_path}")
 
-    ps_mask = io.load_gdal(ps_mask_file, masked=True).astype(bool)
-    full_rows, full_cols = ps_mask.shape
-    ps_mask_looked = utils.take_looks(
-        ps_mask, strides["y"], strides["x"], func_type="any", edge_strategy="pad"
-    )
-    # make sure it's the same size as the MLE result/temp_coh after padding
-    out_rows, out_cols = full_rows // strides["y"], full_cols // strides["x"]
-    ps_mask_looked = ps_mask_looked[:out_rows, :out_cols]
-    ps_mask_looked = ps_mask_looked.astype("uint8").filled(NODATA_VALUES["ps"])
-    io.write_arr(
-        arr=ps_mask_looked,
-        like_filename=ps_mask_file,
-        output_name=out_path,
-        strides=strides,
-        nodata=NODATA_VALUES["ps"],
+    if Path(ps_out_path).exists():
+        logger.info(f"{ps_out_path} exists, skipping.")
+    else:
+        ps_mask = io.load_gdal(ps_mask_file, masked=True).astype(bool)
+        ps_mask_looked = utils.take_looks(
+            ps_mask, strides["y"], strides["x"], func_type="any", edge_strategy="pad"
+        )
+        # make sure it's the same size as the MLE result/temp_coh after padding
+        ps_mask_looked = ps_mask_looked[:out_rows, :out_cols]
+        ps_mask_looked = ps_mask_looked.astype("uint8").filled(NODATA_VALUES["ps"])
+        io.write_arr(
+            arr=ps_mask_looked,
+            like_filename=ps_mask_file,
+            output_name=ps_out_path,
+            strides=strides,
+            nodata=NODATA_VALUES["ps"],
+        )
+
+    amp_disp_suffix = Path(amp_dispersion_file).suffix
+    amp_disp_out_path = Path(
+        str(amp_dispersion_file).replace(amp_disp_suffix, f"_looked{amp_disp_suffix}")
     )
-    return out_path
+    if amp_disp_out_path.exists():
+        logger.info(f"{amp_disp_out_path} exists, skipping.")
+    else:
+        amp_disp = io.load_gdal(amp_dispersion_file, masked=True)
+        # We use `nanmin` assuming that the multilooked PS is using
+        # the strongest PS (the one with the lowest amplitude dispersion)
+        amp_disp_looked = utils.take_looks(
+            amp_disp,
+            strides["y"],
+            strides["x"],
+            func_type="nanmin",
+            edge_strategy="pad",
+        )
+        amp_disp_looked = amp_disp_looked[:out_rows, :out_cols]
+        amp_disp_looked = amp_disp_looked.filled(NODATA_VALUES["amp_dispersion"])
+        io.write_arr(
+            arr=amp_disp_looked,
+            like_filename=amp_dispersion_file,
+            output_name=amp_disp_out_path,
+            strides=strides,
+            nodata=NODATA_VALUES["amp_dispersion"],
+        )
+    return ps_out_path, amp_disp_out_path
diff --git a/src/dolphin/timeseries.py b/src/dolphin/timeseries.py
index e821b0fc..210e4b52 100644
--- a/src/dolphin/timeseries.py
+++ b/src/dolphin/timeseries.py
@@ -2,7 +2,7 @@
 from concurrent.futures import Future, ThreadPoolExecutor
 from pathlib import Path
 from tempfile import NamedTemporaryFile
-from typing import NamedTuple, Protocol, Sequence, TypeVar
+from typing import Callable, Protocol, Sequence, TypeVar
 
 import jax.numpy as jnp
 import numpy as np
@@ -12,17 +12,19 @@
 from tqdm.auto import tqdm
 
 from dolphin import DateOrDatetime, io
-from dolphin._types import PathOrStr
+from dolphin._overviews import ImageType, create_overviews
+from dolphin._types import PathOrStr, ReferencePoint
 from dolphin.utils import DummyProcessPoolExecutor, flatten, format_dates
 
 __all__ = [
-    "weighted_lstsq_single",
     "invert_stack",
     "get_incidence_matrix",
     "estimate_velocity",
     "create_velocity",
     "create_temporal_average",
     "invert_unw_network",
+    "correlation_to_variance",
+    "select_reference_point",
 ]
 
 T = TypeVar("T")
@@ -285,12 +287,11 @@ def process_blocks(
     block_shape: tuple[int, int] = (512, 512),
     num_threads: int = 5,
 ):
-    """Perform block-wise processing over blocks in `readers`, writing to `writer`."""
-    # reader = RasterStackReader.from_file_list(file_list=file_list)
-    # writer = io.GdalStackWriter(
-    #     file_list=output_files, like_filename=like_filename or file_list[0]
-    # )
-    # shape = io.get_raster_xysize(file_list[0])[::-1]
+    """Perform block-wise processing over blocks in `readers`, writing to `writer`.
+
+    Extracts the common functionality from several routines to read and process
+    a stack of rasters in parallel.
+    """
     shape = readers[0].shape[-2:]
     slices = list(io.iter_blocks(shape, block_shape=block_shape))
 
@@ -312,21 +313,82 @@ def write_callback(fut: Future):
 def create_velocity(
     unw_file_list: Sequence[PathOrStr],
     output_file: PathOrStr,
+    reference: ReferencePoint,
     date_list: Sequence[DateOrDatetime] | None = None,
     cor_file_list: Sequence[PathOrStr] | None = None,
     cor_threshold: float = 0.2,
     block_shape: tuple[int, int] = (512, 512),
     num_threads: int = 5,
+    add_overviews: bool = True,
 ) -> None:
-    """Perform pixel-wise (weighted) linear regression to estimate velocity."""
+    """Perform pixel-wise (weighted) linear regression to estimate velocity.
+
+    Parameters
+    ----------
+    unw_file_list : Sequence[PathOrStr]
+        List of unwrapped phase files.
+    output_file : PathOrStr
+        The output file to save the velocity to.
+    reference : ReferencePoint
+        The (row, col) to use as reference before fitting the velocity.
+        This point will be subtracted from all other points before solving.
+    date_list : Sequence[DateOrDatetime], optional
+        List of dates corresponding to the unwrapped phase files.
+        If not provided, will be parsed from filenames in `unw_file_list`.
+    cor_file_list : Sequence[PathOrStr], optional
+        List of correlation files to use for weighting the velocity estimation.
+        If not provided, all weights are set to 1.
+    cor_threshold : float, optional
+        The correlation threshold to use for weighting the velocity estimation.
+        Default is 0.2.
+    block_shape : tuple[int, int], optional
+        The shape of the blocks to process in parallel.
+        Default is (512, 512)
+    num_threads : int, optional
+        The parallel blocks to process at once.
+        Default is 5.
+    add_overviews : bool, optional
+        If True, creates overviews of the new velocity raster.
+        Default is True.
+
+    """
+    if Path(output_file).exists():
+        logger.info(f"Output file {output_file} already exists, skipping velocity")
+        return
+
     if date_list is None:
         date_list = [get_dates(f)[1] for f in unw_file_list]
-        # TODO: need the relative one?
     x_arr = datetime_to_float(date_list)
 
+    # Set up the input data readers
+    out_dir = Path(output_file).parent
+    unw_reader = io.VRTStack(
+        file_list=unw_file_list,
+        outfile=out_dir / "velocity_inputs.vrt",
+        skip_size_check=True,
+    )
+    if cor_file_list is not None:
+        if len(cor_file_list) != len(unw_file_list):
+            msg = "Mismatch in number of input files provided:"
+            msg += f"{len(cor_file_list) = }, but {len(unw_file_list) = }"
+            raise ValueError(msg)
+
+        cor_reader = io.VRTStack(
+            file_list=cor_file_list,
+            outfile=out_dir / "cor_inputs.vrt",
+            skip_size_check=True,
+        )
+    else:
+        cor_reader = None
+
+    # Read in the reference point
+    ref_row, ref_col = reference
+    ref_data = unw_reader[:, ref_row, ref_col].reshape(-1, 1, 1)
+
     def read_and_fit(
         readers: Sequence[io.StackReader], rows: slice, cols: slice
     ) -> tuple[slice, slice, np.ndarray]:
+        # Only use the cor_reader if it's the same shape as the unw_reader
         if len(readers) == 2:
             unw_reader, cor_reader = readers
             unw_stack = unw_reader[:, rows, cols]
@@ -335,45 +397,44 @@ def read_and_fit(
         else:
             unw_stack = readers[0][:, rows, cols]
             weights = np.ones_like(unw_stack)
+        # Reference the data
+        unw_stack = unw_stack - ref_data
         # Fit a line to each pixel with weighted least squares
         return (
             rows,
             cols,
-            estimate_velocity(x_arr=x_arr, stack=unw_stack, weight_stack=weights),
+            estimate_velocity(x_arr=x_arr, unw_stack=unw_stack, weight_stack=weights),
         )
 
-    # Note: For some reason, the `RasterStackReader` is much slower than the VRT:
+    # Note: For some reason, the `RasterStackReader` is much slower than the VRT
+    # for files on S3:
     # ~300 files takes >2 min to open, >2 min to read each block
     # VRTStack seems to take ~30 secs to open, 1 min to read
     # Very possible there's a tuning param/rasterio config to fix, but not sure.
-    # unw_reader = io.RasterStackReader.from_file_list(file_list=unw_file_list)
-    # cor_reader = io.RasterStackReader.from_file_list(file_list=cor_file_list)
-    with NamedTemporaryFile(mode="w", suffix=".vrt") as f1, NamedTemporaryFile(
-        mode="w", suffix=".vrt"
-    ) as f2:
-        unw_reader = io.VRTStack(
-            file_list=unw_file_list, outfile=f1.name, skip_size_check=True
-        )
-        if cor_file_list is not None:
-            cor_reader = io.VRTStack(
-                file_list=cor_file_list, outfile=f2.name, skip_size_check=True
-            )
-            readers = [unw_reader, cor_reader]
-        else:
-            readers = [unw_reader]
+    readers = [unw_reader]
+    if cor_reader is not None:
+        readers.append(cor_reader)
 
-        writer = io.BackgroundRasterWriter(output_file, like_filename=unw_file_list[0])
-        process_blocks(
-            # file_list=unw_file_list,
-            # output_files=[output_file],
-            readers=readers,
-            writer=writer,
-            func=read_and_fit,
-            block_shape=block_shape,
-            num_threads=num_threads,
-        )
+    writer = io.BackgroundRasterWriter(output_file, like_filename=unw_file_list[0])
+    process_blocks(
+        readers=readers,
+        writer=writer,
+        func=read_and_fit,
+        block_shape=block_shape,
+        num_threads=num_threads,
+    )
 
     writer.notify_finished()
+    if add_overviews:
+        logger.info("Creating overviews for velocity image")
+        create_overviews([output_file])
+
+
+class AverageFunc(Protocol):
+    """Protocol for temporally averaging a block of data."""
+
+    def __call__(self, ArrayLike, axis: int) -> ArrayLike:
+        ...
 
 
 def create_temporal_average(
@@ -381,21 +442,48 @@ def create_temporal_average(
     output_file: PathOrStr,
     block_shape: tuple[int, int] = (512, 512),
     num_threads: int = 5,
+    average_func: Callable[[ArrayLike, int], np.ndarray] = np.nanmean,
+    read_masked: bool = False,
 ) -> None:
-    """Average all images in `reader` to create a 2D image in `output_file`."""
+    """Average all images in `reader` to create a 2D image in `output_file`.
+
+    Parameters
+    ----------
+    file_list : Sequence[PathOrStr]
+        List of files to average
+    output_file : PathOrStr
+        The output file to save the average to
+    block_shape : tuple[int, int], optional
+        The shape of the blocks to process in parallel.
+        Default is (512, 512)
+    num_threads : int, optional
+        The parallel blocks to process at once.
+        Default is 5.
+    average_func : Callable[[ArrayLike, int], np.ndarray], optional
+        The function to use to average the images.
+        Default is `np.nanmean`, which calls `np.nanmean(arr, axis=0)` on each block.
+    read_masked : bool, optional
+        If True, reads the data as a masked array based on the rasters' nodata values.
+        Default is False.
+
+    """
 
     def read_and_average(
         readers: Sequence[io.StackReader], rows: slice, cols: slice
     ) -> tuple[slice, slice, np.ndarray]:
-        return rows, cols, np.nanmean(readers[0][:, rows, cols], axis=0)
+        chunk = readers[0][:, rows, cols]
+        return rows, cols, average_func(chunk, 0)
 
     writer = io.BackgroundRasterWriter(output_file, like_filename=file_list[0])
     with NamedTemporaryFile(mode="w", suffix=".vrt") as f:
-        reader = io.VRTStack(file_list=file_list, outfile=f.name, skip_size_check=True)
+        reader = io.VRTStack(
+            file_list=file_list,
+            outfile=f.name,
+            skip_size_check=True,
+            read_masked=read_masked,
+        )
 
         process_blocks(
-            # file_list=file_list,
-            # output_files=[output_file],
             readers=[reader],
             writer=writer,
             func=read_and_average,
@@ -406,20 +494,17 @@ def read_and_average(
     writer.notify_finished()
 
 
-class ReferencePoint(NamedTuple):
-    row: int
-    col: int
-
-
 def invert_unw_network(
     unw_file_list: Sequence[PathOrStr],
     reference: ReferencePoint,
     output_dir: PathOrStr,
     cor_file_list: Sequence[PathOrStr] | None = None,
     cor_threshold: float = 0.2,
+    n_cor_looks: int = 1,
     ifg_date_pairs: Sequence[Sequence[DateOrDatetime]] | None = None,
     block_shape: tuple[int, int] = (512, 512),
     num_threads: int = 5,
+    add_overviews: bool = True,
 ) -> list[Path]:
     """Perform pixel-wise inversion of unwrapped network to get phase per date.
 
@@ -442,8 +527,17 @@ def invert_unw_network(
         List of correlation files to use for weighting the inversion
     cor_threshold : float, optional
         The correlation threshold to use for weighting the inversion
+        Default is 0.2
+    n_cor_looks : int, optional
+        The number of looks used to form the input correlation data, used
+        to convert correlation to phase variance.
+        Default is 1.
     num_threads : int, optional
         The parallel blocks to process at once.
+        Default is 5.
+    add_overviews : bool, optional
+        If True, creates overviews of the new unwrapped phase rasters.
+        Default is True.
 
     Returns
     -------
@@ -462,14 +556,18 @@ def invert_unw_network(
             "Each item in `ifg_date_pairs` must be a sequence of length 2"
         ) from e
 
-    sar_dates = sorted(set(flatten(ifg_tuples)))
-    A = get_incidence_matrix(ifg_pairs=ifg_tuples, sar_idxs=sar_dates)
     # Make the names of the output files from the SAR dates to solve for
+    sar_dates = sorted(set(flatten(ifg_tuples)))
     ref_date = sar_dates[0]
     suffix = ".tif"
     out_paths = [
         Path(output_dir) / (format_dates(ref_date, d) + suffix) for d in sar_dates
     ]
+    if all(p.exists() for p in out_paths):
+        logger.info("All output files already exist, skipping inversion")
+        return out_paths
+
+    A = get_incidence_matrix(ifg_pairs=ifg_tuples, sar_idxs=sar_dates)
 
     out_vrt_name = Path(output_dir) / "unw_network.vrt"
     unw_reader = io.VRTStack(
@@ -487,8 +585,9 @@ def read_and_solve(
         if len(readers) == 2:
             unw_reader, cor_reader = readers
             stack = unw_reader[:, rows, cols]
-            weights = cor_reader[:, rows, cols]
-            weights[weights < cor_threshold] = 0
+            cor = cor_reader[:, rows, cols]
+            weights = correlation_to_variance(cor, n_cor_looks)
+            weights[cor < cor_threshold] = 0
         else:
             stack = readers[0][:, rows, cols]
             weights = np.ones_like(stack)
@@ -520,5 +619,103 @@ def read_and_solve(
         num_threads=num_threads,
     )
     writer.notify_finished()
-    # Return the output files
+
+    if add_overviews:
+        logger.info("Creating overviews for unwrapped images")
+        create_overviews(out_paths, image_type=ImageType.UNWRAPPED)
     return out_paths
+
+
+def correlation_to_variance(correlation: ArrayLike, nlooks: int) -> Array:
+    r"""Convert interferometric correlation to phase variance.
+
+    Uses the CRLB formula from Rodriguez, 1992 [1]_ to get the phase variance,
+    \sigma_{\phi}^2:
+
+    \[
+        \sigma_{\phi}^{2} = \frac{1}{2N_{L}} \frac{1 - \gamma^{2}}{\gamma^{2}}
+    \]
+
+    where \gamma is the correlation and N_L is the effective number of looks.
+
+    References
+    ----------
+    .. [1] Rodriguez, E., and J. M. Martin. "Theory and design of interferometric
+    synthetic aperture radars." IEEE Proceedings of Radar and Signal Processing.
+    Vol. 139. No. 2. IET Digital Library, 1992.
+
+    """
+    return (1 - correlation**2) / (2 * nlooks * correlation**2 + 1e-6)
+
+
+def select_reference_point(
+    ccl_file_list: Sequence[PathOrStr],
+    amp_dispersion_file: PathOrStr,
+    output_dir: Path,
+    block_shape: tuple[int, int] = (512, 512),
+    num_threads: int = 5,
+) -> ReferencePoint:
+    """Automatically select a reference point for a stack of unwrapped interferograms.
+
+    Uses the amplitude dispersion and connected component labels, the point is selected
+    which
+
+    1. is within intersection of all nonzero connected component labels (always valid)
+    2. has the lowest amplitude dispersion
+    """
+    conncomp_intersection_file = Path(output_dir) / "conncomp_intersection.tif"
+
+    def intersect_conncomp(arr: np.ma.MaskedArray, axis: int) -> np.ndarray:
+        # Track where input is nodata
+        any_masked = np.any(arr.mask, axis=axis)
+        # Get the logical AND of all nonzero conncomp labels
+        fillval = arr.fill_value
+        is_valid_conncomp = arr.filled(0) > 0
+        all_are_valid = np.all(is_valid_conncomp, axis=axis).astype(arr.dtype)
+        # Reset nodata
+        all_are_valid[any_masked] = fillval
+        return all_are_valid
+
+    if not conncomp_intersection_file.exists():
+        logger.info("Creating intersection of connected components")
+        create_temporal_average(
+            file_list=ccl_file_list,
+            output_file=conncomp_intersection_file,
+            block_shape=block_shape,
+            num_threads=num_threads,
+            average_func=intersect_conncomp,
+            read_masked=True,
+        )
+
+    logger.info("Selecting reference point")
+    conncomp_intersection = io.load_gdal(conncomp_intersection_file, masked=True)
+
+    # Find the largest conncomp region in the intersection
+    isin_largest_conncomp = get_largest_conncomp(conncomp_intersection)
+
+    amp_dispersion = io.load_gdal(amp_dispersion_file, masked=True)
+    # Mask out where the conncomps aren't equal to the largest
+    amp_dispersion.mask = amp_dispersion.mask | (~isin_largest_conncomp)
+
+    # Pick the (unmasked) point with the lowest amplitude dispersion
+    ref_row, ref_col = np.unravel_index(np.argmin(amp_dispersion), amp_dispersion.shape)
+
+    # Cast to `int` to avoid having `np.int64` types
+    return ReferencePoint(int(ref_row), int(ref_col))
+
+
+def get_largest_conncomp(conncomp_intersection: ArrayLike) -> np.ndarray:
+    """Get a max for the largest nonzero connected component in the intersection."""
+    from scipy import ndimage
+
+    label, nlabels = ndimage.label(
+        conncomp_intersection.filled(0), structure=np.ones((3, 3))
+    )
+    logger.info("Found %d connected components in intersection", nlabels)
+    # Find the label with the most pixels using bincount
+    label_counts = np.bincount(label.ravel())
+    # (ignore the 0 label)
+    largest_idx = np.argmax(label_counts[1:]) + 1
+
+    # Make a mask of the largest conncomp
+    return label == largest_idx
diff --git a/src/dolphin/workflows/config/_common.py b/src/dolphin/workflows/config/_common.py
index f1e4a4ef..0c65ac96 100644
--- a/src/dolphin/workflows/config/_common.py
+++ b/src/dolphin/workflows/config/_common.py
@@ -32,6 +32,7 @@
     "PsOptions",
     "PhaseLinkingOptions",
     "InterferogramNetwork",
+    "TimeseriesOptions",
     "UnwrapOptions",
 ]
 
@@ -214,6 +215,39 @@ def _to_tuple(cls, v):
         return v
 
 
+class TimeseriesOptions(BaseModel, extra="forbid"):
+    """Options for inversion/time series fitting."""
+
+    _directory: Path = PrivateAttr(Path("timeseries"))
+    _velocity_file: Path = PrivateAttr(Path("timeseries/velocity.tif"))
+    run_inversion: bool = Field(
+        True,
+        description=(
+            "Whether to run the inversion step after unwrapping, if more than "
+            " a single-reference network is used."
+        ),
+    )
+    reference_point: Optional[tuple[int, int]] = Field(
+        None,
+        description=(
+            "Reference point (row, col) used if performing a time series inversion. "
+            "If not provided, a point will be selected from a consistent connected "
+            "component with low amplitude dispersion."
+        ),
+    )
+
+    run_velocity: bool = Field(
+        True,
+        description="Run the velocity estimation from the phase time series.",
+    )
+    correlation_threshold: float = Field(
+        0.2,
+        description="Pixels with correlation below this value will be masked out.",
+        ge=0.0,
+        le=1.0,
+    )
+
+
 class WorkerSettings(BaseModel, extra="forbid"):
     """Settings for controlling CPU/GPU settings and parallelism."""
 
diff --git a/src/dolphin/workflows/config/_displacement.py b/src/dolphin/workflows/config/_displacement.py
index 2990af0c..494cf6e4 100644
--- a/src/dolphin/workflows/config/_displacement.py
+++ b/src/dolphin/workflows/config/_displacement.py
@@ -22,6 +22,7 @@
     OutputOptions,
     PhaseLinkingOptions,
     PsOptions,
+    TimeseriesOptions,
     UnwrapOptions,
     WorkflowBase,
     _read_file_list_or_glob,
@@ -133,6 +134,7 @@ class DisplacementWorkflow(WorkflowBase):
         default_factory=InterferogramNetwork
     )
     unwrap_options: UnwrapOptions = Field(default_factory=UnwrapOptions)
+    timeseries_options: TimeseriesOptions = Field(default_factory=TimeseriesOptions)
     correction_options: CorrectionOptions = Field(default_factory=CorrectionOptions)
 
     # internal helpers
@@ -197,6 +199,7 @@ def model_post_init(self, __context: Any) -> None:
             "phase_linking",
             "interferogram_network",
             "unwrap_options",
+            "timeseries_options",
         ]:
             opts = getattr(self, step)
             if isinstance(opts, dict):
@@ -218,6 +221,7 @@ def model_post_init(self, __context: Any) -> None:
             self.phase_linking._directory,
             self.interferogram_network._directory,
             self.unwrap_options._directory,
+            self.timeseries_options._directory,
         ]
         # Add the output PS files we'll create to the `PS` directory, making
         # sure they're inside the work directory
@@ -225,3 +229,7 @@ def model_post_init(self, __context: Any) -> None:
         ps_opts._amp_dispersion_file = work_dir / ps_opts._amp_dispersion_file
         ps_opts._amp_mean_file = work_dir / ps_opts._amp_mean_file
         ps_opts._output_file = work_dir / ps_opts._output_file
+
+        self.timeseries_options._velocity_file = (
+            work_dir / self.timeseries_options._velocity_file
+        )
diff --git a/src/dolphin/workflows/displacement.py b/src/dolphin/workflows/displacement.py
index b3ba3580..0596b3ad 100755
--- a/src/dolphin/workflows/displacement.py
+++ b/src/dolphin/workflows/displacement.py
@@ -1,6 +1,7 @@
 #!/usr/bin/env python
 from __future__ import annotations
 
+import contextlib
 import multiprocessing as mp
 from collections import defaultdict
 from concurrent.futures import ProcessPoolExecutor
@@ -10,24 +11,16 @@
 from pprint import pformat
 from typing import Mapping, NamedTuple, Sequence
 
-from opera_utils import group_by_burst, group_by_date
+from opera_utils import get_dates, group_by_burst, group_by_date
 from tqdm.auto import tqdm
 
-from dolphin import __version__
+from dolphin import __version__, io, timeseries, utils
 from dolphin._log import get_log, log_runtime
 from dolphin.atmosphere import estimate_ionospheric_delay, estimate_tropospheric_delay
-from dolphin.io import get_raster_bounds, get_raster_crs
-from dolphin.utils import (
-    DummyProcessPoolExecutor,
-    disable_gpu,
-    get_max_memory_usage,
-    prepare_geometry,
-    set_num_threads,
-)
 
 from . import stitching_bursts, unwrapping, wrapped_phase
 from ._utils import _create_burst_cfg, _remove_dir_if_empty
-from .config import DisplacementWorkflow
+from .config import DisplacementWorkflow, TimeseriesOptions
 
 logger = get_log(__name__)
 
@@ -40,6 +33,7 @@ class OutputPaths(NamedTuple):
     stitched_cor_paths: list[Path]
     stitched_temp_coh_file: Path
     stitched_ps_file: Path
+    stitched_amp_dispersion_file: Path
     unwrapped_paths: list[Path] | None
     conncomp_paths: list[Path] | None
     tropospheric_corrections: list[Path] | None
@@ -67,8 +61,8 @@ def run(
     logger.debug(pformat(cfg.model_dump()))
 
     if not cfg.worker_settings.gpu_enabled:
-        disable_gpu()
-    set_num_threads(cfg.worker_settings.threads_per_worker)
+        utils.disable_gpu()
+    utils.set_num_threads(cfg.worker_settings.threads_per_worker)
 
     try:
         grouped_slc_files = group_by_burst(cfg.cslc_file_list)
@@ -140,6 +134,7 @@ def run(
     ifg_file_list: list[Path] = []
     temp_coh_file_list: list[Path] = []
     ps_file_list: list[Path] = []
+    amp_dispersion_file_list: list[Path] = []
     # The comp_slc tracking object is a dict, since we'll need to organize
     # multiple comp slcs by burst (they'll have the same filename)
     comp_slc_dict: dict[str, list[Path]] = {}
@@ -148,7 +143,7 @@ def run(
     Executor = (
         ProcessPoolExecutor
         if cfg.worker_settings.n_parallel_bursts > 1
-        else DummyProcessPoolExecutor
+        else utils.DummyProcessPoolExecutor
     )
     mw = cfg.worker_settings.n_parallel_bursts
     ctx = mp.get_context("spawn")
@@ -173,11 +168,12 @@ def run(
         for fut in fut_to_burst:
             burst = fut_to_burst[fut]
 
-            cur_ifg_list, comp_slcs, temp_coh, ps_file = fut.result()
+            cur_ifg_list, comp_slcs, temp_coh, ps_file, amp_disp_file = fut.result()
             ifg_file_list.extend(cur_ifg_list)
             comp_slc_dict[burst] = comp_slcs
             temp_coh_file_list.append(temp_coh)
             ps_file_list.append(ps_file)
+            amp_dispersion_file_list.append(amp_disp_file)
 
     # ###################################
     # 2. Stitch burst-wise interferograms
@@ -192,10 +188,12 @@ def run(
         stitched_cor_paths,
         stitched_temp_coh_file,
         stitched_ps_file,
+        stitched_amp_dispersion_file,
     ) = stitching_bursts.run(
         ifg_file_list=ifg_file_list,
         temp_coh_file_list=temp_coh_file_list,
         ps_file_list=ps_file_list,
+        amp_dispersion_list=amp_dispersion_file_list,
         stitched_ifg_dir=cfg.interferogram_network._directory,
         output_options=cfg.output_options,
         file_date_fmt=cfg.input_options.cslc_date_fmt,
@@ -214,6 +212,7 @@ def run(
             stitched_cor_paths=stitched_cor_paths,
             stitched_temp_coh_file=stitched_temp_coh_file,
             stitched_ps_file=stitched_ps_file,
+            stitched_amp_dispersion_file=stitched_amp_dispersion_file,
             unwrapped_paths=None,
             conncomp_paths=None,
             tropospheric_corrections=None,
@@ -231,7 +230,25 @@ def run(
     )
 
     # ##############################################
-    # 4. Estimate corrections for each interferogram
+    # 4. If a network was unwrapped, invert it
+    # ##############################################
+
+    ts_opts = cfg.timeseries_options
+    if ts_opts.run_inversion or ts_opts.run_velocity:
+        inverted_phase_paths = run_timeseries(
+            ts_opts=ts_opts,
+            unwrapped_paths=unwrapped_paths,
+            conncomp_paths=conncomp_paths,
+            cor_paths=stitched_cor_paths,
+            stitched_amp_dispersion_file=stitched_amp_dispersion_file,
+            # TODO: do i care to configure block shape, or num threads from somewhere?
+            # num_threads=cfg.worker_settings....?
+        )
+    else:
+        inverted_phase_paths = unwrapped_paths
+
+    # ##############################################
+    # 5. Estimate corrections for each interferogram
     # ##############################################
     tropo_paths: list[Path] | None = None
     iono_paths: list[Path] | None = None
@@ -245,10 +262,10 @@ def run(
         # Prepare frame geometry files
         geometry_dir = out_dir / "geometry"
         geometry_dir.mkdir(exist_ok=True)
-        crs = get_raster_crs(ifg_filenames[0])
+        crs = io.get_raster_crs(ifg_filenames[0])
         epsg = crs.to_epsg()
-        out_bounds = get_raster_bounds(ifg_filenames[0])
-        frame_geometry_files = prepare_geometry(
+        out_bounds = io.get_raster_bounds(ifg_filenames[0])
+        frame_geometry_files = utils.prepare_geometry(
             geometry_dir=geometry_dir,
             geo_files=cfg.correction_options.geometry_files,
             matching_file=ifg_filenames[0],
@@ -302,7 +319,9 @@ def run(
         stitched_cor_paths=stitched_cor_paths,
         stitched_temp_coh_file=stitched_temp_coh_file,
         stitched_ps_file=stitched_ps_file,
-        unwrapped_paths=unwrapped_paths,
+        stitched_amp_dispersion_file=stitched_amp_dispersion_file,
+        # unwrapped_paths=unwrapped_paths,
+        unwrapped_paths=inverted_phase_paths,
         conncomp_paths=conncomp_paths,
         tropospheric_corrections=tropo_paths,
         ionospheric_corrections=iono_paths,
@@ -311,7 +330,82 @@ def run(
 
 def _print_summary(cfg):
     """Print the maximum memory usage and version info."""
-    max_mem = get_max_memory_usage(units="GB")
+    max_mem = utils.get_max_memory_usage(units="GB")
     logger.info(f"Maximum memory usage: {max_mem:.2f} GB")
     logger.info(f"Config file dolphin version: {cfg._dolphin_version}")
     logger.info(f"Current running dolphin version: {__version__}")
+
+
+def run_timeseries(
+    ts_opts: TimeseriesOptions,
+    unwrapped_paths: Sequence[Path],
+    conncomp_paths: Sequence[Path],
+    cor_paths: Sequence[Path],
+    stitched_amp_dispersion_file: Path,
+    num_threads: int = 5,
+) -> list[Path]:
+    """Invert the unwrapped interferograms, estimate timeseries and phase velocity."""
+    output_path = ts_opts._directory
+    output_path.mkdir(exist_ok=True, parents=True)
+
+    # First we find the reference point for the unwrapped interferograms
+    reference = timeseries.select_reference_point(
+        conncomp_paths,
+        stitched_amp_dispersion_file,
+        output_dir=output_path,
+    )
+
+    ifg_date_pairs = [get_dates(f) for f in unwrapped_paths]
+    sar_dates = sorted(set(utils.flatten(ifg_date_pairs)))
+    # if we did single-reference interferograms, for `n` sar dates, we will only have
+    # `n-1` interferograms. Any more than n-1 ifgs means we need to invert
+    needs_inversion = len(unwrapped_paths) > len(sar_dates) - 1
+    # check if we even need to invert, or if it was single reference
+    inverted_phase_paths: list[Path] = []
+    if needs_inversion:
+        logger.info("Selecting a reference point for unwrapped interferograms")
+
+        logger.info("Inverting network of %s unwrapped ifgs", len(unwrapped_paths))
+        inverted_phase_paths = timeseries.invert_unw_network(
+            unw_file_list=unwrapped_paths,
+            reference=reference,
+            output_dir=output_path,
+        )
+    else:
+        logger.info(
+            "Skipping inversion step: only single reference interferograms exist."
+        )
+        # Symlink the unwrapped paths to `timeseries/`
+        for p in unwrapped_paths:
+            target = output_path / p.name
+            with contextlib.suppress(FileExistsError):
+                target.symlink_to(p)
+            inverted_phase_paths.append(target)
+        # Make extra "0" raster so that the number of rasters matches len(sar_dates)
+        ref_raster = output_path / (
+            utils.format_dates(sar_dates[0], sar_dates[0]) + ".tif"
+        )
+        io.write_arr(
+            arr=None, output_name=ref_raster, like_filename=inverted_phase_paths[0]
+        )
+        inverted_phase_paths.append(ref_raster)
+
+    if ts_opts.run_velocity:
+        #  We can't pass the correlations after an inversion- the numbers don't match
+        # TODO:
+        # Is there a better weighting then?
+        cor_file_list = (
+            cor_paths if len(cor_paths) == len(inverted_phase_paths) else None
+        )
+        logger.info("Estimating phase velocity")
+        timeseries.create_velocity(
+            unw_file_list=inverted_phase_paths,
+            output_file=ts_opts._velocity_file,
+            reference=reference,
+            date_list=sar_dates,
+            cor_file_list=cor_file_list,
+            cor_threshold=ts_opts.correlation_threshold,
+            num_threads=num_threads,
+        )
+
+    return inverted_phase_paths
diff --git a/src/dolphin/workflows/ps.py b/src/dolphin/workflows/ps.py
index c21e1c09..47f7e81d 100644
--- a/src/dolphin/workflows/ps.py
+++ b/src/dolphin/workflows/ps.py
@@ -100,10 +100,13 @@ def run(
     # Save a looked version of the PS mask too
     strides = cfg.output_options.strides
     if compute_looked:
-        ps_looked_file = dolphin.ps.multilook_ps_mask(
-            strides=strides, ps_mask_file=cfg.ps_options._output_file
+        ps_looked_file, amp_disp_looked = dolphin.ps.multilook_ps_files(
+            strides=strides,
+            ps_mask_file=cfg.ps_options._output_file,
+            amp_dispersion_file=cfg.ps_options._amp_dispersion_file,
         )
         output_file_list.append(ps_looked_file)
+        output_file_list.append(amp_disp_looked)
 
     # Print the maximum memory usage for each worker
     max_mem = get_max_memory_usage(units="GB")
diff --git a/src/dolphin/workflows/stitching_bursts.py b/src/dolphin/workflows/stitching_bursts.py
index 87e2b1e4..63f324f5 100644
--- a/src/dolphin/workflows/stitching_bursts.py
+++ b/src/dolphin/workflows/stitching_bursts.py
@@ -20,11 +20,12 @@ def run(
     ifg_file_list: Sequence[Path],
     temp_coh_file_list: Sequence[Path],
     ps_file_list: Sequence[Path],
+    amp_dispersion_list: Sequence[Path],
     stitched_ifg_dir: Path,
     output_options: OutputOptions,
     file_date_fmt: str = "%Y%m%d",
     corr_window_size: tuple[int, int] = (11, 11),
-) -> tuple[list[Path], list[Path], Path, Path]:
+) -> tuple[list[Path], list[Path], Path, Path, Path]:
     """Run the displacement workflow on a stack of SLCs.
 
     Parameters
@@ -35,6 +36,8 @@ def run(
         Sequence of paths to the burst-wise temporal coherence files.
     ps_file_list : Sequence[Path]
         Sequence of paths to the (looked) burst-wise ps mask files.
+    amp_dispersion_list : Sequence[Path]
+        Sequence of paths to the (looked) burst-wise amplitude dispersion files.
     stitched_ifg_dir : Path
         Location to store the output stitched ifgs and correlations
     output_options : OutputOptions
@@ -57,6 +60,8 @@ def run(
         Path to temporal correlation file created.
     stitched_ps_file : Path
         Path to ps mask file created.
+    stitched_amp_disp_file : Path
+        Path to amplitude dispersion file created.
 
     """
     stitched_ifg_dir.mkdir(exist_ok=True, parents=True)
@@ -101,16 +106,28 @@ def run(
         out_bounds_epsg=output_options.bounds_epsg,
     )
 
+    # Stitch the amp dispersion files
+    stitched_amp_disp_file = stitched_ifg_dir / "amp_dispersion_looked.tif"
+    stitching.merge_images(
+        amp_dispersion_list,
+        outfile=stitched_amp_disp_file,
+        driver="GTiff",
+        out_bounds=out_bounds,
+        out_bounds_epsg=output_options.bounds_epsg,
+    )
+
     if output_options.add_overviews:
         logger.info("Creating overviews for stitched images")
         create_overviews(stitched_ifg_paths, image_type=ImageType.INTERFEROGRAM)
         create_overviews(interferometric_corr_paths, image_type=ImageType.CORRELATION)
         create_image_overviews(stitched_ps_file, image_type=ImageType.PS)
         create_image_overviews(stitched_temp_coh_file, image_type=ImageType.CORRELATION)
+        create_image_overviews(stitched_amp_disp_file, image_type=ImageType.CORRELATION)
 
     return (
         stitched_ifg_paths,
         interferometric_corr_paths,
         stitched_temp_coh_file,
         stitched_ps_file,
+        stitched_amp_disp_file,
     )
diff --git a/src/dolphin/workflows/wrapped_phase.py b/src/dolphin/workflows/wrapped_phase.py
index a7d51c56..411c1448 100644
--- a/src/dolphin/workflows/wrapped_phase.py
+++ b/src/dolphin/workflows/wrapped_phase.py
@@ -18,7 +18,7 @@
 @log_runtime
 def run(
     cfg: DisplacementWorkflow, debug: bool = False, tqdm_kwargs=None
-) -> tuple[list[Path], list[Path], Path, Path]:
+) -> tuple[list[Path], list[Path], Path, Path, Path]:
     """Run the displacement workflow on a stack of SLCs.
 
     Parameters
@@ -44,6 +44,8 @@ def run(
         In the case of sequential phase linking, this is the average of all ministacks.
     ps_looked_file : Path
         The multilooked boolean persistent scatterer file.
+    amp_disp_looked_file : Path
+        The multilooked amplitude dispersion file.
 
     """
     if tqdm_kwargs is None:
@@ -95,8 +97,10 @@ def run(
 
     # Save a looked version of the PS mask too
     strides = cfg.output_options.strides
-    ps_looked_file = ps.multilook_ps_mask(
-        strides=strides, ps_mask_file=cfg.ps_options._output_file
+    ps_looked_file, amp_disp_looked_file = ps.multilook_ps_files(
+        strides=strides,
+        ps_mask_file=cfg.ps_options._output_file,
+        amp_dispersion_file=cfg.ps_options._amp_dispersion_file,
     )
 
     # #########################
@@ -178,16 +182,28 @@ def run(
     # ###################################################
 
     ifg_network = cfg.interferogram_network
-    existing_ifgs = list(ifg_network._directory.glob("*.int.*"))
+    existing_ifgs = list(ifg_network._directory.glob("*.int.vrt"))
     if len(existing_ifgs) > 0:
         logger.info(f"Skipping interferogram step, {len(existing_ifgs)} exists")
-        return existing_ifgs, comp_slc_list, temp_coh_file, ps_looked_file
+        return (
+            existing_ifgs,
+            comp_slc_list,
+            temp_coh_file,
+            ps_looked_file,
+            amp_disp_looked_file,
+        )
 
     logger.info(f"Creating virtual interferograms from {len(phase_linked_slcs)} files")
     ifg_file_list = create_ifgs(
         ifg_network, phase_linked_slcs, any(is_compressed), reference_date
     )
-    return ifg_file_list, comp_slc_list, temp_coh_file, ps_looked_file
+    return (
+        ifg_file_list,
+        comp_slc_list,
+        temp_coh_file,
+        ps_looked_file,
+        amp_disp_looked_file,
+    )
 
 
 def create_ifgs(


Test patch:
diff --git a/tests/test_cli.py b/tests/test_cli.py
index 57133424..9ce14dff 100644
--- a/tests/test_cli.py
+++ b/tests/test_cli.py
@@ -31,18 +31,17 @@ def test_subcommand_help(capsys, sub_cmd, option):
 
 
 def test_cli_config_basic(tmpdir, slc_file_list):
-    with tmpdir.as_cwd():
-        with contextlib.suppress(SystemExit):
-            main(
-                [
-                    "config",
-                    "--threads-per-worker",
-                    "2",
-                    "--strides",
-                    "6",
-                    "3",
-                    "--slc-files",
-                    *list(map(str, slc_file_list)),
-                ]
-            )
+    with tmpdir.as_cwd(), contextlib.suppress(SystemExit):
+        main(
+            [
+                "config",
+                "--threads-per-worker",
+                "2",
+                "--strides",
+                "6",
+                "3",
+                "--slc-files",
+                *list(map(str, slc_file_list)),
+            ]
+        )
         assert Path("dolphin_config.yaml").exists()
diff --git a/tests/test_interferogram.py b/tests/test_interferogram.py
index 83299451..1366911b 100644
--- a/tests/test_interferogram.py
+++ b/tests/test_interferogram.py
@@ -136,119 +136,143 @@ def _get_pair_stems(slc_file_pairs):
     ]
 
 
-def test_single_reference_network(tmp_path, four_slc_files):
-    n = Network(four_slc_files, reference_idx=0, outdir=tmp_path)
-
-    assert n.slc_file_pairs[0][0] == four_slc_files[0]
-    assert n.slc_file_pairs[0][1] == four_slc_files[1]
-
-    assert _get_pair_stems(n.slc_file_pairs) == [
-        ("20220101", "20220102"),
-        ("20220101", "20220103"),
-        ("20220101", "20220104"),
-    ]
-    # check the written out files:
-    assert Path(tmp_path / "20220101_20220102.int.vrt").exists()
-    assert Path(tmp_path / "20220101_20220103.int.vrt").exists()
-    assert Path(tmp_path / "20220101_20220104.int.vrt").exists()
-
-    n = Network(four_slc_files, reference_idx=1, outdir=tmp_path)
-    assert _get_pair_stems(n.slc_file_pairs) == [
-        ("20220101", "20220102"),  # still has the same order (early, late)
-        ("20220102", "20220103"),
-        ("20220102", "20220104"),
-    ]
-    assert Path(tmp_path / "20220101_20220102.int.vrt").exists()
-    assert Path(tmp_path / "20220102_20220103.int.vrt").exists()
-    assert Path(tmp_path / "20220102_20220104.int.vrt").exists()
-
-
-def test_limit_by_bandwidth(tmp_path, four_slc_files):
-    n = Network(four_slc_files, max_bandwidth=1, outdir=tmp_path)
-    assert _get_pair_stems(n.slc_file_pairs) == [
-        ("20220101", "20220102"),
-        ("20220102", "20220103"),
-        ("20220103", "20220104"),
-    ]
-    n = Network(four_slc_files, max_bandwidth=2, outdir=tmp_path)
-    assert _get_pair_stems(n.slc_file_pairs) == [
-        ("20220101", "20220102"),
-        ("20220101", "20220103"),
-        ("20220102", "20220103"),
-        ("20220102", "20220104"),
-        ("20220103", "20220104"),
-    ]
-    n = Network(four_slc_files, max_bandwidth=3, outdir=tmp_path)
-    assert _get_pair_stems(n.slc_file_pairs) == [
-        ("20220101", "20220102"),
-        ("20220101", "20220103"),
-        ("20220101", "20220104"),
-        ("20220102", "20220103"),
-        ("20220102", "20220104"),
-        ("20220103", "20220104"),
-    ]
-
-
-def test_limit_by_temporal_baseline(tmp_path, four_slc_files):
-    n = Network(four_slc_files, max_temporal_baseline=0, outdir=tmp_path)
-    assert _get_pair_stems(n.slc_file_pairs) == []
-
-    n = Network(four_slc_files, max_temporal_baseline=1, outdir=tmp_path)
-    assert _get_pair_stems(n.slc_file_pairs) == [
-        ("20220101", "20220102"),
-        ("20220102", "20220103"),
-        ("20220103", "20220104"),
-    ]
-    n = Network(four_slc_files, max_temporal_baseline=2, outdir=tmp_path)
-    assert _get_pair_stems(n.slc_file_pairs) == [
-        ("20220101", "20220102"),
-        ("20220101", "20220103"),
-        ("20220102", "20220103"),
-        ("20220102", "20220104"),
-        ("20220103", "20220104"),
-    ]
-    n = Network(four_slc_files, max_temporal_baseline=500, outdir=tmp_path)
-    assert _get_pair_stems(n.slc_file_pairs) == [
-        ("20220101", "20220102"),
-        ("20220101", "20220103"),
-        ("20220101", "20220104"),
-        ("20220102", "20220103"),
-        ("20220102", "20220104"),
-        ("20220103", "20220104"),
-    ]
-
-
-def test_annual_ifgs():
-    dates = [
-        datetime(2021, 1, 1),
-        datetime(2021, 2, 1),
-        datetime(2022, 1, 4),
-        datetime(2022, 6, 1),
-    ]
-    slcs = [d.strftime("%Y%m%d") for d in dates]
-    assert Network._find_annuals(slcs, dates) == [
-        (slcs[0], slcs[2]),
-        (slcs[1], slcs[2]),
-    ]
-    assert Network._find_annuals(slcs, dates, buffer_days=10) == [(slcs[0], slcs[2])]
-
-    assert Network._find_annuals(slcs[1:], dates[1:], buffer_days=10) == []
-
-
-def test_manual_indexes(tmp_path, four_slc_files):
-    n = Network(four_slc_files, indexes=[], outdir=tmp_path)
-    assert _get_pair_stems(n.slc_file_pairs) == []
-
-    n = Network(four_slc_files, indexes=[(0, 1)], outdir=tmp_path)
-    assert _get_pair_stems(n.slc_file_pairs) == [
-        ("20220101", "20220102"),
-    ]
-
-    n = Network(four_slc_files, indexes=[(0, -1), (0, 1)], outdir=tmp_path)
-    assert _get_pair_stems(n.slc_file_pairs) == [
-        ("20220101", "20220104"),
-        ("20220101", "20220102"),
-    ]
+class TestNetwork:
+    def test_single_reference_network(self, tmp_path, four_slc_files):
+        n = Network(four_slc_files, reference_idx=0, outdir=tmp_path)
+
+        assert n.slc_file_pairs[0][0] == four_slc_files[0]
+        assert n.slc_file_pairs[0][1] == four_slc_files[1]
+
+        assert _get_pair_stems(n.slc_file_pairs) == [
+            ("20220101", "20220102"),
+            ("20220101", "20220103"),
+            ("20220101", "20220104"),
+        ]
+        # check the written out files:
+        assert Path(tmp_path / "20220101_20220102.int.vrt").exists()
+        assert Path(tmp_path / "20220101_20220103.int.vrt").exists()
+        assert Path(tmp_path / "20220101_20220104.int.vrt").exists()
+
+        n = Network(four_slc_files, reference_idx=1, outdir=tmp_path)
+        assert _get_pair_stems(n.slc_file_pairs) == [
+            ("20220101", "20220102"),  # still has the same order (early, late)
+            ("20220102", "20220103"),
+            ("20220102", "20220104"),
+        ]
+        assert Path(tmp_path / "20220101_20220102.int.vrt").exists()
+        assert Path(tmp_path / "20220102_20220103.int.vrt").exists()
+        assert Path(tmp_path / "20220102_20220104.int.vrt").exists()
+
+    def test_limit_by_bandwidth(self, tmp_path, four_slc_files):
+        n = Network(four_slc_files, max_bandwidth=1, outdir=tmp_path)
+        assert _get_pair_stems(n.slc_file_pairs) == [
+            ("20220101", "20220102"),
+            ("20220102", "20220103"),
+            ("20220103", "20220104"),
+        ]
+        n = Network(four_slc_files, max_bandwidth=2, outdir=tmp_path)
+        assert _get_pair_stems(n.slc_file_pairs) == [
+            ("20220101", "20220102"),
+            ("20220101", "20220103"),
+            ("20220102", "20220103"),
+            ("20220102", "20220104"),
+            ("20220103", "20220104"),
+        ]
+        n = Network(four_slc_files, max_bandwidth=3, outdir=tmp_path)
+        assert _get_pair_stems(n.slc_file_pairs) == [
+            ("20220101", "20220102"),
+            ("20220101", "20220103"),
+            ("20220101", "20220104"),
+            ("20220102", "20220103"),
+            ("20220102", "20220104"),
+            ("20220103", "20220104"),
+        ]
+        with pytest.raises(ValueError):
+            Network(four_slc_files, max_bandwidth=0)
+
+    def test_limit_by_temporal_baseline(self, tmp_path, four_slc_files):
+        with pytest.raises(ValueError):
+            Network(four_slc_files, max_temporal_baseline=0)
+
+        n = Network(four_slc_files, max_temporal_baseline=1, outdir=tmp_path)
+        assert _get_pair_stems(n.slc_file_pairs) == [
+            ("20220101", "20220102"),
+            ("20220102", "20220103"),
+            ("20220103", "20220104"),
+        ]
+        n = Network(four_slc_files, max_temporal_baseline=2, outdir=tmp_path)
+        assert _get_pair_stems(n.slc_file_pairs) == [
+            ("20220101", "20220102"),
+            ("20220101", "20220103"),
+            ("20220102", "20220103"),
+            ("20220102", "20220104"),
+            ("20220103", "20220104"),
+        ]
+        n = Network(four_slc_files, max_temporal_baseline=500, outdir=tmp_path)
+        assert _get_pair_stems(n.slc_file_pairs) == [
+            ("20220101", "20220102"),
+            ("20220101", "20220103"),
+            ("20220101", "20220104"),
+            ("20220102", "20220103"),
+            ("20220102", "20220104"),
+            ("20220103", "20220104"),
+        ]
+
+    def test_annual_ifgs(self):
+        dates = [
+            datetime(2021, 1, 1),
+            datetime(2021, 2, 1),
+            datetime(2022, 1, 4),
+            datetime(2022, 6, 1),
+        ]
+        slcs = [d.strftime("%Y%m%d") for d in dates]
+        assert Network._find_annuals(slcs, dates) == [
+            (slcs[0], slcs[2]),
+            (slcs[1], slcs[2]),
+        ]
+        assert Network._find_annuals(slcs, dates, buffer_days=10) == [
+            (slcs[0], slcs[2])
+        ]
+
+        assert Network._find_annuals(slcs[1:], dates[1:], buffer_days=10) == []
+
+    def test_manual_indexes(self, tmp_path, four_slc_files):
+        n = Network(four_slc_files, indexes=[(0, 1)], outdir=tmp_path)
+        assert _get_pair_stems(n.slc_file_pairs) == [
+            ("20220101", "20220102"),
+        ]
+
+        n = Network(four_slc_files, indexes=[(0, -1), (0, 1)], outdir=tmp_path)
+        assert _get_pair_stems(n.slc_file_pairs) == [
+            # It should still always come back sorted
+            ("20220101", "20220102"),
+            ("20220101", "20220104"),
+        ]
+
+    def test_empty_error(self, four_slc_files):
+        with pytest.raises(ValueError):
+            Network(four_slc_files)
+            Network(four_slc_files, indexes=[])
+            Network(four_slc_files, max_bandwidth=0)
+            Network(four_slc_files, max_temporal_baseline=0)
+
+    def test_combination(self, four_slc_files):
+        n1 = Network(four_slc_files, reference_idx=0, write=False)
+        single_ref_pairs = [
+            ("20220101", "20220102"),
+            ("20220101", "20220103"),
+            ("20220101", "20220104"),
+        ]
+        assert _get_pair_stems(n1.slc_file_pairs) == single_ref_pairs
+
+        n2 = Network(
+            four_slc_files, reference_idx=0, indexes=[(1, 2), (1, 3)], write=False
+        )
+        new_ifgs = [
+            ("20220102", "20220103"),
+            ("20220102", "20220104"),
+        ]
+        assert _get_pair_stems(n2.slc_file_pairs) == sorted(single_ref_pairs + new_ifgs)
 
 
 @pytest.fixture()
@@ -332,7 +356,7 @@ def test_network_no_verify():
 
 def test_network_from_ifgs():
     """Check that the `Network` can work when passing in ifgs"""
-    ifg_files = ["20210101_20210107", "202010101_20210108", "20210101_20210109"]
+    ifg_files = ["20210101_20210107", "20210101_20210108", "20210101_20210109"]
     n = Network(
         ifg_files,
         max_bandwidth=10,
diff --git a/tests/test_ps.py b/tests/test_ps.py
index 64b7e717..68c60b17 100644
--- a/tests/test_ps.py
+++ b/tests/test_ps.py
@@ -98,7 +98,10 @@ def test_multilook_ps_file(tmp_path, vrt_stack):
         output_amp_mean_file=amp_mean_file,
         output_file=ps_mask_file,
     )
-    output_file = dolphin.ps.multilook_ps_mask(
-        strides={"x": 5, "y": 3}, ps_mask_file=ps_mask_file
+    output_ps_file, output_amp_disp_file = dolphin.ps.multilook_ps_files(
+        strides={"x": 5, "y": 3},
+        ps_mask_file=ps_mask_file,
+        amp_dispersion_file=amp_dispersion_file,
     )
-    assert io.get_raster_dtype(output_file) == np.uint8
+    assert io.get_raster_dtype(output_ps_file) == np.uint8
+    assert io.get_raster_dtype(output_amp_disp_file) == np.float32
diff --git a/tests/test_timeseries.py b/tests/test_timeseries.py
index c631b1e8..279935c6 100644
--- a/tests/test_timeseries.py
+++ b/tests/test_timeseries.py
@@ -217,7 +217,6 @@ def test_stack(self, data, x_arr, expected_velo):
 
 
 if __name__ == "__main__":
-    # import the fixtures
     sar_dates = make_sar_dates()
     sar_phases = make_sar_phases(sar_dates)
     ifg_date_pairs = make_ifg_date_pairs(sar_dates)
diff --git a/tests/test_workflows_displacement.py b/tests/test_workflows_displacement.py
index ae855ef2..d02e7e51 100644
--- a/tests/test_workflows_displacement.py
+++ b/tests/test_workflows_displacement.py
@@ -25,6 +25,7 @@ def test_displacement_run_single(
             input_options={"subdataset": "/data/VV"},
             interferogram_network={
                 "indexes": [(0, -1)],
+                "max_bandwidth": 2,
             },
             phase_linking={
                 "ministack_size": 500,
@@ -42,6 +43,7 @@ def test_displacement_run_single(
         assert all(p.exists() for p in paths.stitched_cor_paths)
         assert paths.stitched_temp_coh_file.exists()
         assert paths.stitched_ps_file.exists()
+        assert paths.stitched_amp_dispersion_file.exists()
         assert paths.unwrapped_paths is not None
         assert paths.conncomp_paths is not None
         assert all(p.exists() for p in paths.conncomp_paths)