Apply nan/0 mask to output of goldstein filter (#447)

* Apply nan/0 mask to output of `goldstein`

Closes #251, adds citation to paper.

* remove deprecated `complex_`

also see if `nlooks` is causing whirlwind to freeze

* fix numpy 2.0 deprecations by ruff

* fix sizes

* test spurt branch

.github/workflows/test-build-push.yml

@@ -51,7 +51,7 @@ jobs:
             "opera-utils>=0.4.1" \
             git+https://github.com/isce-framework/tophu@main \
             git+https://github.com/isce-framework/whirlwind@40defb38d2d6deca2819934788ebbc57e418e32d
-          python -m pip install git+https://github.com/isce-framework/spurt@main
+          python -m pip install git+https://github.com/scottstanie/spurt@use-mp-spawn
           python -m pip install --no-deps .
       - name: Install test dependencies
         run: |

docs/references.bib

@@ -165,6 +165,18 @@ @article{Fornaro2015CAESARApproachBased
   keywords = {3-D,4-D and multidimensional (Multi-D) SAR imaging,Covariance matrices,Covariance matrix decomposition,differential SAR tomography,differential synthetic aperture radar (SAR) interferometry (DInSAR),Interferometry,Monitoring,principal component analysis (PCA),SAR interferometry (InSAR),SAR tomography,Scattering,Spatial resolution,Synthetic aperture radar,Tomography}
 }
 
+@article{Goldstein1998RadarInterferogramFiltering,
+  title = {Radar Interferogram Filtering for Geophysical Applications},
+  author = {Goldstein, Richard M. and Werner, Charles L.},
+  year = {1998},
+  journal = {Geophysical Research Letters},
+  volume = {25},
+  number = {21},
+  pages = {4035--4038},
+  issn = {1944-8007},
+  doi = {10.1029/1998GL900033}
+}
+
 @article{Guarnieri2008ExploitationTargetStatistics,
   title = {On the {{Exploitation}} of {{Target Statistics}} for {{SAR Interferometry Applications}}},
   author = {Guarnieri, A. M. and Tebaldini, S.},

pyproject.toml

@@ -64,6 +64,7 @@ select = [
   "I",   # isort
   "ISC", # flake8-implicit-str-concat
   "N",   # pep8-naming
+  "NPY201", # numpy 2.0 depgrcations
   "PGH", # pygrep-hooks
   "PIE", # flake8-pie
   "PL",  # Pylint

src/dolphin/atmosphere/weather_model.py

@@ -342,10 +342,10 @@ def _get_ztd(self):
         # Get the integrated ZTD
         wet_total, hydro_total = np.zeros(wet.shape), np.zeros(hydro.shape)
         for level in range(wet.shape[2]):
-            wet_total[..., level] = 1e-6 * np.trapz(
+            wet_total[..., level] = 1e-6 * np.trapezoid(
                 wet[..., level:], x=self._zs[level:], axis=2
             )
-            hydro_total[..., level] = 1e-6 * np.trapz(
+            hydro_total[..., level] = 1e-6 * np.trapezoid(
                 hydro[..., level:], x=self._zs[level:], axis=2
             )
         self._hydrostatic_ztd = hydro_total

src/dolphin/goldstein.py

@@ -1,14 +1,16 @@
 import numpy as np
-from numpy.typing import ArrayLike
+from numpy.typing import NDArray
 
 
-def goldstein(phase: ArrayLike, alpha: float, psize: int = 32) -> np.ndarray:
+def goldstein(
+    phase: NDArray[np.complex64] | NDArray[np.float64], alpha: float, psize: int = 32
+) -> np.ndarray:
     """Apply the Goldstein adaptive filter to the given data.
 
     Parameters
     ----------
     phase : np.ndarray
-        2D complex array containing the data to be filtered.
+        2D array of floating point phase, or complex data, to be filtered.
     alpha : float
         Filtering parameter for Goldstein algorithm
         Must be between 0 (no filtering) and 1 (maximum filtering)
@@ -20,43 +22,47 @@ def goldstein(phase: ArrayLike, alpha: float, psize: int = 32) -> np.ndarray:
     -------
         2D numpy array of filtered data.
 
+    References
+    ----------
+    [@Goldstein1998RadarInterferogramFiltering]
+
     """
 
-    def apply_pspec(data: ArrayLike):
+    def apply_pspec(data: NDArray[np.complex64]) -> np.ndarray:
         # NaN is allowed value
         if alpha < 0:
             raise ValueError(f"alpha must be >= 0, got {alpha = }")
 
-        wgt = np.power(np.abs(data) ** 2, alpha / 2)
-        data = wgt * data
+        weight = np.power(np.abs(data) ** 2, alpha / 2)
+        data = weight * data
         return data
 
-    def make_wgt(nxp: int, nyp: int) -> np.ndarray:
+    def make_weight(nxp: int, nyp: int) -> np.ndarray:
         # Create arrays of horizontal and vertical weights
         wx = 1.0 - np.abs(np.arange(nxp // 2) - (nxp / 2.0 - 1.0)) / (nxp / 2.0 - 1.0)
         wy = 1.0 - np.abs(np.arange(nyp // 2) - (nyp / 2.0 - 1.0)) / (nyp / 2.0 - 1.0)
         # Compute the outer product of wx and wy to create
         # the top-left quadrant of the weight matrix
         quadrant = np.outer(wy, wx)
         # Create a full weight matrix by mirroring the quadrant along both axes
-        wgt = np.block(
+        weight = np.block(
             [
                 [quadrant, np.flip(quadrant, axis=1)],
                 [np.flip(quadrant, axis=0), np.flip(np.flip(quadrant, axis=0), axis=1)],
             ]
         )
-        return wgt
+        return weight
 
     def patch_goldstein_filter(
-        data: ArrayLike, wgt: ArrayLike, psize: int
+        data: NDArray[np.complex64], weight: NDArray[np.float64], psize: int
     ) -> np.ndarray:
         """Apply the filter to a single patch of data.
 
         Parameters
         ----------
         data : np.ndarray
             2D complex array containing the data to be filtered.
-        wgt : np.ndarray
+        weight : np.ndarray
             weight matrix for summing neighboring data
         psize : int
             edge length of square FFT area
@@ -70,30 +76,39 @@ def patch_goldstein_filter(
         data = np.fft.fft2(data, s=(psize, psize))
         data = apply_pspec(data)
         data = np.fft.ifft2(data, s=(psize, psize))
-        return wgt * data
+        return weight * data
 
-    def apply_goldstein_filter(data: ArrayLike) -> np.ndarray:
+    def apply_goldstein_filter(data: NDArray[np.complex64]) -> np.ndarray:
         # Create an empty array for the output
         out = np.zeros(data.shape, dtype=np.complex64)
-        # ignore processing for empty chunks
-        if np.all(np.isnan(data)):
+        empty_mask = np.isnan(data) | (np.angle(data) == 0)
+        # ignore processing for a chunks
+        if np.all(empty_mask):
             return data
         # Create the weight matrix
-        wgt_matrix = make_wgt(psize, psize)
+        weight_matrix = make_weight(psize, psize)
         # Iterate over windows of the data
         for i in range(0, data.shape[0] - psize, psize // 2):
             for j in range(0, data.shape[1] - psize, psize // 2):
-                # Create proocessing windows
+                # Create processing windows
                 data_window = data[i : i + psize, j : j + psize]
-                wgt_window = wgt_matrix[: data_window.shape[0], : data_window.shape[1]]
+                weight_window = weight_matrix[
+                    : data_window.shape[0], : data_window.shape[1]
+                ]
                 # Apply the filter to the window
-                filtered_window = patch_goldstein_filter(data_window, wgt_window, psize)
+                filtered_window = patch_goldstein_filter(
+                    data_window, weight_window, psize
+                )
                 # Add the result to the output array
                 slice_i = slice(i, min(i + psize, out.shape[0]))
                 slice_j = slice(j, min(j + psize, out.shape[1]))
                 out[slice_i, slice_j] += filtered_window[
                     : slice_i.stop - slice_i.start, : slice_j.stop - slice_j.start
                 ]
+        out[empty_mask] = 0
         return out
 
-    return apply_goldstein_filter(phase)
+    if np.iscomplexobj(phase):
+        return apply_goldstein_filter(phase)
+    else:
+        return apply_goldstein_filter(np.exp(1j * phase))

tests/make_netcdf.py

@@ -15,7 +15,7 @@ def _add_complex_type(h5_root_group):
     if "complex64" in h5_root_group:
         return
     ctype = h5py.h5t.py_create(np.complex64)
-    ctype.commit(h5_root_group.id, np.string_("complex64"))
+    ctype.commit(h5_root_group.id, np.bytes_("complex64"))
 
 
 def create_test_nc(

tests/test_unwrap.py

@@ -322,18 +322,18 @@ def test_unwrap_multiscale_callback_given(
 
 class TestSpurt:
     @pytest.fixture()
-    def ifg_file_list(self, tmp_path, slc_stack, slc_date_list):
+    def ifg_file_list(self, tmp_path, slc_date_list):
         from dolphin import io
         from dolphin.phase_link import simulate
 
-        slc_stack = np.exp(
-            1j * simulate.make_defo_stack((20, 100, 200), sigma=1)
-        ).astype("complex64")
+        slc_stack = np.exp(1j * simulate.make_defo_stack((10, 20, 30), sigma=1)).astype(
+            "complex64"
+        )
         ifg_stack = slc_stack[1:] * slc_stack[[0]].conj()
         # Write to a file
         d = tmp_path / "gtiff"
         d.mkdir()
-        name_template = d / f"{slc_date_list[0].strftime('%Y%m%d')}_{{date}}.slc.tif"
+        name_template = d / f"{slc_date_list[0].strftime('%Y%m%d')}_{{date}}.int.tif"
 
         file_list = []
         for cur_date, cur_ifg in zip(slc_date_list[1:], ifg_stack):
@@ -343,16 +343,27 @@ def ifg_file_list(self, tmp_path, slc_stack, slc_date_list):
 
         return file_list
 
+    @pytest.fixture()
+    def temp_coh_raster(self, ifg_file_list):
+        d = Path(ifg_file_list[0]).parent
+        coh_raster = d / "temporal_coherence.tif"
+        io.write_arr(
+            arr=np.random.rand(20, 30).astype(np.float32),
+            output_name=coh_raster,
+            driver="GTiff",
+        )
+        return coh_raster
+
     @pytest.mark.skipif(
         not SPURT_INSTALLED, reason="spurt not installed for 3d unwrapping"
     )
-    def test_unwrap_spurt(self, tmp_path, ifg_file_list, corr_raster):
+    def test_unwrap_spurt(self, tmp_path, ifg_file_list, temp_coh_raster):
         opts = SpurtOptions()
         unwrap_options = UnwrapOptions(unwrap_method="spurt", spurt_options=opts)
         out_paths, conncomp_paths = dolphin.unwrap.run(
             ifg_filenames=ifg_file_list,
             cor_filenames=ifg_file_list,  # NOT USED... but required for `run`?
-            temporal_coherence_file=corr_raster,
+            temporal_coherence_file=temp_coh_raster,
             unwrap_options=unwrap_options,
             output_path=tmp_path,
             nlooks=5,