Goldstein filtering step for unwrapping workflow

[rtburns-jpl]

No description provided.

[scottstanie]

1. i'll add an issue about the nan block problem to fix in a later PR
2. let's add a simple test like this (from my git diff

diff --git a/tests/test_unwrap.py b/tests/test_unwrap.py
index 8789820..10f16c4 100644
--- a/tests/test_unwrap.py
+++ b/tests/test_unwrap.py
@@ -97,6 +97,21 @@ class TestUnwrapSingle:
         assert u_path.exists()
         assert c_path.exists()

+    @pytest.mark.parametrize("method", [UnwrapMethod.SNAPHU, UnwrapMethod.PHASS])
+    def test_goldstein(self, tmp_path, list_of_gtiff_ifgs, corr_raster, method):
+        # test other init_method
+        unw_filename = tmp_path / "unwrapped.unw.tif"
+        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_method=method,
+            run_goldstein=True,
+        )
+        assert unw_path.exists()
+        assert conncomp_path.exists()
+

[scottstanie]

we might end up changing this name slightly to avoid confusion with the SHP alpha parameter but for now it's fine and different

[scottstanie]

Was it causing any problems to have this imported at the top of the file? Or was this a mistake to move the import?

[rtburns-jpl]

Oops, yeah I moved this here because I was testing in an environment without snaphu. It was very convenient to not have to install snaphu simply by moving the import to where it's needed, but I can definitely move this back if that's an antipattern.

[scottstanie]

I suppose we can leave it for now if someone really want's to use another unwrapper that's not snaphu py, we're not currently putting snaphu-py in the environment file

[scottstanie]

Note for @mirzaees and @taliboliver - this is where it could make sense to add the interpolation logic.
If this section gets unwieldy, perhaps we could break it into a helper function that returns the unwrapper ifg_filename and unw_filename

[gmgunter]

Sorry for the late review. I started taking a look yesterday but was too slow to finish editing my comments. I'll post the comments for posterity but overall I think the code looks good and mostly just want to complain about the lack of unit testing.

I posted some suggestions for unit tests below without providing much detail. If you ever want to talk through the details of how we might implement these, I'd be happy to tag up.

[gmgunter]

Naming the input argument "phase" might be a little misleading. The Goldstein filter is intended to operate on complex-valued interferogram data -- not on the wrapped phase directly.

I typically name this parameter "igram".

Suggested change

	def goldstein(phase: ArrayLike, alpha: float, psize: int = 32) -> np.ndarray:
	def goldstein(igram: ArrayLike, alpha: float, psize: int = 32) -> np.ndarray:

[scottstanie]

ah I lost this same comment in a closed window- let's definitely change it for the next PR because it confused me several times

[gmgunter]

Let's add a reference to the Goldstein/Werner 1998 paper here.

@scottstanie has made it really easy to cite references in docstrings in this repo--

You could add a bibtex entry like this

@article{Goldstein1998RadarInterferogramFiltering,
    title = {Radar interferogram filtering for geophysical applications},
    author = {Goldstein, Richard M. and Werner, Charles L.},
    year = {1998},
    month = nov,
    journal = {Geophysical Research Letters},
    volume = {25},
    number = {21},
    pages = {4035–4038}
}

to this file: docs/references.bib

Then you could cite the paper in a docstring like this:

Suggested change

	"""Apply the Goldstein adaptive filter to the given data.
	"""Apply the Goldstein adaptive filter [@Goldstein1998RadarInterferogramFiltering] to the given data.

[gmgunter]

The annotation above says this argument is of type ArrayLike. Usually types like this are documented like this in the NumPy docs:

Suggested change

	phase : np.ndarray
	phase : array_like

where "array_like" is a keyword in the NumPy glossary: https://numpy.org/doc/stable/glossary.html#term-array_like

[gmgunter]

I prefer to use full sentences in docstring descriptions (with periods and proper capitalization), since I think it looks more professional.

Suggested change

	Filtering parameter for Goldstein algorithm
	Must be between 0 (no filtering) and 1 (maximum filtering)
	psize : int, optional
	edge length of square patch
	Default = 32
	Filtering parameter for the Goldstein algorithm.
	Must be between 0 (no filtering) and 1 (maximum filtering).
	psize : int, optional
	Edge length of square patch. Defaults to 32.

[gmgunter]

Must be between 0 (no filtering) and 1 (maximum filtering)

Should we enforce this in the function body? E.g. something like this:

if (alpha < 0.0) or (alpha > 1.0):
    raise ValueError(f"alpha must be between 0 and 1, instead got {alpha}")

[gmgunter]

Similarly, you could consider adding a check that psize has a valid value.

