diff --git a/python/packages/nisar/antenna/pattern.py b/python/packages/nisar/antenna/pattern.py index 9b6b65b44..0ab3acee5 100644 --- a/python/packages/nisar/antenna/pattern.py +++ b/python/packages/nisar/antenna/pattern.py @@ -6,7 +6,7 @@ from nisar.mixed_mode.logic import PolChannelSet from nisar.products.readers.antenna import AntennaParser from nisar.products.readers.instrument import InstrumentParser -from nisar.products.readers.Raw import Raw +from nisar.products.readers.Raw import Raw, chirpcorrelator_caltype_from_raw from nisar.antenna import TxTrmInfo, RxTrmInfo, TxBMF, RxDBF from nisar.antenna.beamformer import get_pulse_index import numpy as np @@ -114,8 +114,10 @@ def build_tx_trm(raw: Raw, pulse_times: np.ndarray, freq_band: str, """Build TxTrmInfo object """ # Parse Tx-related Cal stuff used for Tx BMF tx_chanl = raw.getListOfTxTRMs(freq_band, tx_pol) - corr_tap2 = raw.getChirpCorrelator(freq_band, tx_pol)[..., 1] - cal_type = raw.getCalType(freq_band, tx_pol) + # get chirp correlator and cal type for co-pol product + chp_corr, cal_type = chirpcorrelator_caltype_from_raw( + raw, txrx_pol=2 * tx_pol) + corr_tap2 = chp_corr[..., 1] # build TxTRM from Tx Cal stuff w/o optional "tx_phase" return TxTrmInfo(pulse_times, tx_chanl, corr_tap2, cal_type) diff --git a/python/packages/nisar/noise/noise_estimation_from_raw.py b/python/packages/nisar/noise/noise_estimation_from_raw.py index d31df2855..fa24cb1b2 100644 --- a/python/packages/nisar/noise/noise_estimation_from_raw.py +++ b/python/packages/nisar/noise/noise_estimation_from_raw.py @@ -16,6 +16,12 @@ from nisar.antenna import get_calib_range_line_idx from nisar.log import set_logger from isce3.core import DateTime, TimeDelta +from nisar.products.readers.Raw import ( + chirpcorrelator_caltype_from_raw, + is_raw_quad_pol, + first_tx_pol_for_quad, + opposite_pol +) # Global Noise-related Constants # Min number of range bins recommended per noise range block @@ -203,7 +209,7 @@ def extract_noise_only_lines(raw, freq_band, txrx_pol, max_lines=18944): # special RCIDs (NISAR mode numbers) to be treated as noise-only product rcid_special = (1, 2, 3) # get noise-only range lines if any - cal_path_mask = raw.getCalType(freq_band, tx=txrx_pol[0]) + _, cal_path_mask = chirpcorrelator_caltype_from_raw(raw, txrx_pol=txrx_pol) _, _, _, noise_index = get_calib_range_line_idx(cal_path_mask) # check if it is a special case with RCID=1,2,3 where there is no # noise-only range line and TX=OFF. @@ -456,7 +462,11 @@ def est_noise_power_from_raw( f'Number of range blocks is smaller than min {n_rg_blk_min}' ) logger.info(f'Number of range blocks -> {num_rng_block}') - + # check if product is quad + is_quad_pol = is_raw_quad_pol(raw) + if is_quad_pol: + first_tx_pol = first_tx_pol_for_quad(raw) + logger.info(f'Quad pol product w/ first {first_tx_pol}-pol TX!') # container for all noise products noise_prods = [] # if quad pol, then do MVE or MEE @@ -467,8 +477,9 @@ def est_noise_power_from_raw( # L-band NISAR. # loop over freq bands for freq_band in frq_pol: - # check if it is QP and product differentiation is set to True - if dif_quad and _is_quad_pol(frq_pol[freq_band]): + # check if it contains all linear pol combination and + # product differentiation is set to True + if dif_quad and _contains_all_linear_pols(frq_pol[freq_band]): logger.info('The difference of co-pol and cx-pol with' ' the same RX pol will be used in Noise est!') # let's combine datasets with the same RX Pol @@ -559,7 +570,9 @@ def est_noise_power_from_raw( ) noise_prods.append(ns_prod) - else: # other pol types than QP + else: # no polarimetric diff! + # For qaud pol, use noise range lines of the + # first TX pol for the other TX pol. for txrx_pol in frq_pol[freq_band]: logger.info( 'Processing individually frequency band ' @@ -575,12 +588,20 @@ def est_noise_power_from_raw( # calculate approximate ENBW for relatively white noise! enbw = enbw_from_raw(raw, freq_band, txrx_pol[0]) logger.info(f'Approximate ENBW in (MHz) -> {enbw * 1e-6}') + # check if quad pol per telemetry then use the oppsoite + # TX pol for noise range lines if that pol is not the + # first TX pol. + txrx_p = txrx_pol + if is_quad_pol and txrx_pol[0] != first_tx_pol: + txrx_p = opposite_pol(txrx_pol[0]) + txrx_pol[1] + logger.warning(f'Use noise-only range lines from {txrx_p} ' + f'for {txrx_pol}!') # parse one noise dataset # loop over several AZ blocks of noise-only range lines noise_power_azblk = [] az_dt_utc = [] for (dset_noise, idx_rgl_ns) in extract_noise_only_lines( - raw, freq_band, txrx_pol, max_lines): + raw, freq_band, txrx_p, max_lines): if exclude_first_last: logger.info( 'Exclude the first and last noise range lines.') @@ -641,10 +662,10 @@ def _pow2db(p: float) -> float: return 10 * np.log10(p) -def _is_quad_pol(txrx_pols): +def _contains_all_linear_pols(txrx_pols): """ - Whether the list of two-char TxRx Pols represents linear quad - polarization or not. + Whether the list of two-char TxRx Pols represents all + combinations of linear polarizations or not. Parameters ---------- diff --git a/python/packages/nisar/products/readers/Raw/Raw.py b/python/packages/nisar/products/readers/Raw/Raw.py index c68113dbb..16c68f9c6 100644 --- a/python/packages/nisar/products/readers/Raw/Raw.py +++ b/python/packages/nisar/products/readers/Raw/Raw.py @@ -1,3 +1,4 @@ +from __future__ import annotations from .DataDecoder import DataDecoder import h5py import isce3 @@ -9,8 +10,9 @@ import journal import re from warnings import warn -from scipy.interpolate import interp1d -from nisar.antenna import CalPath +from typing import Tuple +from enum import IntEnum, unique +from nisar.antenna import CalPath, get_calib_range_line_idx # TODO some CSV logger log = logging.getLogger("Raw") @@ -399,6 +401,140 @@ def getRangeLineIndex(self, frequency: str = 'A', tx: str = None): return f[path]["rangeLineIndex"][()] + def _parse_chirpcorrelator_from_hrt_qfsp( + self, + txrx_pol: str) -> np.ndarray | None: + """ + Parse three-tap chirp correlator array with shape (lines, 12, 3) + as well ass cal type with shape (lines,) from HRT QFSP. + + Parameters + ---------- + txrx_pol : str + TxRx polarization such as HH, VH, etc + + Returns + ------- + np.ndarray(complex) or None + 3-D complex array of chirp correlator with shape (Lines, channels, 3) + If the field does not exist None will be returned. + + """ + # get HRT path + hrt_path = self.TelemetryPath.replace('low', 'high') + qfsp_path = f'{hrt_path}/tx{txrx_pol[0]}/rx{txrx_pol[1]}/QFSP' + with h5py.File(self.filename, mode='r', swmr=True) as f5: + # loop over three qfsp + for i_qfsp in range(3): + p_qfsp = f'{qfsp_path}{i_qfsp}' + # loop over 4 channels per qfsp: + for nn in range(4): + i_chn = nn + i_qfsp * 4 + n_rx = i_chn + 1 + # loop over 3 taps per channel + for i_tap in range(3): + n_tap = i_tap + 1 + # form the path to the dataset per I and Q + # use RX pol! + p_ds_i = (f'{p_qfsp}/CHIRP_CORRELATOR_I{n_tap}_' + f'{txrx_pol[1]}{n_rx:02d}') + p_ds_q = (f'{p_qfsp}/CHIRP_CORRELATOR_Q{n_tap}_' + f'{txrx_pol[1]}{n_rx:02d}') + try: + ds_i = f5[p_ds_i] + except KeyError as err: + warn( + f'Missing dataset {p_ds_i} in {self.filename}.' + f' Detailed error -> {err}' + ) + return None + else: + # initialize the 3-D array, lines by 12 by 3 + if i_qfsp == nn == i_tap == 0: + # initialize the 3-D array for chirp correlator + num_lines = ds_i.size + chp_cor = np.ones((num_lines, 12, 3), dtype='c8') + chp_cor[:, i_chn, i_tap].real = ds_i[()] + chp_cor[:, i_chn, i_tap].imag = f5[p_ds_q][()] + return chp_cor + + + def _parse_caltype_from_hrt_qfsp( + self, + txrx_pol: str) -> np.ndarray | None: + """ + Parse cal type with shape (lines,) from HRT QFSP. + + Parameters + ---------- + txrx_pol : str + TxRx polarization such as HH, VH, etc + + Returns + ------- + np.ndarray(uint8) or None + 1-D array of cal type w/ values HPA=0, LNA=1, BYPASS=2, and + INVALID=255. If the field does not exist None will be returned. + + """ + # get HRT path + hrt_path = self.TelemetryPath.replace('low', 'high') + qfsp_path = f'{hrt_path}/tx{txrx_pol[0]}/rx{txrx_pol[1]}/QFSP' + with h5py.File(self.filename, mode='r', swmr=True) as f5: + # XXX get caltype from the very first qFSP assuming + # it is qFSP independent! + i_qfsp = 0 + p_qfsp = f'{qfsp_path}{i_qfsp}' + p_type = f'{p_qfsp}/CP_CAL_TYPE_{txrx_pol[1]}{i_qfsp}' + # XXX Following Try/exception block is added to + # support old sim L0B products lacking HRT! + try: + ds_cal_type = f5[p_type] + except KeyError as err: + warn(f'Missing dataset "{p_type}" in ' + f'"{self.filename}". Detailed error -> {err}') + return None + else: + return ds_cal_type[()].astype(CalPath) + + + def _parse_rangeline_index_from_hrt( + self, + txrx_pol: str = None) -> np.ndarray | None: + """ + Get range line index over all range lines from + HRT if exists otherwise None! + + Parameters + ---------- + txrx_pol : str + TxRx polarization such as HH, VH, etc + + Returns + ------- + np.ndarray(uint) or None + If not available in L0b, None will be returned. + + """ + hrt_path = self.TelemetryPath.replace('low', 'high') + freq_band = sorted(self.frequencies)[0] + pols = self.polarizations[freq_band] + if txrx_pol is None: + txrx_pol = pols[0] + elif txrx_pol not in pols: + raise ValueError(f'Available pols {pols} but got {txrx_pol}!') + rgl_idx_path = (f'{hrt_path}/tx{txrx_pol[0]}/rx{txrx_pol[1]}/' + 'RangeLine/RH_RANGELINE_INDEX') + with h5py.File(self.filename, mode='r', swmr=True) as f5: + try: + ds_rgl_idx = f5[rgl_idx_path] + except KeyError as err: + warn(f'Can not parse range line index from HRT. Error -> {err}') + return None + else: + return ds_rgl_idx[()] + + def getCalType(self, frequency: str = 'A', tx: str = None): """ Extract Tx Calibration mask for each range line. @@ -918,3 +1054,131 @@ def open_rrsd(filename) -> RawBase: if "/science/LSAR/RRSD/telemetry" in f: return LegacyRaw(hdf5file=filename) return Raw(hdf5file=filename) + + + +@unique +class PolarizationTypeId(IntEnum): + """Enumeration for polarization types of L-band NISAR""" + single_h = 0 + single_v = 1 + dual_h = 2 + dual_v = 3 + quad = 4 + compact = 5 + none = 6 + quasi_quad = 7 + quasi_dual = 8 + + +def polarization_type_from_drt(raw: Raw) -> PolarizationTypeId: + """Get polarization ID and type from L0B DRT""" + pol_path = f'{raw.TelemetryPath}/DRT/MISC/CP_IFSW_POLARIZATION' + with h5py.File(raw.filename, mode='r', swmr=True) as f5: + try: + ds_pol = f5[pol_path] + except KeyError: + warn(f'Missing dataset "{pol_path}" in "{raw.filename}"') + id_pol = 6 + else: + i_pol = ds_pol[()] + id_pol = np.nanmedian(i_pol) + return PolarizationTypeId(id_pol) + + +def is_raw_quad_pol(raw: Raw) -> bool: + """Determine whether raw L0B is Quad or not""" + return polarization_type_from_drt(raw) == PolarizationTypeId.quad + + +def first_tx_pol_for_quad(raw: Raw) -> str: + """Get first TX polarization, H or V, from only Quad pol product""" + if not is_raw_quad_pol(raw): + raise ValueError('Not a quad pol!') + idx_rgl = raw._parse_rangeline_index_from_hrt()[0] + # if not in HRT parse single-pol version from swath path + if idx_rgl is None: + idx_rgl_h = raw.getRangeLineIndex('A', 'H')[0] + idx_rgl_v = raw.getRangeLineIndex('A', 'V')[0] + if idx_rgl_v < idx_rgl_h: + return 'V' + return 'H' + else: # odd range line is V pol first and even is H pol first! + return {0: 'H', 1: 'V'}.get(idx_rgl % 2) + + +def opposite_pol(pol: str) -> str: + """Get the oppsoite pol""" + if pol == 'H': + return 'V' + elif pol == 'V': + return 'H' + else: + return pol + + +def chirpcorrelator_caltype_from_raw( + raw: Raw, + txrx_pol: str +) -> Tuple[np.ndarray, np.ndarray]: + """ + Parse three-tap chirp correlator array with shape (lines, 12, 3) + as well ass cal type with shape (lines,) from Raw L0B for a certain + TxRX pol + + Parameters + ---------- + raw : nisar.products.readers.Raw + txrx_pol : str + TxRx polarization such as HH, VH, etc + + Returns + ------- + np.ndarray(complex) + 3-D complex array of chirp correlator with shape (Lines, channels, 3) + np.ndarray(uint8) + 1-D array of cal type w/ values HPA=0, LNA=1, BYPASS=2, and INVALID=255 + + """ + chp_cor = raw._parse_chirpcorrelator_from_hrt_qfsp(txrx_pol=txrx_pol) + cal_type = raw._parse_caltype_from_hrt_qfsp(txrx_pol=txrx_pol) + # XXX if the respective field does not exist then use co-pol under + # swath in L0B for the sake of backward compatibility + if chp_cor is None or cal_type is None: + freq_band = [f for f in raw.frequencies if + txrx_pol in raw.polarizations[f]][0] + chp_cor = raw.getChirpCorrelator(freq_band, txrx_pol[0]) + cal_type = raw.getCalType(freq_band, txrx_pol[0]) + return chp_cor, cal_type + # Quad pol case + if is_raw_quad_pol(raw): + tx_pol_first = first_tx_pol_for_quad(raw) + if txrx_pol[0] == tx_pol_first: + chp_cor = chp_cor[::2] + cal_type = cal_type[::2] + else: # the second TX pol + # get data from the opssoite TX pol + x_pol = opposite_pol(txrx_pol[0]) + txrx_pol[1] + chp_cor_x, cal_type_x = chirpcorrelator_caltype_from_raw( + raw, txrx_pol=x_pol) + # if co-pol get HPA value from same TX but + # fill in LNA/BYP from oppsoite TX + if txrx_pol[0] == txrx_pol[1]: + chp_cor = chp_cor[1::2] + cal_type = cal_type[1::2] + _, idx_byp, idx_lna, _ = get_calib_range_line_idx(cal_type_x) + chp_cor[idx_byp] = chp_cor_x[idx_byp] + chp_cor[idx_lna] = chp_cor_x[idx_lna] + cal_type[idx_byp] = CalPath.BYPASS + cal_type[idx_lna] = CalPath.LNA + else: # x-pol product + chp_cor = chp_cor_x + cal_type = cal_type_x + # set x-pol HPA to INVALID given they are the mix of + # LNA from co-pol and HPA from x-pol! + if txrx_pol in ('HV', 'VH'): + idx_hpa, _, _, _ = get_calib_range_line_idx(cal_type) + if idx_hpa.size > 0: + warn(f'Set HPA cal type for x-pol {txrx_pol} to INVALID!') + cal_type[idx_hpa] = CalPath.INVALID + return chp_cor, cal_type \ No newline at end of file diff --git a/python/packages/nisar/products/readers/Raw/__init__.py b/python/packages/nisar/products/readers/Raw/__init__.py index 67b7cb788..a073287e7 100644 --- a/python/packages/nisar/products/readers/Raw/__init__.py +++ b/python/packages/nisar/products/readers/Raw/__init__.py @@ -1,2 +1,10 @@ -from .Raw import Raw, open_rrsd +from .Raw import ( + Raw, + open_rrsd, + chirpcorrelator_caltype_from_raw, + PolarizationTypeId, + is_raw_quad_pol, + first_tx_pol_for_quad, + opposite_pol +) from .DataDecoder import complex32, DataDecoder diff --git a/python/packages/nisar/workflows/focus.py b/python/packages/nisar/workflows/focus.py index dde3c990a..d66f25d23 100644 --- a/python/packages/nisar/workflows/focus.py +++ b/python/packages/nisar/workflows/focus.py @@ -16,7 +16,14 @@ find_overlapping_channel) from nisar.products.readers.antenna import AntennaParser from nisar.products.readers.instrument import InstrumentParser -from nisar.products.readers.Raw import Raw, open_rrsd +from nisar.products.readers.Raw import ( + Raw, + open_rrsd, + chirpcorrelator_caltype_from_raw, + is_raw_quad_pol, + first_tx_pol_for_quad, + opposite_pol +) from nisar.products.readers.rslc_cal import (RslcCalibration, parse_rslc_calibration, get_scale_and_delay, check_cal_validity_dates) from nisar.products.writers import SLC @@ -1931,9 +1938,10 @@ def temp(suffix): # Compute NESZ if there exist noise-only range lines # get noise only range line indexes within processing interval - cal_path_mask = raw.getCalType( - channel_in.freq_id, pol[0])[pulse_begin:pulse_end] - _, _, _, idx_noise = get_calib_range_line_idx(cal_path_mask) + _, cal_path_mask = chirpcorrelator_caltype_from_raw( + raw, txrx_pol=pol) + _, _, _, idx_noise = get_calib_range_line_idx( + cal_path_mask[pulse_begin:pulse_end]) # form output slant range vector for all noise products if cfg.processing.noise_equivalent_backscatter.fill_nan_ends: @@ -1962,7 +1970,29 @@ def temp(suffix): data_noise = np.memmap( fid_noise, mode='w+', shape=(nrgl_noise, rc.output_size), dtype=np.complex64) - rc.rangecompress(data_noise, raw_clean[idx_noise]) + # Check if raw is quad pol and the TX pol is not the + # first TX pol. Then extract noise only range line + # from the opposite TX pol w/ the same RX pol. + # XXX No RFI/caltone clean up of noise-only range lines + # for second TX pol products of quad pol! + raw_ns = np.copy(raw_clean[idx_noise]) + if is_raw_quad_pol(raw): + first_tx_pol = first_tx_pol_for_quad(raw) + log.info(f'Quad pol w/ first {first_tx_pol} pol!') + if pol[0] != first_tx_pol: + pol_ns = opposite_pol(pol[0]) + pol[1] + log.warning('Get noise-only range lines from ' + f'{pol_ns} for {pol} of quad pol!') + ds_ns = raw.getRawDataset(channel_in.freq_id, pol_ns) + idx_ns = np.arange(pulse_begin, pulse_end)[idx_noise] + # decode simply noise-only range lines and + # thus no need for memmap + raw_ns = ds_ns[idx_ns] + raw_ns *= bb_phasor[idx_noise, np.newaxis] + raw_ns[np.isnan(raw_ns)] = 0.0 + if cfg.processing.zero_fill_gaps: + fill_gaps(raw_ns, swaths[:, idx_noise, :], 0.0) + rc.rangecompress(data_noise, raw_ns) # build and apply antenna pattern correction for noise # pulses if EAP is True if cfg.processing.is_enabled.eap: