From 1012e7c1248e37b5625c8b4cfd6977bd73805b7c Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Tue, 12 Mar 2024 21:07:20 -0400 Subject: [PATCH 1/5] add reference, convert float weights once --- docs/references.bib | 15 ++++++++++++++- src/dolphin/interpolation.py | 26 +++++++++++++------------- 2 files changed, 27 insertions(+), 14 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index c7557bcd5..67f70d1c8 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -55,8 +55,21 @@ @article{Chen2012IonosphericArtifactsSimultaneous keywords = {Advanced Land Observation Satellite (ALOS) PhasedArray L-band Synthetic Aperture Radar (PALSAR),Azimuth,California,correlation signatures,Delay,dual-frequency GPS carrier phase data,global positioning system,Global Positioning System,global positioning system (GPS),GPS observations,Hawaii,Iceland,InSAR images,interferograms,interferometric synthetic aperture radar,interpretability,Ionosphere,ionospheric artifacts,ionospheric delay,ionospheric electromagnetic wave propagation,ionospheric propagation delays,ionospheric total electron content,ionospheric variability,L-band InSAR data,L-band InSAR observations,L-band SAR interferometry,neutral atmospheric delays,phase artifacts,pixel misregistration,radar imaging,radar interferometry,radiowave propagation,SAR acquisition times,Satellites,Spaceborne radar,synthetic aperture length scales,synthetic aperture radar,terrain,total electron content (TEC)} } +@article{Chen2015PersistentScattererInterpolation, + title = {A Persistent Scatterer Interpolation for Retrieving Accurate Ground Deformation over {{InSAR-decorrelated}} Agricultural Fields}, + author = {Chen, Jingyi and Zebker, Howard A. and Knight, Rosemary}, + year = {2015}, + journal = {Geophysical Research Letters}, + volume = {42}, + number = {21}, + pages = {9294--9301}, + issn = {19448007}, + doi = {10.1002/2015GL065031}, + keywords = {decorrelation,groundwater,InSAR deformation map,persistent scatterers} +} + @article{Fornaro2015CAESARApproachBased, - title = {{{CAESAR}}: {{An Approach Based}} on {{Covariance Matrix Decomposition}} to {{Improve Multibaseline}}{\textendash}{{Multitemporal Interferometric SAR Processing}}}, + title = {{{CAESAR}}: {{An Approach Based}} on {{Covariance Matrix Decomposition}} to {{Improve Multibaseline}}--{{Multitemporal Interferometric SAR Processing}}}, shorttitle = {{{CAESAR}}}, author = {Fornaro, Gianfranco and Verde, Simona and Reale, Diego and Pauciullo, Antonio}, year = {2015}, diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py index 267ed2eed..ad4cc9833 100644 --- a/src/dolphin/interpolation.py +++ b/src/dolphin/interpolation.py @@ -1,5 +1,6 @@ import numba import numpy as np +from numpy.typing import ArrayLike from dolphin._log import get_log @@ -9,15 +10,20 @@ def interpolate( - ifg: np.ndarray, - weights: np.ndarray, + ifg: ArrayLike, + weights: ArrayLike, num_neighbors: int = 20, max_radius: int = 51, min_radius: int = 0, alpha: float = 0.75, weight_cutoff: float = 0.0, ) -> np.ndarray: - r"""Interpolation of persistent scatterers. + """Interpolate a complex interferogram based on pixel weights. + + Build upon persistent scatterer interpolation used in + [@Chen2015PersistentScattererInterpolation] and + [@Wang2022AccuratePersistentScatterer] by allowing floating-point weights + instead of 0/1 PS weights. Parameters ---------- @@ -58,12 +64,6 @@ def interpolate( interpolated interferogram with the same amplitude, but different wrapped phase at non-ps pixels. - References - ---------- - "A persistent scatterer interpolation for retrieving accurate ground - deformation over InSAR\-decorrelated agricultural fields" - Chen et al., 2015, https://doi.org/10.1002/2015GL065031 - """ nrow, ncol = weights.shape @@ -79,13 +79,13 @@ def interpolate( f"Range of values as weights, using weight_cutoff = {weight_cutoff}" ) + weights_float = np.clip(weights.astype(np.float32), 0, 1) # Ensure weights are between 0 and 1 - if np.any(weights.astype(np.float32) > 1): + if np.any(weights_float > 1): logger.warning("weights array has values greater than 1. Clipping to 1.") - if np.any(weights.astype(np.float32) < 0): + if np.any(weights_float < 0): logger.warning("weights array has negative values. Clipping to 0.") - # Make shared versions of the input arrays to avoid copying in each thread - weights_float = np.clip(weights.astype(np.float32), 0, 1) + weights_float = np.clip(weights_float, 0, 1) interpolated_ifg = np.zeros((nrow, ncol), dtype=np.complex64) From 97246065bcfd1e56e6a4c5911b47360e17ea7918 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Tue, 12 Mar 2024 21:10:49 -0400 Subject: [PATCH 2/5] remove the extra all check --- src/dolphin/interpolation.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py index ad4cc9833..269ff3814 100644 --- a/src/dolphin/interpolation.py +++ b/src/dolphin/interpolation.py @@ -67,18 +67,6 @@ def interpolate( """ nrow, ncol = weights.shape - if np.all( - np.logical_or( - np.logical_or(weights == 0, weights == 1), - np.logical_or(weights is True, weights is False), - ) - ): - logger.info("Binary weights, using PS-like interpolation.") - else: - logger.info( - f"Range of values as weights, using weight_cutoff = {weight_cutoff}" - ) - weights_float = np.clip(weights.astype(np.float32), 0, 1) # Ensure weights are between 0 and 1 if np.any(weights_float > 1): From 1102da17a187f8b015d4ef61e97da9d14ea7f067 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Tue, 12 Mar 2024 21:11:11 -0400 Subject: [PATCH 3/5] Remove pymp reference --- src/dolphin/interpolation.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py index 269ff3814..3086a9869 100644 --- a/src/dolphin/interpolation.py +++ b/src/dolphin/interpolation.py @@ -98,7 +98,6 @@ def _interp_loop( nrow, ncol = weights.shape nindices = len(indices) for r0 in numba.prange(nrow): - # convert linear idx to row, col for c0 in range(ncol): if weights[r0, c0] >= weight_cutoff: interpolated_ifg[r0, c0] = ifg[r0, c0] From 1ef74c5df742d9ac4b4efc3cc5025cea07d9cab5 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Tue, 12 Mar 2024 21:13:00 -0400 Subject: [PATCH 4/5] remove todo --- src/dolphin/interpolation.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py index 3086a9869..9ea073507 100644 --- a/src/dolphin/interpolation.py +++ b/src/dolphin/interpolation.py @@ -128,7 +128,8 @@ def _interp_loop( if counter >= num_neighbors: break - # TODO : why use the "counter - 1" here to normalize? + # `counter` got up to one more than the number of elements + # The last one will be the largest radius r2_norm = (r2[counter - 1] ** alpha) / 2 for i in range(counter): csum += np.exp(-r2[i] / r2_norm) * cphase[i] From ab26f93794b84ce1fb416d5a6c2adaf10b9ecdae Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Tue, 12 Mar 2024 21:54:10 -0400 Subject: [PATCH 5/5] Pass through the correlation threshold for interpolation The last implementation had a default of 0.0 Since every weight was >=0, no pixels were interpolated. --- src/dolphin/interpolation.py | 16 +++++++--------- src/dolphin/unwrap/_unwrap.py | 22 +++++++++++++++++++--- src/dolphin/workflows/config/_common.py | 9 +++++++++ src/dolphin/workflows/unwrapping.py | 1 + 4 files changed, 36 insertions(+), 12 deletions(-) diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py index 9ea073507..4162b6fe6 100644 --- a/src/dolphin/interpolation.py +++ b/src/dolphin/interpolation.py @@ -12,11 +12,11 @@ def interpolate( ifg: ArrayLike, weights: ArrayLike, + weight_cutoff: float = 0.5, num_neighbors: int = 20, max_radius: int = 51, min_radius: int = 0, alpha: float = 0.75, - weight_cutoff: float = 0.0, ) -> np.ndarray: """Interpolate a complex interferogram based on pixel weights. @@ -36,6 +36,12 @@ def interpolate( weights[i,j] = True if radar pixel (i,j) is a PS weights[i,j] = False if radar pixel (i,j) is not a PS Can also pass a coherence image to use as weights. + weight_cutoff: float + Threshold to use on `weights` so that pixels where + `weight[i, j] < weight_cutoff` have phase values replaced by + an interpolated value. + The default is 0.5: pixels with weight less than 0.5 are replaced with a + smoothed version of the surrounding pixels. num_neighbors: int (optional) number of nearest PS pixels used for interpolation num_neighbors = 20 by default @@ -49,14 +55,6 @@ def interpolate( hyperparameter controlling the weight of PS in interpolation: smaller alpha means more weight is assigned to PS closer to the center pixel. alpha = 0.75 by default - weight_cutoff: float - Threshold to use on `weights` so that pixels where - `weight[i, j] < weight_cutoff` have phase values replaced by - an interpolated value. - If `weight_cutoff = 0` (default), All pixels are replaced with a - smoothed version of the surrounding pixels. - If `weight_cutoff = 1`, only pixels with exactly weight=1 - are kept, and the rest are replaced with an interpolated value. Returns ------- diff --git a/src/dolphin/unwrap/_unwrap.py b/src/dolphin/unwrap/_unwrap.py index b96791db5..02074a4cf 100644 --- a/src/dolphin/unwrap/_unwrap.py +++ b/src/dolphin/unwrap/_unwrap.py @@ -51,6 +51,7 @@ def run( alpha: float = 0.5, run_interpolation: bool = False, max_radius: int = 51, + interpolation_cor_threshold: float = 0.5, ) -> tuple[list[Path], list[Path]]: """Run snaphu on all interferograms in a directory. @@ -109,6 +110,10 @@ def run( Whether to run interpolation on interferogram max_radius : int, optional, default = 51 maximum radius (in pixel) for scatterer searching for interpolation + interpolation_cor_threshold : float, optional, default = 0.5 + Threshold on the correlation raster to use for interpolation. + Pixels with less than this value are replaced by a weighted + combination of neighboring pixels. Returns ------- @@ -175,6 +180,7 @@ def run( alpha=alpha, run_interpolation=run_interpolation, max_radius=max_radius, + interpolation_cor_threshold=interpolation_cor_threshold, ) for ifg_file, out_file, cor_file in zip(in_files, out_files, cor_filenames) ] @@ -210,6 +216,7 @@ def unwrap( alpha: float = 0.5, run_interpolation: bool = False, max_radius: int = 51, + interpolation_cor_threshold: float = 0.5, ) -> tuple[Path, Path]: """Unwrap a single interferogram using snaphu, isce3, or tophu. @@ -272,6 +279,10 @@ def unwrap( Whether to run interpolation on interferogram max_radius : int, optional, default = 51 maximum radius (in pixel) for scatterer searching for interpolation + interpolation_cor_threshold : float, optional, default = 0.5 + Threshold on the correlation raster to use for interpolation. + Pixels with less than this value are replaced by a weighted + combination of neighboring pixels. Returns ------- @@ -354,11 +365,16 @@ def unwrap( ifg = io.load_gdal(ifg_filename) corr = io.load_gdal(corr_filename) - logger.info("Masking pixels with correlation above 0.7") - coherent_pixel_mask = corr[:] >= 0.7 + logger.info( + f"Masking pixels with correlation below {interpolation_cor_threshold}" + ) + coherent_pixel_mask = corr[:] >= interpolation_cor_threshold logger.info(f"Interpolating {ifg_filename} -> {interp_ifg_filename}") modified_ifg = interpolate( - ifg=ifg, weights=coherent_pixel_mask, max_radius=max_radius + ifg=ifg, + weights=coherent_pixel_mask, + weight_cutoff=interpolation_cor_threshold, + max_radius=max_radius, ) logger.info(f"Writing interpolated output to {interp_ifg_filename}") io.write_arr( diff --git a/src/dolphin/workflows/config/_common.py b/src/dolphin/workflows/config/_common.py index 6b951abec..aded917dc 100644 --- a/src/dolphin/workflows/config/_common.py +++ b/src/dolphin/workflows/config/_common.py @@ -222,8 +222,17 @@ class UnwrapOptions(BaseModel, extra="forbid"): ) max_radius: int = Field( 51, + ge=0.0, description=("(for interpolation) maximum radius to find scatterers."), ) + interpolation_cor_threshold: float = Field( + 0.5, + description=" Threshold on the correlation raster to use for interpolation. " + "Pixels with less than this value are replaced by a weighted " + "combination of neighboring pixels.", + ge=0.0, + le=1.0, + ) @field_validator("ntiles", "downsample_factor", mode="before") @classmethod diff --git a/src/dolphin/workflows/unwrapping.py b/src/dolphin/workflows/unwrapping.py index 0f34ab9c5..babd08dee 100644 --- a/src/dolphin/workflows/unwrapping.py +++ b/src/dolphin/workflows/unwrapping.py @@ -91,6 +91,7 @@ def run( alpha=unwrap_options.alpha, run_interpolation=unwrap_options.run_interpolation, max_radius=unwrap_options.max_radius, + interpolation_cor_threshold=unwrap_options.interpolation_cor_threshold, ) if add_overviews: