Skip nodata values of interferograms when interpolating (#476)

* Skip nodata values of interferograms when interpolating

Currently the outside, nodata values will all fall below the threshold, leading to the most possible pixels summed for the area we dont care about.

* Move interpolation test to separate test module

src/dolphin/interpolation.py

@@ -64,6 +64,7 @@ def interpolate(
 
     """
     nrow, ncol = weights.shape
+    ifg_is_valid_mask = ifg != 0
 
     weights_float = np.clip(weights.astype(np.float32), 0, 1)
     # Ensure weights are between 0 and 1
@@ -83,6 +84,7 @@ def interpolate(
         ifg,
         weights_float,
         weight_cutoff,
+        ifg_is_valid_mask,
         num_neighbors,
         alpha,
         indices,
@@ -93,12 +95,21 @@ def interpolate(
 
 @numba.njit(parallel=True)
 def _interp_loop(
-    ifg, weights, weight_cutoff, num_neighbors, alpha, indices, interpolated_ifg
+    ifg,
+    weights,
+    weight_cutoff,
+    ifg_is_valid_mask,
+    num_neighbors,
+    alpha,
+    indices,
+    interpolated_ifg,
 ):
     nrow, ncol = weights.shape
     nindices = len(indices)
     for r0 in numba.prange(nrow):
         for c0 in range(ncol):
+            if not ifg_is_valid_mask[r0, c0]:
+                continue
             if weights[r0, c0] >= weight_cutoff:
                 interpolated_ifg[r0, c0] = ifg[r0, c0]
                 continue

src/dolphin/unwrap/_tophu.py

@@ -186,7 +186,7 @@ def _get_cb_and_nodata(unwrap_method, unwrap_callback, nodata):
         logger.info(f"Zeroing unw/conncomp of pixels masked in {mask_file}")
         return _zero_from_mask(unw_filename, conncomp_filename, mask_file)
     elif unwrap_method == UnwrapMethod.PHASS:
-        # Fill in the nan pixels with the nearest amgibuities
+        # Fill in the nan pixels with the nearest ambiguities
         from ._post_process import interpolate_masked_gaps
 
         with rio.open(unw_filename, mode="r+") as u_src, rio.open(

src/dolphin/unwrap/_unwrap.py

@@ -367,7 +367,8 @@ def unwrap(
             str(pre_interp_unw_filename).split(".")[0] + (name_change + "unw" + suf)
         )
 
-        pre_interp_ifg = io.load_gdal(pre_interp_ifg_filename)
+        pre_interp_ifg = io.load_gdal(pre_interp_ifg_filename, masked=True).filled(0)
+
         corr = io.load_gdal(corr_filename)
         if similarity_filename and preproc_options.interpolation_similarity_threshold:
             cutoff = preproc_options.interpolation_similarity_threshold

tests/test_interpolation.py

@@ -0,0 +1,34 @@
+import numpy as np
+
+from dolphin import interpolation
+
+
+def test_interpolate():
+    x, y = np.meshgrid(np.arange(200), np.arange(100))
+    # simulate a simple phase ramp
+    phase = 0.003 * x + 0.002 * y
+    # interferogram with the simulated phase ramp and with a constant amplitude
+    ifg = np.exp(1j * phase)
+    corr = np.ones(ifg.shape)
+    # mask out the ifg/corr at a given pixel for example at pixel 50,40
+    # index of the pixel of interest
+    x_idx = 50
+    y_idx = 40
+    corr[y_idx, x_idx] = 0
+
+    interpolated_ifg = np.zeros((100, 200), dtype=np.complex64)
+    interpolation.interpolate(
+        ifg,
+        weights=corr,
+        weight_cutoff=0.5,
+        num_neighbors=20,
+        alpha=0.75,
+        max_radius=51,
+        min_radius=0,
+    )
+    # expected phase based on the model above used for simulation
+    expected_phase = 0.003 * x_idx + 0.002 * y_idx
+    phase_error = np.angle(
+        interpolated_ifg[y_idx, x_idx] * np.exp(-1j * expected_phase)
+    )
+    assert np.allclose(phase_error, 0.0, atol=1e-3)

tests/test_unwrap.py

@@ -192,42 +192,6 @@ def test_multiple_preprocess(
         unw_dir = Path(unw_path).parent
         assert set(unw_dir.glob("*.unw.tif")) == {unw_path}
 
-    def test_interp_loop(self):
-        x, y = np.meshgrid(np.arange(200), np.arange(100))
-        # simulate a simple phase ramp
-        phase = 0.003 * x + 0.002 * y
-        # interferogram with the simulated phase ramp and with a constant amplitude
-        ifg = np.exp(1j * phase)
-        corr = np.ones(ifg.shape)
-        # mask out the ifg/corr at a given pixel for example at pixel 50,40
-        # index of the pixel of interest
-        x_idx = 50
-        y_idx = 40
-        corr[y_idx, x_idx] = 0
-        # generate indices for the pixels to be used in interpolation
-        indices = np.array(
-            dolphin.similarity.get_circle_idxs(
-                max_radius=51, min_radius=0, sort_output=False
-            )
-        )
-        # interpolate pixels with zero in the corr and write to interpolated_ifg
-        interpolated_ifg = np.zeros((100, 200), dtype=np.complex64)
-        dolphin.interpolation._interp_loop(
-            ifg,
-            corr,
-            weight_cutoff=0.5,
-            num_neighbors=20,
-            alpha=0.75,
-            indices=indices,
-            interpolated_ifg=interpolated_ifg,
-        )
-        # expected phase based on the model above used for simulation
-        expected_phase = 0.003 * x_idx + 0.002 * y_idx
-        phase_error = np.angle(
-            interpolated_ifg[y_idx, x_idx] * np.exp(-1j * expected_phase)
-        )
-        assert np.allclose(phase_error, 0.0, atol=1e-3)
-
     @pytest.mark.parametrize("method", [UnwrapMethod.SNAPHU, UnwrapMethod.PHASS])
     def test_interpolation(self, tmp_path, list_of_gtiff_ifgs, corr_raster, method):
         # test other init_method