filtering

filter_long_wavelength(unwrapped_phase, bad_pixel_mask, wavelength_cutoff=25 * 1000.0, pixel_spacing=30, workers=1, fill_value=None, scratch_dir=None)

Filter out signals with spatial wavelength longer than a threshold.

Parameters:

Name	Type	Description	Default
unwrapped_phase	ArrayLike	Unwrapped interferogram phase to filter.	required
bad_pixel_mask	ArrayLike	Boolean array with same shape as unwrapped_phase where True indicates a pixel to ignore during missing-data filling.	required
wavelength_cutoff	float	Spatial wavelength threshold to filter the unwrapped phase. Signals with wavelength longer than 'wavelength_cutoff' are filtered out. The default is 25*1e3 (meters).	25 * 1000.0
pixel_spacing	float	Pixel spatial spacing. Assume same spacing for x, y axes. The default is 30 (meters).	30
workers	int	Number of fft workers to use for scipy.fft.fft2. Default is 1.	1
fill_value	float	Value to place in output pixels which were masked. If None, masked pixels are filled with interpolated values using gdal_fillnodata before filtering to suppress outliers.	None
scratch_dir	Path	Directory to use for temporary files. If not provided, uses system default for Python's tempfile module.	None

Returns:

Name	Type	Description
filtered_ifg	2D complex array	filtered interferogram that does not contain signals with spatial wavelength longer than a threshold.

Raises:

Type	Description
ValueError	If wavelength_cutoff too large for image size/pixel spacing.

Source code in src/dolphin/filtering.py

def filter_long_wavelength(
    unwrapped_phase: ArrayLike,
    bad_pixel_mask: ArrayLike,
    wavelength_cutoff: float = 25 * 1e3,
    pixel_spacing: float = 30,
    workers: int = 1,
    fill_value: float | None = None,
    scratch_dir: Path | None = None,
) -> np.ndarray:
    """Filter out signals with spatial wavelength longer than a threshold.

    Parameters
    ----------
    unwrapped_phase : ArrayLike
        Unwrapped interferogram phase to filter.
    bad_pixel_mask : ArrayLike
        Boolean array with same shape as `unwrapped_phase` where `True` indicates a
        pixel to ignore during missing-data filling.
    wavelength_cutoff : float
        Spatial wavelength threshold to filter the unwrapped phase.
        Signals with wavelength longer than 'wavelength_cutoff' are filtered out.
        The default is 25*1e3 (meters).
    pixel_spacing : float
        Pixel spatial spacing. Assume same spacing for x, y axes.
        The default is 30 (meters).
    workers : int
        Number of `fft` workers to use for `scipy.fft.fft2`.
        Default is 1.
    fill_value : float, optional
        Value to place in output pixels which were masked.
        If `None`, masked pixels are filled with interpolated values
        using `gdal_fillnodata` before filtering to suppress outliers.
    scratch_dir : Path, optional
        Directory to use for temporary files. If not provided, uses system default
        for Python's tempfile module.

    Returns
    -------
    filtered_ifg : 2D complex array
        filtered interferogram that does not contain signals with spatial wavelength
        longer than a threshold.

    Raises
    ------
    ValueError
        If wavelength_cutoff too large for image size/pixel spacing.

    """
    from osgeo_utils.gdal_fillnodata import gdal_fillnodata
    from scipy import fft, ndimage

    from dolphin import io

    # Find the filter `sigma` which gives the correct cutoff in meters
    sigma = _compute_filter_sigma(wavelength_cutoff, pixel_spacing, cutoff_value=0.5)
    rows, cols = unwrapped_phase.shape

    if sigma > rows or sigma > cols:
        msg = f"{wavelength_cutoff = } too large for image."
        msg += f"Shape = {(rows, cols)}, and {pixel_spacing = }"
        raise ValueError(msg)

    # Work on a copy of the displacement field to avoid modifying the input
    displacement = np.nan_to_num(unwrapped_phase)

    # Take either nan or 0 pixels in `unwrapped_phase` to be nodata
    nodata_mask = displacement == 0
    in_bounds_pixels = np.logical_not(nodata_mask)

    # Convert bad pixels to NaN for GDAL fillnodata
    displacement[bad_pixel_mask] = np.nan

    # Calculate the number of pixels for max_distance based on wavelength
    max_distance_pixels = int((wavelength_cutoff / 2) / pixel_spacing)

    scratch_dir = Path(scratch_dir) if scratch_dir is not None else None

    # Create temporary directory in the specified location or system default
    with tempfile.TemporaryDirectory(dir=scratch_dir) as temp_dir:
        tmp_path = Path(temp_dir)
        # Create paths for temporary files in the temp directory
        temp_src = tmp_path / "src.tif"
        temp_dst = tmp_path / "filled.tif"

        # Save the array to a temporary GeoTIFF
        io.write_arr(
            arr=displacement,
            nodata=np.nan,
            output_name=temp_src,
            shape=(rows, cols),
            dtype=displacement.dtype,
        )

        # Fill nodata using GDAL
        gdal_fillnodata(
            src_filename=str(temp_src),
            dst_filename=str(temp_dst),
            max_distance=max_distance_pixels,
            smoothing_iterations=0,
            interpolation="nearest",
            quiet=True,
        )

        filled_data = io.load_gdal(temp_dst, masked=True).filled(0)

    # Apply Gaussian filter
    lowpass_filtered = fft.fft2(filled_data, workers=workers)
    lowpass_filtered = ndimage.fourier_gaussian(lowpass_filtered, sigma=sigma)
    # Make sure to only take the real part (ifft returns complex)
    lowpass_filtered = fft.ifft2(lowpass_filtered, workers=workers).real.astype(
        unwrapped_phase.dtype
    )

    filtered_ifg = (filled_data - lowpass_filtered) * in_bounds_pixels
    if fill_value is not None:
        good_pixel_mask = np.logical_not(bad_pixel_mask)
        total_valid_mask = in_bounds_pixels & good_pixel_mask
        return np.where(total_valid_mask, filtered_ifg, fill_value)
    else:
        return filtered_ifg