I think psize should be at least >=2, since if psize is 1, then the block stride (psize // 2) will be 0.

[gmgunter]

Do we really need any special handling for this case?

Let's say hypothetically we remove lines 77 & 78. In that case, if the input data contained all NaN values, I expect that the output filtered data will still contain all NaNs. So lines 77 & 78 are just optimizing that case -- the results should still be the same but we get to skip some unnecessary computation.

But should we really be optimizing for the case where the input data is all NaN-valued? Note that this is (slightly) pessimizing the usual case where the input contains at least one non-NaN value, since np.all(np.isnan(data)) takes relatively small but nonzero time.

I would suggest removing these lines since I think encountering all-NaN-valued datasets seems unlikely in normal operation so, on net, this is probably a pessimization.

Suggested change

	# ignore processing for empty chunks
	if np.all(np.isnan(data)):
	return data

[scottstanie]

This could be a good thing to test- For the stitched, geocoded interferograms, The tilt means it's about 30% nodata. For single bursts, it's about 60% nodata (which is why adding a check like this speeds up the wrapped phase part by ~2x).

On the other hand, having this hardcoded is very geocoded-input-centric. An RIFG has no need for this and only gets it slowed down

[rtburns-jpl]

Yeah will remove this. I think this is a remnant from when I was nanning out zero values, which can't be done with tophu at the moment anyway.

[gmgunter]

For the stitched, geocoded interferograms, The tilt means it's about 30% nodata. For single bursts, it's about 60% nodata (which is why adding a check like this speeds up the wrapped phase part by ~2x).

Right now, the check for all NaN values is outside the loop over patches. So this shortcut will only trigger if the dataset contains 100% nodata. IIUC, there shouldn't be any speed up for interferograms that contain 60% nodata.

Maybe we should consider adding a similar check inside the loop over patches? If an individual patch contains all NaNs, there's no need to compute & apply the filter there.

[scottstanie]

oh right, yes this should be removed then, the comment # ignore processing for empty chunks fooled me

[gmgunter]

Noting another potential issue with odd-valued psize here--

If psize is even-valued, then each pixel will be covered by exactly 4 overlapping blocks (except for pixels near the borders of the data array.

But if psize is odd-valued, then the last row & column of pixels in each block will each be covered by an additional block.

[gmgunter]

Do we care about the amplitude values of the filtered data at all? Maybe a question for @scottstanie.

Currently, the filter coefficients and block weights are not normalized, so the processing is effectively applying a scaling factor to the output. It shouldn't affect the phase, but the amplitude values will not be comparable to the input amplitudes.

[scottstanie]

I don't believe so- snaphu just ignores the amplitude if you also pass the correlation, right?

[gmgunter]

Yeah, I believe you're right about SNAPHU. The ICU algorithm does care of about interferogram amplitudes a bit (for computing "intensity neutrons"), but I don't think we're seriously considering ICU in our comparison of unwrappers. I'm not aware of other unwrapping algorithms that use amplitude information, so it seems reasonably safe to not worry about normalizing amplitudes.

[gmgunter]

No unit tests ☹️

If you make this a standalone function (so that its importable from the test suite), it should be pretty trivial to implement a test that checks that the phase values are unchanged by filtering with alpha equal to 0. I think that would be a nice test.

Another good test would be to simulate a signal with a known magnitude spectrum (e.g. by constructing the signal in the frequency domain and then taking the IFFT) and then applying patch_goldstein_filter() to it. Then you could compute the magnitude spectrum of the filtered data and compare it to the expected spectrum (which should be the original spectrum to the power 1+alpha).

It might also be worthwhile to exercise the blocking logic in a test. Common issues in block processing include off-by-one indexing errors and mishandling the remainder samples past the end of the last full-sized block. I think you could guard against most of these common issues in a test by passing an array of ones to the goldstein() function and checking that the output array values are constant (i.e. are all the same) to within some tolerance.

[scottstanie]

related to the future nan change #261 , we could check that the number of input nodata values are similar/the same as the output nodata values.

[gmgunter]

Looks like this is baking in the assumption that we will only ever want to process GeoTiff files or ENVI files.

It'd be nice if dolphin had a function like

def get_raster_driver(filename: str | os.PathLike[str]) -> str:
    ...

that abstracted the logic for inferring the driver name of a particular raster, so that any future changes to this logic could be contained in a single function.

Similarly, it'd be nice to see a function like

def get_default_driver_options(driver: str) -> list[str]:
    ...

as the single source of truth for getting our default driver-specific raster creation options.

[scottstanie]

I also meant to comment here- get_raster_driver is a good one to add to the other similar io functions, but here I might actually move to cut the ENVI one and just produce geotiff outputs at this stage. I think Ryan got this ENVI check back when I was half way out of only setting up ENVI files, way before the wrapper/tophu was set up.

If we do want to accomodate more than GTiff, that would probably make more sense as a larger refactor to pass around DatasetWriters instead of Filenames. It would just require more thought how you would swap the user's DatasetWriter for a temporary one to hold the goldstein output