Set the final_arr nodata to be same as original ifg mask (#356)

* Set the `final_arr` nodata to be same as original ifg mask

closes #353

* move unwrapped, not ifg

* add more tests for combinations of goldstein/interp

src/dolphin/unwrap/_unwrap.py

@@ -270,6 +270,7 @@ def unwrap(
     unwrapper_unw_filename = Path(unw_filename)
     name_change = "."
 
+    ifg = io.load_gdal(ifg_filename, masked=True)
     if unwrap_options.run_goldstein:
         suf = Path(unw_filename).suffix
         if suf == ".tif":
@@ -289,9 +290,8 @@ def unwrap(
             str(unw_filename).split(".")[0] + (name_change + "unw" + suf)
         )
 
-        ifg = io.load_gdal(ifg_filename)
         logger.info(f"Goldstein filtering {ifg_filename} -> {filt_ifg_filename}")
-        modified_ifg = goldstein(ifg, alpha=preproc_options.alpha)
+        modified_ifg = goldstein(ifg.filled(0), alpha=preproc_options.alpha)
         logger.info(f"Writing filtered output to {filt_ifg_filename}")
         io.write_arr(
             arr=modified_ifg,
@@ -324,15 +324,15 @@ def unwrap(
             str(pre_interp_unw_filename).split(".")[0] + (name_change + "unw" + suf)
         )
 
-        ifg = io.load_gdal(pre_interp_ifg_filename)
+        pre_interp_ifg = io.load_gdal(pre_interp_ifg_filename)
         corr = io.load_gdal(corr_filename)
         cutoff = preproc_options.interpolation_cor_threshold
         logger.info(f"Masking pixels with correlation below {cutoff}")
         coherent_pixel_mask = corr[:] >= cutoff
 
         logger.info(f"Interpolating {pre_interp_ifg_filename} -> {interp_ifg_filename}")
         modified_ifg = interpolate(
-            ifg=ifg,
+            ifg=pre_interp_ifg,
             weights=coherent_pixel_mask,
             weight_cutoff=cutoff,
             max_radius=preproc_options.max_radius,
@@ -405,26 +405,17 @@ def unwrap(
 
     # post-processing steps go here:
 
-    # Reset the input nodata values to be nodata in the unwrapped and CCL
-    logger.info(f"Setting nodata values of {unw_path} file")
-    set_nodata_values(
-        filename=unw_path, output_nodata=unw_nodata, like_filename=ifg_filename
-    )
-    logger.info(f"Setting nodata values of {conncomp_path} file")
-    set_nodata_values(
-        filename=conncomp_path, output_nodata=ccl_nodata, like_filename=ifg_filename
-    )
-
     # Transfer ambiguity numbers from filtered/interpolated unwrapped interferogram
     # back to original interferogram
     if unwrap_options.run_goldstein or unwrap_options.run_interpolation:
         logger.info(
             "Transferring ambiguity numbers from filtered/interpolated"
             f" ifg {unwrapper_unw_filename}"
         )
-        unw_arr = io.load_gdal(unwrapper_unw_filename)
+        unw_arr = io.load_gdal(unwrapper_unw_filename, masked=True).filled(unw_nodata)
 
         final_arr = np.angle(ifg) + (unw_arr - np.angle(modified_ifg))
+        final_arr[ifg.mask] = unw_nodata
 
         io.write_arr(
             arr=final_arr,
@@ -450,6 +441,20 @@ def unwrap(
             scratchdir=scratchdir,
         )
 
+        # Move the intermediate ".interp" or ".goldstein" into the scratch directory
+        if scratchdir is not None:
+            shutil.move(unwrapper_unw_filename, scratchdir)
+
+    # Reset the input nodata values to be nodata in the unwrapped and CCL
+    logger.info(f"Setting nodata values of {unw_path} file")
+    set_nodata_values(
+        filename=unw_filename, output_nodata=unw_nodata, like_filename=ifg_filename
+    )
+    logger.info(f"Setting nodata values of {conncomp_path} file")
+    set_nodata_values(
+        filename=conncomp_path, output_nodata=ccl_nodata, like_filename=ifg_filename
+    )
+
     if delete_scratch:
         assert scratchdir is not None
         shutil.rmtree(scratchdir)

tests/test_unwrap.py

@@ -131,15 +131,18 @@ def test_unwrap_snaphu_nodata(
         assert np.isnan(io.get_raster_nodata(unw_path))
         assert io.get_raster_nodata(conncomp_path) == 123
 
-    @pytest.mark.parametrize("method", [UnwrapMethod.SNAPHU, UnwrapMethod.PHASS])
     def test_goldstein(
-        self, tmp_path, list_of_gtiff_ifgs, corr_raster, unwrap_options, method
+        self,
+        tmp_path,
+        list_of_gtiff_ifgs,
+        corr_raster,
+        unwrap_options,
     ):
         # test other init_method
         unw_filename = tmp_path / "unwrapped.unw.tif"
         scratch_dir = tmp_path / "scratch"
         scratch_dir.mkdir()
-        unwrap_options.unwrap_method = method
+
         unwrap_options.run_goldstein = True
         unw_path, conncomp_path = dolphin.unwrap.unwrap(
             ifg_filename=list_of_gtiff_ifgs[0],
@@ -151,6 +154,43 @@ def test_goldstein(
         )
         assert unw_path.exists()
         assert conncomp_path.exists()
+        # Check there are no extraneous ".filter" files
+        unw_dir = Path(unw_path).parent
+        assert set(unw_dir.glob("*.unw.tif")) == {unw_path}
+
+    @pytest.mark.parametrize("method", [UnwrapMethod.SNAPHU, UnwrapMethod.PHASS])
+    @pytest.mark.parametrize("run_interpolation", [False, True])
+    @pytest.mark.parametrize("run_goldstein", [False, True])
+    def test_multiple_preprocess(
+        self,
+        tmp_path,
+        list_of_gtiff_ifgs,
+        corr_raster,
+        unwrap_options,
+        method,
+        run_interpolation,
+        run_goldstein,
+    ):
+        # test other init_method
+        unw_filename = tmp_path / "unwrapped.unw.tif"
+        scratch_dir = tmp_path / "scratch"
+        scratch_dir.mkdir()
+        unwrap_options.unwrap_method = method
+        unwrap_options.run_interpolation = run_interpolation
+        unwrap_options.run_goldstein = run_goldstein
+        unw_path, conncomp_path = dolphin.unwrap.unwrap(
+            ifg_filename=list_of_gtiff_ifgs[0],
+            corr_filename=corr_raster,
+            unw_filename=unw_filename,
+            nlooks=1,
+            unwrap_options=unwrap_options,
+            scratchdir=scratch_dir,
+        )
+        assert unw_path.exists()
+        assert conncomp_path.exists()
+        # Check there are no extraneous ".interp" files
+        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))
@@ -205,6 +245,9 @@ def test_interpolation(self, tmp_path, list_of_gtiff_ifgs, corr_raster, method):
         )
         assert unw_path.exists()
         assert conncomp_path.exists()
+        # Check there are no extraneous ".interp" files
+        unw_dir = Path(unw_path).parent
+        assert set(unw_dir.glob("*.unw.tif")) == {unw_path}
 
 
 class TestUnwrapRun:
@@ -222,7 +265,7 @@ def test_run_gtiff(self, list_of_gtiff_ifgs, corr_raster, unwrap_options):
 
     def test_run_envi(self, list_of_envi_ifgs, corr_raster, unwrap_options):
         ifg_path = list_of_envi_ifgs[0].parent
-        unwrap_options.snaphu_options.init_method = "mst"
+        unwrap_options.snaphu_options.init_method = "mcf"
         u_paths, c_paths = dolphin.unwrap.run(
             ifg_filenames=list_of_envi_ifgs,
             cor_filenames=[corr_raster] * len(list_of_envi_ifgs),