filter_rasters(unw_filenames, cor_filenames=None, conncomp_filenames=None, temporal_coherence_filename=None, wavelength_cutoff=25000, correlation_cutoff=0.5, output_dir=None, max_workers=4)

Filter a list of unwrapped interferogram files using a long-wavelength filter.

Parameters:

Name	Type	Description	Default
unw_filenames	list[Path]	List of paths to unwrapped interferogram files to be filtered.	required
cor_filenames	list[Path] | None	List of paths to correlation files Passing None skips filtering on correlation.	None
conncomp_filenames	list[Path] | None	List of paths to connected component files, filters any 0 labeled pixels. Passing None skips filtering on connected component labels.	None
temporal_coherence_filename	Path | None	Path to the temporal coherence file for masking. Passing None skips filtering on temporal coherence.	None
wavelength_cutoff	float	Spatial wavelength cutoff (in meters) for the filter. Default is 50,000 meters.	25000
correlation_cutoff	float	Threshold of correlation (if passing cor_filenames) to use to ignore pixels during filtering.	0.5
output_dir	Path | None	Directory to save the filtered results. If None, saves in the same location as inputs with .filt.tif extension.	None
max_workers	int	Number of parallel images to process. Default is 4.	4

Returns:

Type	Description
list[Path]	Output filtered rasters.

Notes

Remove long-wavelength components from each unwrapped interferogram. It can optionally use temporal coherence, correlation, and connected component information for masking.

Source code in src/dolphin/filtering.py

