PTA: Add chip to outputs and estimate signal-to-clutter ratio.

python/packages/isce3/cal/point_target_info.py

@@ -591,7 +591,7 @@ def analyze_point_target(
 
     chip, chip_min_i, chip_min_j = generate_chip_on_slc(slc, i, j, chipsize=chipsize)
 
-    return analyze_point_target_chip(
+    results = analyze_point_target_chip(
         chip=chip,
         chip_min_i=chip_min_i,
         chip_min_j=chip_min_j,
@@ -610,6 +610,21 @@ def analyze_point_target(
         pixel_spacing=pixel_spacing,
     )
 
+    # XXX: The current method of estimating SCR does not work for geolocated products
+    # which may have skewed sidelobes. Doing this would require resampling these
+    # products to acquire a chip usable by the SCR estimator, or else finding some new
+    # means of estimating and masking the locations of the side lobes on the chip.
+    if geo_heading is None:
+
+        scr = estimate_scr(
+            chip=chip,
+            peak_magnitude=results[0]["magnitude"],
+        )
+
+        results[0]["signal clutter ratio"] = scr
+
+    return results
+
 
 def generate_chip_on_slc(
     slc: DatasetReader,
@@ -846,6 +861,8 @@ def analyze_point_target_chip(
     return_dict = {
         "magnitude": np.abs(chipmax),
         "phase": np.angle(chipmax),
+        "chip magnitude": np.abs(chip).tolist(),
+        "chip phase": np.angle(chip).tolist(),
         "azimuth": {
             "ISLR": azimuth_islr_db,
             "PSLR": azimuth_pslr_db,
@@ -914,6 +931,92 @@ def get_slice_indices(n: int) -> np.ndarray:
     return return_dict, None
 
 
+def estimate_scr(
+    chip,
+    clutter_half_width: int = 12,
+    peak_magnitude: float | None = None,
+) -> float:
+    """
+    Estimate the signal-to-clutter ratio (SCR), in dB, for a chip of data around a point
+    target.
+
+    Parameters
+    ----------
+    chip : np.ndarray of complex64
+        The data chip.
+    clutter_half_width : int, optional
+        An integer greater than 3 and smaller than half the width or length of `chip`.
+        Defines half the region used for estimating SCR. Defaults to 12.
+    peak_magnitude : float, optional
+        The identified peak magnitude of the data, or None for a rough estimate.
+        Defaults to None.
+
+    Returns
+    -------
+    float
+        The chip SCR, in dB.
+    """
+    if clutter_half_width < 3:
+        raise ValueError("clutter_half_width must be a positive number greater than 2.")
+    if clutter_half_width > chip.shape[0] / 2 or clutter_half_width > chip.shape[1] / 2:
+        raise ValueError(
+            f"clutter_half_width of {clutter_half_width} must be less than half of "
+            f"chip dimensions {chip.shape[0]} and {chip.shape[1]}."
+        )
+
+    # If the peak magnitude was not passed in, estimate it.
+    if peak_magnitude is None:
+        peak_magnitude = np.nanmax(np.abs(chip))
+
+    # Get the location of the peak magnitude.
+    k = np.nanargmax(np.abs(chip))
+    ichip, jchip = np.unravel_index(k, chip.shape)
+
+    # Make sure that the point target is far enough from the edge of the chip that
+    # an estimation region can be selected.
+    if (
+        ichip < clutter_half_width
+        or jchip < clutter_half_width
+        or ichip > chip.shape[0] - clutter_half_width
+        or jchip > chip.shape[1] - clutter_half_width
+    ):
+        warn(
+            "PTA Warning: Detected peak location too close to edge of chip. SCR could "
+            "not be calculated. SCR value will be returned as NaN."
+        )
+        return np.nan
+
+    # Create an estimation chip that is the column and row of the peak magnitude plus
+    # clutter_half_width pixels in either direction, converted to units of linear power.
+    scr_est_chip = np.abs(
+        chip[
+            ichip - clutter_half_width : ichip + clutter_half_width + 1,
+            jchip - clutter_half_width : jchip + clutter_half_width + 1,
+        ]
+    ) ** 2
+
+    # Set the row and column of peak power, plus one row and column on either side, to
+    # 0. This creates a "+" on the image of zeroed pixels corresponding to the IRF peak
+    # and its side lobes. The remaining nonzero pixels are the clutter region.
+    scr_est_chip[clutter_half_width - 1:clutter_half_width + 2, :] = 0
+    scr_est_chip[:, clutter_half_width - 1:clutter_half_width + 2] = 0
+
+    # Calculate the number of pixels that were not zeroed out, which is four regions
+    # of clutter_half_width -1 pixels squared.
+    sampled_area = (clutter_half_width - 1) ** 2 * 4
+
+    # The clutter power is the total power of the clutter region divided by the total
+    # number of pixels constituting that region.
+    clutter = np.sum(scr_est_chip) / sampled_area
+
+    # Calculate the peak power divided by average clutter power to get signal to clutter
+    # ratio.
+    scr = peak_magnitude ** 2 / clutter
+
+    # Convert to dB and then return.
+    return 10 * np.log10(scr)
+
+
 def sample_geocoded_side_lobe(
     chip: np.ndarray,
     heading: float,
@@ -1031,11 +1134,20 @@ def tofloatvals(x):
 
     Modifies the dictionary in-place and returns None.
     """
+    def list2floats(my_list):
+        """
+        Convert a list of any dimensions into a list of floats. The list must either
+        contain only other lists or only values that can be converted into floats.
+        """
+        if all(isinstance(xi, list) for xi in my_list):
+            return [list2floats(xi) for xi in my_list]
+        return [float(xi) for xi in my_list]
+
     for k in x:
         if type(x[k]) == dict:
             tofloatvals(x[k])
         elif type(x[k]) == list:
-            x[k] = [float(xki) for xki in x[k]]
+            x[k] = list2floats(x[k])
         else:
             x[k] = float(x[k])
 

tests/python/packages/nisar/workflows/point_target_analysis.py

@@ -137,6 +137,9 @@ def test_nisar_csv():
             "survey_date",
             "validity",
             "velocity",
+            "chip magnitude",
+            "chip phase",
+            "signal clutter ratio",
         }
         assert set(cr_info.keys()) == expected_keys
 
@@ -246,6 +249,9 @@ def test_uavsar_csv():
             "phase",
             "range",
             "azimuth",
+            "chip magnitude",
+            "chip phase",
+            "signal clutter ratio",
         }
         assert set(cr_info.keys()) == expected_keys