diff --git a/python/packages/isce3/signal/qfsp_slip_insar.py b/python/packages/isce3/signal/qfsp_slip_insar.py new file mode 100644 index 000000000..53a01656e --- /dev/null +++ b/python/packages/isce3/signal/qfsp_slip_insar.py @@ -0,0 +1,323 @@ +import numpy as np +import h5py +from scipy.ndimage import distance_transform_edt + + +def check_qfsp_flag(slc_path): + qFSP_path = '/science/LSAR/identification/hasInputDataException/' + with h5py.File(slc_path) as src: + qfsp_flag = src[qFSP_path][()] + return qfsp_flag + + +def _fit_2d_polynomial_surface(data, valid_mask, order=2): + """ + Fit a global 2D polynomial surface to data using valid_mask pixels only. + """ + nrows, ncols = data.shape + yy, xx = np.indices((nrows, ncols)) + + good = valid_mask & np.isfinite(data) + if np.count_nonzero(good) < 10: + raise ValueError("Not enough valid pixels for 2D polynomial fit.") + + x = xx[good].astype(float) + y = yy[good].astype(float) + z = data[good].astype(float) + + x_mean = x.mean() + y_mean = y.mean() + x_std = x.std() if x.std() > 0 else 1.0 + y_std = y.std() if y.std() > 0 else 1.0 + + xn = (x - x_mean) / x_std + yn = (y - y_mean) / y_std + + terms = [np.ones_like(xn)] + if order >= 1: + terms += [xn, yn] + if order >= 2: + terms += [xn**2, xn * yn, yn**2] + if order >= 3: + terms += [xn**3, (xn**2) * yn, xn * (yn**2), yn**3] + + A = np.vstack(terms).T + coeff, _, _, _ = np.linalg.lstsq(A, z, rcond=None) + + yy_full, xx_full = np.indices((nrows, ncols)) + xn_full = (xx_full.astype(float) - x_mean) / x_std + yn_full = (yy_full.astype(float) - y_mean) / y_std + + full_terms = [np.ones_like(xn_full)] + if order >= 1: + full_terms += [xn_full, yn_full] + if order >= 2: + full_terms += [xn_full**2, xn_full * yn_full, yn_full**2] + if order >= 3: + full_terms += [xn_full**3, (xn_full**2) * yn_full, xn_full * (yn_full**2), yn_full**3] + + surface = np.zeros((nrows, ncols), dtype=float) + for c, t in zip(coeff, full_terms): + surface += c * t + + return surface + + +def _find_column_groups(artifact_mask, min_fraction_rows=0.3, min_group_width=2): + """ + Find contiguous column groups affected in many rows based on artifact_mask. + """ + nrows, ncols = artifact_mask.shape + frac = artifact_mask.sum(axis=0) / max(nrows, 1) + affected_cols = frac >= min_fraction_rows + + groups = [] + in_group = False + start = None + + for j in range(ncols): + if affected_cols[j] and not in_group: + start = j + in_group = True + elif not affected_cols[j] and in_group: + end = j + if end - start >= min_group_width: + groups.append((start, end)) + in_group = False + + if in_group: + end = ncols + if end - start >= min_group_width: + groups.append((start, end)) + + return groups + + +def _moving_average_1d(x, win): + """ + NaN-aware moving average for 1D array. + """ + x = np.asarray(x, dtype=float) + if win <= 1: + return x.copy() + + valid = np.isfinite(x).astype(float) + x0 = np.where(np.isfinite(x), x, 0.0) + kernel = np.ones(win, dtype=float) + + num = np.convolve(x0, kernel, mode="same") + den = np.convolve(valid, kernel, mode="same") + + out = np.full_like(x, np.nan, dtype=float) + ok = den > 0 + out[ok] = num[ok] / den[ok] + return out + + +def make_feather_weight( + artifact_mask, + inner_shrink=2, + outer_feather=10, +): + """ + Create a smooth feathering weight for artifact correction. + + Parameters + ---------- + artifact_mask : 2D bool array + True where the artifact correction should be strongest. + + inner_shrink : int + Pixels inside the artifact edge where weight begins tapering. + Larger values make the full-strength region smaller. + + outer_feather : int + Number of pixels outside artifact_mask over which correction tapers to 0. + + Returns + ------- + weight : 2D float array + Smooth weight from 0 to 1. + """ + + artifact_mask = np.asarray(artifact_mask, dtype=bool) + + # Distance inside artifact region + dist_inside = distance_transform_edt(artifact_mask) + + # Distance outside artifact region + dist_outside = distance_transform_edt(~artifact_mask) + + weight = np.zeros(artifact_mask.shape, dtype=np.float32) + + # Inside artifact area: mostly full correction + inside = artifact_mask + weight[inside] = 1.0 + + # Taper near inner boundary + if inner_shrink > 0: + inner_edge = inside & (dist_inside <= inner_shrink) + weight[inner_edge] = dist_inside[inner_edge] / float(inner_shrink) + + # Outside artifact area: apply small tapered correction into neighborhood + if outer_feather > 0: + outside_feather = (~artifact_mask) & (dist_outside <= outer_feather) + weight[outside_feather] = ( + 1.0 - dist_outside[outside_feather] / float(outer_feather) + ) + + # Smoothstep curve: smoother transition than linear + weight = weight * weight * (3.0 - 2.0 * weight) + + return weight + + +def correct_qfsp_phase_artifact( + phase, + fit_background_mask, + artifact_mask, + background_order=2, + min_fraction_rows=0.3, + min_group_width=2, + template_smooth_win=3, + inner_shrink=2, + outer_feather=10, + fill_value=np.nan, +): + """ + Detrend phase, estimate smooth qFSP artifact from artifact pixels, + then subtract the smooth artifact from the original phase. + + Parameters + ---------- + phase : np.ndarray + 2D differential interferogram phase. + + fit_background_mask : np.ndarray + Valid clean pixels used to estimate the smooth background. + This should usually exclude artifact_mask. + + artifact_mask : np.ndarray + Pixels affected by qFSP artifact. The artifact template is estimated + from these pixels, and correction is applied to these pixels with + feathering near the boundary. + """ + + phase = np.asarray(phase, dtype=float) + fit_background_mask = np.asarray(fit_background_mask, dtype=bool) + artifact_mask = np.asarray(artifact_mask, dtype=bool) + + if phase.ndim != 2: + raise ValueError("phase must be 2D") + + nrows, ncols = phase.shape + + for name, arr in [ + ("fit_background_mask", fit_background_mask), + ("artifact_mask", artifact_mask), + ]: + if arr.shape != phase.shape: + raise ValueError(f"{name} must have same shape as phase") + + valid_phase = np.isfinite(phase) & (phase != 0) + + # Pixels used to fit the background: clean, valid, non-artifact pixels. + fit_mask = ( + fit_background_mask + & (~artifact_mask) + & valid_phase + ) + + # Pixels used to estimate the artifact template: valid artifact pixels. + template_estimation_mask = ( + artifact_mask + & valid_phase + ) + + background = _fit_2d_polynomial_surface( + phase, + fit_mask, + order=background_order, + ) + + residual = phase - background + + groups = _find_column_groups( + artifact_mask, + min_fraction_rows=min_fraction_rows, + min_group_width=min_group_width, + ) + + artifact_2d = np.full_like(phase, fill_value, dtype=float) + templates = [] + + for c0, c1 in groups: + group_residual = residual[:, c0:c1] + + valid_for_template = ( + template_estimation_mask[:, c0:c1] + & np.isfinite(group_residual) + ) + + tmp = np.where(valid_for_template, group_residual, np.nan) + template = np.nanmean(tmp, axis=0) + + if np.all(~np.isfinite(template)): + templates.append(template) + continue + + if template_smooth_win > 1: + template = _moving_average_1d(template, template_smooth_win) + + templates.append(template) + + # Extend the artifact model slightly outside the artifact columns + # so feathering can smoothly taper the correction. + c0_ext = max(0, c0 - outer_feather) + c1_ext = min(ncols, c1 + outer_feather) + + template_ext = np.interp( + np.arange(c0_ext, c1_ext), + np.arange(c0, c1), + template, + left=template[0], + right=template[-1], + ) + + artifact_group_ext = np.tile(template_ext[None, :], (nrows, 1)) + + valid_ext = np.isfinite(artifact_group_ext) + artifact_2d[:, c0_ext:c1_ext][valid_ext] = ( + artifact_group_ext[valid_ext] + ) + + correction_weight = make_feather_weight( + artifact_mask, + inner_shrink=inner_shrink, + outer_feather=outer_feather, + ) + + corrected_phase = phase.copy() + + valid_corr = ( + (correction_weight > 0) + & np.isfinite(artifact_2d) + & np.isfinite(corrected_phase) + ) + + corrected_phase[valid_corr] -= ( + correction_weight[valid_corr] * artifact_2d[valid_corr] + ) + + return { + "corrected_phase": corrected_phase, + "artifact_2d": artifact_2d, + "background": background, + "residual": residual, + "fit_mask": fit_mask, + "template_estimation_mask": template_estimation_mask, + "groups": groups, + "templates": templates, + "correction_weight": correction_weight, + "valid_corr": valid_corr, + } diff --git a/python/packages/isce3/unwrap/preprocess.py b/python/packages/isce3/unwrap/preprocess.py index ab88651e9..878f38b0c 100644 --- a/python/packages/isce3/unwrap/preprocess.py +++ b/python/packages/isce3/unwrap/preprocess.py @@ -644,3 +644,135 @@ def interpret_subswath_mask(subswath_mask, nodata=255): water = np.where(nd, False, water) return reference_valid, secondary_valid, water + + +def parse_insar_mask(mask): + """ + Parse the NISAR InSAR 2D mask generated by generate_insar_mask(). + + Encoding from generate_insar_mask() + ----------------------------------- + mask_id = ( + subswath_mask_id + | (sec_input_exception_mask << 8) + | (ref_input_exception_mask << 16) + ) + + bits 0-7: + subswath_mask_id + + This is produced by _compute_subswath_mask_id(). + In the current convention, it is decimal-style encoded: + + hundreds digit : water flag + tens digit : reference subswath ID + ones digit : secondary subswath ID + + Example: + 121 means: + water flag = 1 + reference subswath = 2 + secondary subswath = 1 + + 20 means: + water flag = 0 + reference subswath = 2 + secondary subswath = 0 + + bits 8-15: + secondary RSLC inputDataExceptionMask value. + + bits 16-23: + reference RSLC inputDataExceptionMask value. + + bit 24: + optional custom ionosphere/qFSP mask bit, if added later. + + Parameters + ---------- + mask : np.ndarray + 2D uint32 InSAR mask array. + + Returns + ------- + dict + Dictionary containing parsed 2D arrays. + """ + + mask = np.asarray(mask, dtype=np.uint32) + + if mask.ndim != 2: + raise ValueError(f"mask must be 2D. Got shape {mask.shape}") + + # bits 0-7: subswath mask id + subswath_mask_id = mask & np.uint32(0xFF) + + secondary_subswath = subswath_mask_id % np.uint32(10) + reference_subswath = ( + subswath_mask_id // np.uint32(10) + ) % np.uint32(10) + water_flag = ( + subswath_mask_id // np.uint32(100) + ) % np.uint32(10) + + secondary_valid = secondary_subswath != 0 + reference_valid = reference_subswath != 0 + valid_subswath = reference_valid & secondary_valid + + water = water_flag != 0 + + # bits 8-15: secondary inputDataExceptionMask + secondary_input_exception = ( + mask >> np.uint32(8) + ) & np.uint32(0xFF) + + # bits 16-23: reference inputDataExceptionMask + reference_input_exception = ( + mask >> np.uint32(16) + ) & np.uint32(0xFF) + + secondary_has_exception = secondary_input_exception != 0 + reference_has_exception = reference_input_exception != 0 + has_exception = secondary_has_exception | reference_has_exception + + # bit 24: optional custom ionosphere/qFSP mask + ionosphere_mask = ( + (mask >> np.uint32(24)) & np.uint32(1) + ) != 0 + + # bits 25-31: reserved + reserved = ( + mask >> np.uint32(25) + ) & np.uint32(0x7F) + + return { + # original + "mask": mask, + + # low byte / subswath information + "subswath_mask_id": subswath_mask_id, + "reference_subswath": reference_subswath, + "secondary_subswath": secondary_subswath, + "water_flag": water_flag, + "water": water, + + # validity from subswath id + "reference_valid": reference_valid, + "secondary_valid": secondary_valid, + "valid_subswath": valid_subswath, + + # inputDataExceptionMask values + "reference_input_exception": reference_input_exception, + "secondary_input_exception": secondary_input_exception, + + # exception booleans + "ref_has_anomaly": reference_has_exception, + "sec_has_anomaly": secondary_has_exception, + "has_exception": has_exception, + + # optional custom bit + "ionosphere_mask": ionosphere_mask, + + # reserved upper bits + "reserved": reserved, + } diff --git a/python/packages/nisar/products/insar/utils.py b/python/packages/nisar/products/insar/utils.py index 0643bde0c..578af736e 100644 --- a/python/packages/nisar/products/insar/utils.py +++ b/python/packages/nisar/products/insar/utils.py @@ -589,15 +589,18 @@ def _load_exception_mask(h5_obj, rslc_obj, swath): sec_subswaths) # reference RSLC input exception mask id - ref_input_exception_mask_id = ref_input_exception_mask[int(i),int(j)] << 16 + ref_input_exception_mask_id = int( + ref_input_exception_mask[int(i), int(j)] + ) << 16 # secondary RSLC input exception mask id sec_i = round(i + azimuth_off[0,int(j)]) sec_j = round(j + range_off[0,int(j)]) if ((sec_i >=0 and sec_i < sec_swath.lines) and (sec_j >=0 and sec_j < sec_swath.samples)): - sec_input_exception_mask_id = sec_input_exception_mask[sec_i,sec_j] << 8 - + sec_input_exception_mask_id = int( + sec_input_exception_mask[sec_i, sec_j] + ) << 8 # mask id mask_id = subswath_mask_id | ref_input_exception_mask_id | sec_input_exception_mask_id diff --git a/python/packages/nisar/workflows/ionosphere.py b/python/packages/nisar/workflows/ionosphere.py index ce7c218b9..70e0390bc 100644 --- a/python/packages/nisar/workflows/ionosphere.py +++ b/python/packages/nisar/workflows/ionosphere.py @@ -23,9 +23,12 @@ from isce3.io import HDF5OptimizedReader from isce3.signal.interpolate_by_range import (decimate_freq_a_array, interpolate_freq_b_array) +from isce3.signal.qfsp_slip_insar import (check_qfsp_flag, + correct_qfsp_phase_artifact) from isce3.splitspectrum import splitspectrum from isce3.unwrap.bridge_phase import bridge_unwrapped_phase -from isce3.unwrap.preprocess import project_map_to_radar +from isce3.unwrap.preprocess import (parse_insar_mask, + project_map_to_radar) from nisar.products.insar.product_paths import (CommonPaths, RIFGGroupsPaths, RUNWGroupsPaths) from nisar.products.readers import SLC @@ -1146,7 +1149,7 @@ def run(cfg: dict, runw_hdf5: str): filling_guide_median_size = filter_cfg['filling_guide_median_size'] filling_outlier_threshold = filter_cfg['filling_outlier_threshold'] filling_outlier_min_scale = filter_cfg['filling_outlier_min_scale'] - filling_outlier_mad_scale_factor = filter_cfg['filling_outlier_mad_scale_factor'] + filling_outlier_mad_scale_factor = filter_cfg['filling_outlier_mad_scale_factor'] filter_iterations = filter_cfg['filter_iterations'] median_filter_size = filter_cfg['median_filter_size'] @@ -1167,10 +1170,36 @@ def run(cfg: dict, runw_hdf5: str): unwrap_rg_looks = cfg['processing']['phase_unwrap']['range_looks'] unwrap_az_looks = cfg['processing']['phase_unwrap']['azimuth_looks'] + qfsp_cfg = iono_args.get("qfsp_correction", {}) + + iono_qfsp_correction_flag = qfsp_cfg.get("enabled", False) + + qfsp_background_order = qfsp_cfg.get("background_order", 2) + qfsp_min_fraction_rows = qfsp_cfg.get("min_fraction_rows", 0.05) + qfsp_min_group_width = qfsp_cfg.get("min_group_width", 1) + qfsp_template_smooth_win = qfsp_cfg.get("template_smooth_win", 3) + qfsp_inner_shrink = qfsp_cfg.get("inner_shrink", 2) + qfsp_outer_feather = qfsp_cfg.get("outer_feather", 10) + if unwrap_rg_looks != 1 or unwrap_az_looks != 1: rg_looks = unwrap_rg_looks az_looks = unwrap_az_looks + iono_qfsp_correction_flag = True # iono_args['qFSP_correction'] + ref_qfsp_flag = check_qfsp_flag(cfg["input_file_group"]["reference_rslc_file"]) + sec_qfsp_flag = check_qfsp_flag(cfg["input_file_group"]["secondary_rslc_file"]) + if iono_qfsp_correction_flag: + if not ref_qfsp_flag and not sec_qfsp_flag: + info_channel.log( + "qFSP correction is enabled but none of input " + "does not qFSP slip" + ) + iono_qfsp_correction_flag = False + else: + info_channel.log(f"reference qFSP: {ref_qfsp_flag}") + info_channel.log(f"secondary qFSP: {sec_qfsp_flag}") + need_insar_mask = iono_qfsp_correction_flag or ("subswath_mask" in mask_type) + # set paths for ionosphere and split spectrum iono_path = os.path.join(scratch_path, 'ionosphere') split_slc_path = os.path.join(iono_path, 'split_spectrum') @@ -1462,7 +1491,7 @@ def run(cfg: dict, runw_hdf5: str): [block_rows_data, cols_main], dtype=float) - if "subswath_mask" in mask_type: + if need_insar_mask: subswath_mask_image = np.empty( [block_rows_data, cols_main], dtype=int) @@ -1520,7 +1549,7 @@ def run(cfg: dict, runw_hdf5: str): sub_high_conn_image, np.s_[row_start:row_start + block_rows_data, :]) - if "subswath_mask" in mask_type: + if need_insar_mask: src_low_h5[subswath_mask_freq_a_path].read_direct( subswath_mask_image, np.s_[row_start:row_start + block_rows_data, :]) @@ -1629,7 +1658,7 @@ def run(cfg: dict, runw_hdf5: str): [block_rows_data, cols_side], dtype=float) - if "subswath_mask" in mask_type: + if need_insar_mask: subswath_mask_main_image = np.empty( [block_rows_data, cols_main], dtype=int) @@ -1678,7 +1707,7 @@ def run(cfg: dict, runw_hdf5: str): side_conn_image, np.s_[row_start:row_start + block_rows_data, :]) - if "subswath_mask" in mask_type: + if need_insar_mask: src_main_h5[subswath_mask_freq_a_path].read_direct( subswath_mask_main_image, np.s_[row_start:row_start + block_rows_data, :]) @@ -1772,6 +1801,168 @@ def run(cfg: dict, runw_hdf5: str): iterations=filter_iterations, filter_method='convolution') + if iono_qfsp_correction_flag or filter_bool: + + info_channel.log(f'{mask_type} is used for mask construction') + + available_mask = np.ones([block_rows_data, cols_output], dtype=bool) + + if "coherence" in mask_type: + mask_image = iono_phase_obj.get_coherence_mask_array( + main_array=main_coh_image, + side_array=side_coh_image, + diff_ms_array=diff_ms_coh_image, + low_band_array=sub_low_coh_image, + high_band_array=sub_high_coh_image, + diff_low_high_band_array=diff_coh_image, + slant_main=main_slant, + slant_side=side_slant, + threshold=filter_coh_thresh) + available_mask &= np.array(mask_image, dtype=bool) + + if "connected_components" in mask_type: + mask_image = iono_phase_obj.get_conn_component_mask_array( + main_array=main_conn_image, + side_array=side_conn_image, + diff_ms_array=diff_ms_conn_image, + low_band_array=sub_low_conn_image, + high_band_array=sub_high_conn_image, + diff_low_high_band_array=diff_subband_conn_image, + slant_main=main_slant, + slant_side=side_slant) + available_mask &= np.array(mask_image, dtype=bool) + + if "subswath_mask" in mask_type: + mask_subswath = iono_phase_obj.get_subswath_mask_array( + main_array=subswath_mask_main_image, + side_array=subswath_mask_side_image, + low_band_array=subswath_mask_image, + high_band_array=subswath_mask_image, + slant_main=main_slant, + slant_side=side_slant) + available_mask &= np.array(mask_subswath, dtype=bool) + + if "water" in mask_type: + # Extract preprocessing dictionary and open arrays + # water_mask_file is expected to have distance from the + # boundary of the water bodies. The values 0-100 represent + # the distance from the coastline and values from 101-200 + # represent the distance from inland water boundaries. + if iono_method in iono_method_sideband: + mask_image = water_mask_b_blk + else: + mask_image = water_mask_a_blk + available_mask &= np.array(mask_image, dtype=bool) + + valid_area = iono_phase_obj.get_valid_area( + main_array=main_image, + side_array=side_image, + low_band_array=sub_low_image, + diff_ms_array=diff_ms_image, + diff_low_high_band_array=diff_subband_image, + high_band_array=sub_high_image, + slant_main=main_slant, + slant_side=side_slant, + invalid_value=0) + + valid_area_coh = iono_phase_obj.get_valid_area( + main_array=main_coh_image, + side_array=side_coh_image, + diff_ms_array=diff_ms_coh_image, + low_band_array=sub_low_coh_image, + high_band_array=sub_high_coh_image, + diff_low_high_band_array=diff_coh_image, + slant_main=main_slant, + slant_side=side_slant, + invalid_value=0) + available_mask &= np.array(valid_area & valid_area_coh, dtype=bool) + + # qFSP correction should be applied before ionosphere estimation + # so compute_disp_nondisp() uses the corrected differential phase. + if iono_qfsp_correction_flag: + if iono_method == "main_diff_low_high_subband": + diff_phase = diff_subband_image + insar_mask_block = subswath_mask_image + parsed_mask = parse_insar_mask(insar_mask_block) + sec_has_anomaly = parsed_mask["sec_has_anomaly"] + ref_has_anomaly = parsed_mask["ref_has_anomaly"] + + # True where either reference or secondary has an input anomaly + mask_anomaly_array = np.asarray( + sec_has_anomaly | ref_has_anomaly, dtype=bool + ) + elif iono_method == "main_diff_ms_band": + diff_phase = diff_ms_image + insar_mask_block = subswath_mask_side_image + parsed_mask_side = parse_insar_mask(subswath_mask_side_image) + parsed_mask_main = parse_insar_mask(subswath_mask_main_image) + sec_has_anomaly_main = parsed_mask_main["sec_has_anomaly"] + ref_has_anomaly_main = parsed_mask_main["ref_has_anomaly"] + sec_has_anomaly_side = parsed_mask_side["sec_has_anomaly"] + ref_has_anomaly_side = parsed_mask_side["ref_has_anomaly"] + + sec_has_anomaly_main_deci = decimate_freq_a_array( + main_slant, + side_slant, + sec_has_anomaly_main) + ref_has_anomaly_main_deci = decimate_freq_a_array( + main_slant, + side_slant, + ref_has_anomaly_main) + mask_anomaly_array = np.asarray( + sec_has_anomaly_main_deci | ref_has_anomaly_main_deci | + sec_has_anomaly_side | ref_has_anomaly_side, dtype=bool + ) + else: + diff_phase = None + insar_mask_block = None + + if diff_phase is not None and insar_mask_block is not None: + diff_phase_original = diff_phase.copy() + + qfsp_out_dir = os.path.join( + iono_path, iono_method, pol_comb_str) + pathlib.Path(qfsp_out_dir).mkdir(parents=True, exist_ok=True) + + qfsp_original_path = os.path.join( + qfsp_out_dir, "qFSP_original") + + qfsp_corrected_path = os.path.join( + qfsp_out_dir, "qFSP_corrected") + write_array( + qfsp_original_path, + diff_phase_original, + data_type=gdal.GDT_Float32, + block_row=row_start, + data_shape=[rows_output, cols_output]) + # Use only non-anomaly pixels as background support. + # If mask_array already exists, combine it with non-anomaly pixels. + qfsp_background_mask = ( + available_mask + & (~mask_anomaly_array) + ) + + output = correct_qfsp_phase_artifact( + diff_phase, + fit_background_mask=qfsp_background_mask, + artifact_mask=mask_anomaly_array, + background_order=qfsp_background_order, + min_fraction_rows=qfsp_min_fraction_rows, + min_group_width=qfsp_min_group_width, + template_smooth_win=qfsp_template_smooth_win, + inner_shrink=qfsp_inner_shrink, + outer_feather=qfsp_outer_feather, + fill_value=np.nan, + ) + + diff_phase = output["corrected_phase"] + + if iono_method == "main_diff_low_high_subband": + diff_subband_image = diff_phase + + elif iono_method == "main_diff_ms_band": + diff_ms_image = diff_phase + # Estimate dispersive and non-dispersive phase dispersive, non_dispersive = iono_phase_obj.compute_disp_nondisp( phi_sub_low=sub_low_image, @@ -1875,86 +2066,16 @@ def run(cfg: dict, runw_hdf5: str): slant_side=side_slant) else: info_channel.log(f'{mask_type} is used for mask construction') - mask_array = np.ones([block_rows_data, cols_output], - dtype=bool) - if "coherence" in mask_type: - mask_image = iono_phase_obj.get_coherence_mask_array( - main_array=main_coh_image, - side_array=side_coh_image, - diff_ms_array=diff_ms_coh_image, - low_band_array=sub_low_coh_image, - high_band_array=sub_high_coh_image, - diff_low_high_band_array=diff_coh_image, - slant_main=main_slant, - slant_side=side_slant, - threshold=filter_coh_thresh) - mask_array = mask_array & mask_image - if "connected_components" in mask_type: - mask_image = iono_phase_obj.get_conn_component_mask_array( - main_array=main_conn_image, - side_array=side_conn_image, - diff_ms_array=diff_ms_conn_image, - low_band_array=sub_low_conn_image, - high_band_array=sub_high_conn_image, - diff_low_high_band_array=diff_subband_conn_image, - slant_main=main_slant, - slant_side=side_slant) - mask_array = mask_array & mask_image + final_iono_mask = available_mask.copy() if "median_filter" in mask_type: mask_image = iono_phase_obj.get_mask_median_filter( disp=dispersive, looks=number_looks, threshold=median_filter_threshold, - median_filter_size=median_filter_size, - ) - mask_array = mask_array & mask_image - - if "subswath_mask" in mask_type: - mask_subswath = iono_phase_obj.get_subswath_mask_array( - main_array=subswath_mask_main_image, - side_array=subswath_mask_side_image, - low_band_array=subswath_mask_image, - high_band_array=subswath_mask_image, - slant_main=main_slant, - slant_side=side_slant) - mask_array &= mask_subswath - - if "water" in mask_type: - # Extract preprocessing dictionary and open arrays - # water_mask_file is expected to have distance from the - # boundary of the water bodies. The values 0-100 represent - # the distance from the coastline and values from 101-200 - # represent the distance from inland water boundaries. - if iono_method in iono_method_sideband: - mask_image = water_mask_b_blk - else: - mask_image = water_mask_a_blk - mask_array = mask_array & mask_image - - valid_area = iono_phase_obj.get_valid_area( - main_array=main_image, - side_array=side_image, - low_band_array=sub_low_image, - diff_ms_array=diff_ms_image, - diff_low_high_band_array=diff_subband_image, - high_band_array=sub_high_image, - slant_main=main_slant, - slant_side=side_slant, - invalid_value=0) - - valid_area_coh = iono_phase_obj.get_valid_area( - main_array=main_coh_image, - side_array=side_coh_image, - diff_ms_array=diff_ms_coh_image, - low_band_array=sub_low_coh_image, - high_band_array=sub_high_coh_image, - diff_low_high_band_array=diff_coh_image, - slant_main=main_slant, - slant_side=side_slant, - invalid_value=0) - final_iono_mask = mask_array & valid_area & valid_area_coh + median_filter_size=median_filter_size) + final_iono_mask &= mask_image mask_path = os.path.join( iono_path, iono_method, pol_comb_str, 'mask_array') @@ -2127,9 +2248,19 @@ def run(cfg: dict, runw_hdf5: str): main_image = read_block_array( src_main_h5[runw_path_freq_a], block_parm, 0) - diff_subband_image = read_block_array( - src_diff_h5[runw_path_freq_a], - block_parm, 0) + if iono_qfsp_correction_flag: + qfsp_corrected_path = os.path.join( + iono_path, iono_method, pol_comb_str, + "qFSP_corrected") + + diff_subband_image = read_block_array( + qfsp_corrected_path, + block_parm, + 0) + else: + diff_subband_image = read_block_array( + src_diff_h5[runw_path_freq_a], + block_parm, 0) main_image[mask_image == 0] = 0 diff_subband_image[mask_image == 0] = 0 @@ -2189,13 +2320,23 @@ def run(cfg: dict, runw_hdf5: str): block_parm_side, 0) if iono_method == 'main_diff_ms_band': - with HDF5OptimizedReader(name=runw_diff_str, - mode='r', - libver='latest', - swmr=True) as src_diff_h5: + + if iono_qfsp_correction_flag: + qfsp_corrected_path = os.path.join( + iono_path, iono_method, pol_comb_str, "qFSP_corrected") + diff_ms_image = read_block_array( - src_diff_h5[runw_path_freq_b], - block_parm_side, 0) + qfsp_corrected_path, + block_parm_side, + 0) + else: + with HDF5OptimizedReader(name=runw_diff_str, + mode='r', + libver='latest', + swmr=True) as src_diff_h5: + diff_ms_image = read_block_array( + src_diff_h5[runw_path_freq_b], + block_parm_side, 0) if bridge_algorithm_bool: main_image = bridge_unwrapped_phase( diff --git a/share/nisar/defaults/insar.yaml b/share/nisar/defaults/insar.yaml index a6d2c5409..a82c19c8a 100644 --- a/share/nisar/defaults/insar.yaml +++ b/share/nisar/defaults/insar.yaml @@ -345,6 +345,30 @@ runconfig: # During deramping, uniform sampling is applied # when number of samples is larger than maximum pixel bridge_ramp_maximum_pixel: 1e6 + qfsp_correction: + # Apply qFSP phase-slip artifact correction. + enabled: True + # Polynomial order used to estimate the smooth background + # phase from non-anomaly pixels. + # 0: constant, 1: linear, 2: quadratic. + background_order: 2 + # Minimum fraction of rows that must contain usable artifact + # samples to estimate a qFSP artifact template. + # If the valid artifact support is smaller than this value, + # correction may be skipped or considered unreliable. + min_fraction_rows: 0.05 + # Minimum width, in pixels, of an artifact group to correct. + # Groups narrower than this value are ignored. + min_group_width: 1 + # Window size, in pixels, used to smooth the estimated qFSP + # artifact template. Use an odd integer when possible. + template_smooth_win: 3 + # Number of pixels to shrink from the inner artifact region + # before estimating/applying the correction. + inner_shrink: 2 + # Number of pixels used to feather the correction outside + # the artifact boundary. + outer_feather: 10 troposphere_delay: # Boolean flag to activate (True) troposphere delay computation diff --git a/share/nisar/schemas/insar.yaml b/share/nisar/schemas/insar.yaml index 8006d324a..cbe743435 100644 --- a/share/nisar/schemas/insar.yaml +++ b/share/nisar/schemas/insar.yaml @@ -773,6 +773,8 @@ ionosphere_phase_correction_options: # Unwrapping correction bool. If False, only filtering # is applied unwrap_correction: include('unwrap_correction_options', required=False) + # qFSP phase-slip artifact correction options + qfsp_correction: include('qfsp_correction_options', required=False) dispersive_filter_options: # Boolean flag to enable/disable ionosphere phase correction @@ -844,6 +846,22 @@ dispersive_filter_options: # when number of samples is larger than maximum pixel bridge_ramp_maximum_pixel: num(min=0, required=False) +qfsp_correction_options: + # Boolean flag to enable/disable qFSP phase-slip artifact correction + enabled: bool(required=False) + # Polynomial order used to estimate the background phase + background_order: int(min=0, required=False) + # Minimum fraction of valid rows required for estimating the qFSP artifact + min_fraction_rows: num(min=0.0, max=1.0, required=False) + # Minimum width, in pixels, of a qFSP artifact group + min_group_width: int(min=1, required=False) + # Smoothing window size used for qFSP template estimation + template_smooth_win: int(min=1, required=False) + # Number of pixels to shrink from the inner artifact region + inner_shrink: int(min=0, required=False) + # Number of pixels used for feathering around the artifact boundary + outer_feather: int(min=0, required=False) + split_range_spectrum_options: # Number of lines per block for block processing) lines_per_block: int(min=1, required=False)