def filter_rasters(
    unw_filenames: list[Path],
    cor_filenames: list[Path] | None = None,
    conncomp_filenames: list[Path] | None = None,
    temporal_coherence_filename: Path | None = None,
    wavelength_cutoff: float = 25_000,
    correlation_cutoff: float = 0.5,
    output_dir: Path | None = None,
    max_workers: int = 4,
) -> list[Path]:
    """Filter a list of unwrapped interferogram files using a long-wavelength filter.

    Parameters
    ----------
    unw_filenames : list[Path]
        List of paths to unwrapped interferogram files to be filtered.
    cor_filenames : list[Path] | None
        List of paths to correlation files
        Passing None skips filtering on correlation.
    conncomp_filenames : list[Path] | None
        List of paths to connected component files, filters any 0 labeled pixels.
        Passing None skips filtering on connected component labels.
    temporal_coherence_filename : Path | None
        Path to the temporal coherence file for masking.
        Passing None skips filtering on temporal coherence.
    wavelength_cutoff : float, optional
        Spatial wavelength cutoff (in meters) for the filter. Default is 50,000 meters.
    correlation_cutoff : float, optional
        Threshold of correlation (if passing `cor_filenames`) to use to ignore pixels
        during filtering.
    output_dir : Path | None, optional
        Directory to save the filtered results.
        If None, saves in the same location as inputs with .filt.tif extension.
    max_workers : int, optional
        Number of parallel images to process. Default is 4.

    Returns
    -------
    list[Path]
        Output filtered rasters.

    Notes
    -----
    Remove long-wavelength components from each unwrapped interferogram.
    It can optionally use temporal coherence, correlation, and connected component
    information for masking.


    """
    from dolphin import io

    bad_pixel_mask = np.zeros(
        io.get_raster_xysize(unw_filenames[0])[::-1], dtype="bool"
    )
    if temporal_coherence_filename:
        bad_pixel_mask = bad_pixel_mask | (
            io.load_gdal(temporal_coherence_filename) < 0.5
        )

    if output_dir is None:
        assert unw_filenames
        output_dir = unw_filenames[0].parent
    output_dir.mkdir(exist_ok=True)
    ctx = mp.get_context("spawn")

    with ProcessPoolExecutor(max_workers, mp_context=ctx) as pool:
        return list(
            pool.map(
                _filter_and_save,
                unw_filenames,
                cor_filenames or repeat(None),
                conncomp_filenames or repeat(None),
                repeat(output_dir),
                repeat(wavelength_cutoff),
                repeat(bad_pixel_mask),
                repeat(correlation_cutoff),
            )
        )

gaussian_filter_nan(image, sigma, mode='constant', **kwargs)

Apply a gaussian filter to an image with NaNs (avoiding all nans).

The scipy.ndimage gaussian_filter will make the output all NaNs if any of the pixels in the input that touches the kernel is NaN

Source: https://stackoverflow.com/a/36307291

Parameters:

Name	Type	Description	Default
image	ndarray	Image with nans to filter	required
sigma	float	Size of filter kernel. passed into gaussian_filter	required
mode	str	Boundary mode for [scipy.ndimage.gaussian_filter][]	= "constant"
**kwargs	Any	Passed into [scipy.ndimage.gaussian_filter][]	{}

Returns:

Type	Description
ndarray	Filtered version of image.

Source code in src/dolphin/filtering.py

def gaussian_filter_nan(
    image: ArrayLike, sigma: float | Sequence[float], mode="constant", **kwargs
) -> np.ndarray:
    """Apply a gaussian filter to an image with NaNs (avoiding all nans).

    The scipy.ndimage `gaussian_filter` will make the output all NaNs if
    any of the pixels in the input that touches the kernel is NaN

    Source:
    https://stackoverflow.com/a/36307291

    Parameters
    ----------
    image : ndarray
        Image with nans to filter
    sigma : float
        Size of filter kernel. passed into `gaussian_filter`
    mode : str, default = "constant"
        Boundary mode for `[scipy.ndimage.gaussian_filter][]`
    **kwargs : Any
        Passed into `[scipy.ndimage.gaussian_filter][]`

    Returns
    -------
    ndarray
        Filtered version of `image`.

    """
    from scipy.ndimage import gaussian_filter

    if np.sum(np.isnan(image)) == 0:
        return gaussian_filter(image, sigma=sigma, mode=mode, **kwargs)

    V = image.copy()
    nan_idxs = np.isnan(image)
    V[nan_idxs] = 0
    V_filt = gaussian_filter(V, sigma, **kwargs)

    W = np.ones(image.shape)
    W[nan_idxs] = 0
    W_filt = gaussian_filter(W, sigma, **kwargs)

    return V_filt / W_filt