diff --git a/python/packages/isce3/core/__init__.py b/python/packages/isce3/core/__init__.py index bf99adf26..79dcce20c 100644 --- a/python/packages/isce3/core/__init__.py +++ b/python/packages/isce3/core/__init__.py @@ -19,3 +19,4 @@ from .projections import is_utm from .serialization import load_orbit_from_h5_group from . import types +from .decibel import pow2db, amp2db, abs2 diff --git a/python/packages/isce3/core/decibel.py b/python/packages/isce3/core/decibel.py new file mode 100644 index 000000000..41647dfd4 --- /dev/null +++ b/python/packages/isce3/core/decibel.py @@ -0,0 +1,10 @@ +import numpy as np + +def pow2db(x): + return 10 * np.log10(x) + +def abs2(z): + return z.real**2 + z.imag**2 + +def amp2db(z): + return pow2db(abs2(z)) \ No newline at end of file diff --git a/python/packages/nisar/products/readers/SLC/RSLC.py b/python/packages/nisar/products/readers/SLC/RSLC.py index b4f882cb0..dcac91809 100644 --- a/python/packages/nisar/products/readers/SLC/RSLC.py +++ b/python/packages/nisar/products/readers/SLC/RSLC.py @@ -1,6 +1,7 @@ # -*- coding: utf-8 -*- from __future__ import annotations +from functools import cached_property import h5py import journal import logging @@ -389,6 +390,51 @@ def getResampledNoiseEquivalentBackscatter( noise_product.txrx_pol) + @cached_property + def rangeChirpWeighting(self): + """ + Get the range spectral weights. + + Returns + ------- + values : numpy.ndarray + Expected shape of amplitude spectrum in range. Typically 256 + frequency bins, shifted so that the carrier frequency is in the + middle. + name : str + Name of the weighting function. + shape : float + Shape parameter of the window function. + """ + path = _h5join(self.ProcessingInformationPath, "parameters", + "rangeChirpWeighting") + with h5py.File(self.filename, 'r', libver='latest', swmr=True) as h5: + dset = h5[path] + name = dset.attrs["window_name"].decode() + shape = float(dset.attrs["window_shape"]) + values = dset[:] + return values, name, shape + + + @cached_property + def azimuthChirpWeighting(self): + """ + Get the azimuth spectral weights (antenna pattern). + + Returns + ------- + values : numpy.ndarray + Expected shape of amplitude spectrum in azimuth. Typically 256 + frequency bins, shifted so that the Doppler centroid is in the + middle. + """ + path = _h5join(self.ProcessingInformationPath, "parameters", + "azimuthChirpWeighting") + with h5py.File(self.filename, 'r', libver='latest', swmr=True) as h5: + values = h5[path][:] + return values + + def _h5join(*paths: str) -> str: """Join two paths to be used in HDF5""" # avoid repeated path separators diff --git a/python/packages/nisar/workflows/estimate_abscal_factor.py b/python/packages/nisar/workflows/estimate_abscal_factor.py index b9936df6e..3bbfbf549 100644 --- a/python/packages/nisar/workflows/estimate_abscal_factor.py +++ b/python/packages/nisar/workflows/estimate_abscal_factor.py @@ -3,19 +3,127 @@ import argparse import json +import logging import os import traceback import warnings -from collections.abc import Iterable +from collections.abc import Iterable, Callable from pathlib import Path from typing import Any, List, Optional, Union import numpy as np import shapely +from scipy.integrate import quad +from scipy.optimize import root_scalar import isce3 +from isce3.core import abs2 import nisar +log = logging.getLogger("nisar.workflows.estimate_abscal_factor") + + +def make_irf(weights, normalize=True): + """ + Create an impulse response function (IRF) from frequency domain weights. + + Parameters + ---------- + weights : array_like + Frequency domain weights (e.g., window or filter coefficients). + normalize : bool, optional + If True, normalize the weights to sum to 1. Defaults to True. + + Returns + ------- + callable + A function that evaluates the impulse response at a given time offset t. + """ + w = np.fft.fftshift(weights) / len(weights) + if normalize: + w *= 1.0 / np.sum(w) + f = np.fft.fftfreq(len(w)) + return lambda t: np.exp(1j * 2 * np.pi * f * t).dot(w) + +def get_irf_width_area(irf: Callable[[float], float], t_max=np.inf): + """ + Compute the full width at half maximum (FWHM) and integrated area of an + impulse response. + + Parameters + ---------- + irf : callable + Impulse response function that takes a time offset and returns a + complex value. + t_max : float, optional + Maximum integration limit for computing the area. Defaults to infinity. + + Returns + ------- + width : float + Full width at half maximum (FWHM) of the impulse response power. + area : float + Integrated area under the impulse response power curve. + """ + # Find half-power width using bracketing root-finding algorithm. + hw = root_scalar(lambda t: abs2(irf(t)) - 0.5, x0=0.0, x1=1.0).root + # full width = 2 * half width + width = float(2 * hw) + + # Find area with numerical integration (quadrature). + # Integrate from [0, inf) and double the result. + area = 2 * quad(lambda t: abs2(irf(t)), 0, t_max, + limit=10_000, epsabs=1e-6)[0] + return width, area + +def get_window_correction(weights): + """ + Compute the radiometric correction factor for a windowing function. + + The correction factor accounts for the power loss due to apodization + (windowing) relative to a rectangular window. It is computed as the ratio of + the IRF mainlobe area density (area/width) for the windowed response to that + of a rectangular window. + + Parameters + ---------- + weights : array_like + Frequency domain weights representing the window function. + + Returns + ------- + float + Window correction factor (dimensionless, typically < 1). + """ + # Note that DTFT is periodic, so limit integration to one period. + t_max = len(weights) / 2 + width_win, area_win = get_irf_width_area(make_irf(weights), t_max=t_max) + width_box, area_box = get_irf_width_area(make_irf(weights > 0.0), t_max=t_max) + return area_win / width_win / (area_box / width_box) + +def get_abscal_correction(rslc): + """ + Compute the combined window correction factor for absolute calibration. + + Calculates the 2-D radiometric correction by multiplying the range and + azimuth window correction factors derived from the chirp weighting functions + in the RSLC product. + + Parameters + ---------- + rslc : nisar.products.readers.SLC + The input RSLC product containing chirp weighting information. + + Returns + ------- + float + Combined window correction factor for range and azimuth (dimensionless). + """ + weights_az = rslc.azimuthChirpWeighting + weights_rg, _, _ = rslc.rangeChirpWeighting + corr_az = get_window_correction(weights_az) + corr_rg = get_window_correction(weights_rg) + return corr_rg * corr_az class DateTimeEncoder(json.JSONEncoder): def default(self, obj): @@ -116,6 +224,11 @@ def estimate_abscal_factor( function. The total power is estimated by multiplying the peak power by the 3dB response widths in along-track and cross-track directions. + 'adjusted_box': + Similar to 'box' but adjusted to account for the effects of + apodization windows in range and azimuth. This should give results + more comparable with methods that integrate sidelobes (e.g., ESA). + 'integrated': Measures power using the integrated power method. The total power is measured by summing the power of bins whose power exceeds a predefined minimum power @@ -245,6 +358,13 @@ def estimate_abscal_factor( # Get platform attitude data. attitude = rslc.getAttitude() + window_correction = 1.0 + meas_power_method = power_method + if power_method == "adjusted_box": + meas_power_method = "box" + window_correction = get_abscal_correction(rslc) + log.info(f"Will adjust 'box' RCS estimates by {window_correction = }") + # Estimate the absolute calibration error (the ratio of the measured RCS to the # predicted RCS) for a single corner reflector. def estimate_abscal_error( @@ -269,11 +389,11 @@ def estimate_abscal_error( upsample_factor=upsample_factor, peak_find_domain=peak_find_domain, nfit=nfit, - power_method=power_method, + power_method=meas_power_method, pthresh=pthresh, ) - return measured_rcs / predicted_rcs + return measured_rcs / predicted_rcs * window_correction # Estimate the absolute radiometric calibration error of each corner reflector, and # format the results into an object that can be easily JSON-ified. @@ -457,7 +577,7 @@ def parse_cmdline_args() -> dict[str, Any]: parser.add_argument( "--power-method", type=str, - choices=["box", "integrated"], + choices=["box", "adjusted_box", "integrated"], default="box", help=( "The method for estimating the target signal power (rectangular box method" diff --git a/share/nisar/defaults/focus.yaml b/share/nisar/defaults/focus.yaml index 8bbac52dc..d88b755d6 100644 --- a/share/nisar/defaults/focus.yaml +++ b/share/nisar/defaults/focus.yaml @@ -963,12 +963,18 @@ runconfig: # multiplying the peak power by the 3dB response widths in # along-track and cross-track directions. + # 'adjusted_box': + # Similar to 'box' but adjusted to account for the effects of + # apodization windows in range and azimuth. This should give + # results more comparable with methods that integrate sidelobes + # (e.g., ESA). + # 'integrated': # Measures power using the integrated power method. The total # power is measured by summing the power of bins whose power # exceeds a predefined minimum power threshold. - # Default: box - power_method: box + # Default: adjusted_box + power_method: adjusted_box # The minimum power threshold, measured in dB below the # peak power, for estimating the target signal power using the diff --git a/share/nisar/schemas/focus.yaml b/share/nisar/schemas/focus.yaml index 76d5a64e6..1f0a14094 100644 --- a/share/nisar/schemas/focus.yaml +++ b/share/nisar/schemas/focus.yaml @@ -607,7 +607,7 @@ absolute_radiometric_calibration_options: upsample_factor: int(min=1, required=False) peak_find_domain: enum('time', 'freq', required=False) nfit: int(min=3, required=False) - power_method: enum('box', 'integrated', required=False) + power_method: enum('box', 'adjusted_box', 'integrated', required=False) power_threshold: num(required=False) point_target_analyzer_options: