diff --git a/python/packages/isce3/atmosphere/main_band_estimation.py b/python/packages/isce3/atmosphere/main_band_estimation.py index 5d6e5cc40..b8c58967d 100644 --- a/python/packages/isce3/atmosphere/main_band_estimation.py +++ b/python/packages/isce3/atmosphere/main_band_estimation.py @@ -421,6 +421,7 @@ def estimate_iono_std( self, main_coh=None, side_coh=None, + diff_ms_coh=None, low_band_coh=None, high_band_coh=None, diff_low_high_coh=None, @@ -437,6 +438,8 @@ def estimate_iono_std( coherence of main-band interferogram side_coh : numpy.ndarray coherence of side-band interferogram + diff_ms_coh : numpy.ndarray + coherence of difference (main-side) interferogram low_band_coh : numpy.ndarray coherence of low subband interferogram high_band_coh : numpy.ndarray @@ -458,9 +461,22 @@ def estimate_iono_std( sig_nondisp : numpy.ndarray phase standard deviation of non-dispersive """ + error_channel = journal.error( + 'MainBandIonosphereEstimation.estimate_iono_std') + if side_coh is None and diff_ms_coh is None: + err_str = "estimate_iono_std requires at least one of "\ + "`diff_ms_coh` or `side_coh` to be provided." + error_channel.log(err_str) + raise ValueError(err_str) + if main_coh is None: + err_str = "estimate_iono_std requires `main_coh` to be provided." + error_channel.log(err_str) + raise ValueError(err_str) + side_used_coh = diff_ms_coh if diff_ms_coh is not None else side_coh + # resample coherences array of frequency A to # frequency B grid - if (side_coh is not None) and (resample_flag): + if resample_flag: if slant_main is None: slant_main = self.slant_main if slant_side is None: @@ -471,18 +487,18 @@ def estimate_iono_std( slant_side, main_coh) + sig_phi_main = np.divide( + np.sqrt(1 - main_coh**2), + main_coh * np.sqrt(2 * number_looks), + out=np.zeros_like(main_coh), + where=main_coh != 0) + # estimate sigma from main- and side- band coherences - if (main_coh is not None) & (side_coh is not None): - sig_phi_main = np.divide( - np.sqrt(1 - main_coh**2), - main_coh / np.sqrt(2 * number_looks), - out=np.zeros_like(main_coh), - where=main_coh != 0) - sig_phi_side = np.divide( - np.sqrt(1 - side_coh**2), - side_coh / np.sqrt(2 * number_looks), - out=np.zeros_like(side_coh), - where=side_coh != 0) + sig_phi_side = np.divide( + np.sqrt(1 - side_used_coh**2), + side_used_coh * np.sqrt(2 * number_looks), + out=np.zeros_like(side_used_coh), + where=side_used_coh != 0) sig_phi_iono, sig_nondisp = \ self.estimate_sigma( @@ -491,37 +507,6 @@ def estimate_iono_std( return sig_phi_iono, sig_nondisp - def estimate_sigma_main_diff( - self, - sig_phi0, - sig_phi1): - """Estimate sigma from coherence for main_diff_ms method - - Parameters - ---------- - sig_phi0 : numpy.ndarray - phase standard deviation of main-band interferogram - sig_phi1 : numpy.ndarray - phase standard deviation of side-band interferogram - - Returns - ------- - sig_iono : numpy.ndarray - 2D array of phase standard deviation of dispersive - sig_nondisp : numpy.ndarray - 2D array of phase standard deviation of non-dispersive - """ - - a = self.f1 / (self.f1 + self.f0) - b = (self.f0 * self.f1) / (self.f0**2 - self.f1**2) - sig_phi01 = np.sqrt(sig_phi0**2 + sig_phi1**2) - sig_iono = np.sqrt(a**2 * sig_phi0**2 + b**2 * sig_phi01**2) - - c = self.f0**2 / (self.f0**2 - self.f1**2) - d = self.f0 * self.f1 / (self.f0**2 - self.f1**2) - sig_nondisp = np.sqrt(c**2 * sig_phi0**2 + d**2 * sig_phi01**2) - - return sig_iono, sig_nondisp def estimate_sigma_main_side( self, @@ -669,7 +654,7 @@ def __init__(self, high_center_freq, method) self.estimate_iono = estimate_iono_main_diff - self.estimate_sigma = self.estimate_sigma_main_diff + self.estimate_sigma = self.estimate_sigma_main_side self.compute_unwrap_err = compute_unwrapp_error_main_diff_ms_band diff --git a/python/packages/isce3/atmosphere/split_band_estimation.py b/python/packages/isce3/atmosphere/split_band_estimation.py index f93ccf53e..2f4cdee6b 100644 --- a/python/packages/isce3/atmosphere/split_band_estimation.py +++ b/python/packages/isce3/atmosphere/split_band_estimation.py @@ -461,6 +461,7 @@ def estimate_iono_std( self, main_coh=None, side_coh=None, + diff_ms_coh=None, low_band_coh=None, high_band_coh=None, diff_low_high_coh=None, @@ -477,6 +478,8 @@ def estimate_iono_std( coherence of main-band interferogram side_coh : numpy.ndarray coherence of side-band interferogram + diff_ms_coh : numpy.ndarray + coherence of difference (main-side) interferogram low_band_coh : numpy.ndarray coherence of low subband interferogram high_band_coh : numpy.ndarray @@ -817,19 +820,15 @@ def compute_unwrapp_error_main_diff_low_high( """ freq_diff = freq_high - freq_low - freq_sum = freq_high + freq_low - freq_multi = freq_high * freq_low - x_coeff = freq_multi / f0 / freq_sum - z_coeff = - freq_multi / 2 / f0 / freq_diff + diff_phase = diff_low_high_runw - \ + (freq_diff / f0 * nondisp_array + f0 * + (1/freq_high - 1/freq_low) * disp_array) - diff_unw_coeff = np.round((2 * nondisp_array / z_coeff - - disp_array / z_coeff * x_coeff - - diff_low_high_runw) / 2 / np.pi) - com_unw_coeff = np.round(( - (nondisp_array + disp_array) / (2 * x_coeff + 1) - - main_runw / 2) / np.pi) + comm_phase = main_runw - (nondisp_array + disp_array) + com_unw_coeff = np.rint(comm_phase / (2*np.pi)).astype(np.int32) + diff_unw_coeff = np.rint(diff_phase / (2*np.pi)).astype(np.int32) return com_unw_coeff, diff_unw_coeff diff --git a/python/packages/nisar/workflows/ionosphere.py b/python/packages/nisar/workflows/ionosphere.py index f9d1d5eb4..62e5e1dbe 100644 --- a/python/packages/nisar/workflows/ionosphere.py +++ b/python/packages/nisar/workflows/ionosphere.py @@ -1262,6 +1262,7 @@ def run(cfg: dict, runw_hdf5: str): main_conn_image = None side_conn_image = None diff_ms_conn_image = None + diff_subband_conn_image = None # ionosphere method using sub-bands uses subswath mask in freq A # ionosphere method using side-band uses subswath masks in freq A and B @@ -1416,7 +1417,7 @@ def run(cfg: dict, runw_hdf5: str): np.s_[ row_start:row_start + block_rows_data, :]) - src_high_h5[rcom_path_freq_a].read_direct( + src_diff_h5[rcom_path_freq_a].read_direct( diff_subband_conn_image, np.s_[ row_start:row_start + block_rows_data, @@ -1446,9 +1447,10 @@ def run(cfg: dict, runw_hdf5: str): side_conn_image = np.empty( [block_rows_data, cols_side], dtype=float) - diff_ms_conn_image = np.empty( - [block_rows_data, cols_side], - dtype=float) + if iono_method == 'main_diff_ms_band': + diff_ms_conn_image = np.empty( + [block_rows_data, cols_side], + dtype=float) if "subswath_mask" in mask_type: subswath_mask_main_image = np.empty( @@ -1525,7 +1527,7 @@ def run(cfg: dict, runw_hdf5: str): :] ) if "connected_components" in mask_type: - src_side_h5[rcom_path_freq_b].read_direct( + src_diff_h5[rcom_path_freq_b].read_direct( diff_ms_conn_image, np.s_[row_start:row_start + block_rows_data, :]) @@ -1537,13 +1539,23 @@ def run(cfg: dict, runw_hdf5: str): erosion_size=bridge_erosion_size, ramp_type=bridge_deramp_type, deramp_max_num_sample=bridge_ramp_maximum_pixel) - side_image = bridge_unwrapped_phase( - side_image, - radius=bridge_radius, - min_num_pixel=bridge_minimum_samples, - erosion_size=bridge_erosion_size, - ramp_type=bridge_deramp_type, - deramp_max_num_sample=bridge_ramp_maximum_pixel) + + if iono_method == 'main_side_band': + side_image = bridge_unwrapped_phase( + side_image, + radius=bridge_radius, + min_num_pixel=bridge_minimum_samples, + erosion_size=bridge_erosion_size, + ramp_type=bridge_deramp_type, + deramp_max_num_sample=bridge_ramp_maximum_pixel) + elif iono_method == 'main_diff_ms_band': + diff_ms_image = bridge_unwrapped_phase( + diff_ms_image, + radius=bridge_radius, + min_num_pixel=bridge_minimum_samples, + erosion_size=bridge_erosion_size, + ramp_type=bridge_deramp_type, + deramp_max_num_sample=bridge_ramp_maximum_pixel) if unwrap_correction_bool: main_image = unwrapping_correction_with_filter( @@ -1554,14 +1566,25 @@ def run(cfg: dict, runw_hdf5: str): sig_kernel_y=kernel_sigma_azimuth, iterations=filter_iterations, filter_method='convolution') - side_image = unwrapping_correction_with_filter( - side_image, - kernel_width=kernel_range_size, - kernel_length=kernel_azimuth_size, - sig_kernel_x=kernel_sigma_range, - sig_kernel_y=kernel_sigma_azimuth, - iterations=filter_iterations, - filter_method='convolution') + + if iono_method == 'main_side_band': + side_image = unwrapping_correction_with_filter( + side_image, + kernel_width=kernel_range_size, + kernel_length=kernel_azimuth_size, + sig_kernel_x=kernel_sigma_range, + sig_kernel_y=kernel_sigma_azimuth, + iterations=filter_iterations, + filter_method='convolution') + elif iono_method == 'main_diff_ms_band': + diff_ms_image = unwrapping_correction_with_filter( + diff_ms_image, + kernel_width=kernel_range_size, + kernel_length=kernel_azimuth_size, + sig_kernel_x=kernel_sigma_range, + sig_kernel_y=kernel_sigma_azimuth, + iterations=filter_iterations, + filter_method='convolution') # Estimate dispersive and non-dispersive phase dispersive, non_dispersive = iono_phase_obj.compute_disp_nondisp( @@ -1611,6 +1634,7 @@ def run(cfg: dict, runw_hdf5: str): iono_std, nondisp_std = iono_phase_obj.estimate_iono_std( main_coh=main_coh_image, side_coh=side_coh_image, + diff_ms_coh=diff_ms_coh_image, low_band_coh=sub_low_coh_image, high_band_coh=sub_high_coh_image, diff_low_high_coh=diff_coh_image, @@ -1687,6 +1711,7 @@ def run(cfg: dict, runw_hdf5: str): 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 @@ -1701,14 +1726,14 @@ def run(cfg: dict, runw_hdf5: str): mask_array = mask_array & mask_image if "subswath_mask" in mask_type: - mask_array = iono_phase_obj.get_subswath_mask_array( + 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_array & mask_image + mask_array &= mask_subswath if "water" in mask_type: # Extract preprocessing dictionary and open arrays @@ -1844,7 +1869,7 @@ def run(cfg: dict, runw_hdf5: str): zip(block_params, block_params_side)): block_rows_data = block_parm.read_length - row_start = block * blocksize + row_start = block_ind * blocksize if iono_method in iono_method_sideband: block_parm_iono = block_parm_side @@ -1977,25 +2002,43 @@ def run(cfg: dict, runw_hdf5: str): erosion_size=bridge_erosion_size, ramp_type=bridge_deramp_type, deramp_max_num_sample=bridge_ramp_maximum_pixel) - side_image = bridge_unwrapped_phase( - side_image, - radius=bridge_radius, - min_num_pixel=bridge_minimum_samples, - erosion_size=bridge_erosion_size, - ramp_type=bridge_deramp_type, - deramp_max_num_sample=bridge_ramp_maximum_pixel) + + if iono_method == 'main_side_band': + side_image = bridge_unwrapped_phase( + side_image, + radius=bridge_radius, + min_num_pixel=bridge_minimum_samples, + erosion_size=bridge_erosion_size, + ramp_type=bridge_deramp_type, + deramp_max_num_sample=bridge_ramp_maximum_pixel) + + elif iono_method == 'main_diff_ms_band': + diff_ms_image = bridge_unwrapped_phase( + diff_ms_image, + radius=bridge_radius, + min_num_pixel=bridge_minimum_samples, + erosion_size=bridge_erosion_size, + ramp_type=bridge_deramp_type, + deramp_max_num_sample=bridge_ramp_maximum_pixel) if block_ind > 0: main_image, _ = compute_phase_jump( previous_main_with_pad, main_image, half_pad_length) - side_image, _ = compute_phase_jump( - previous_side_with_pad, - side_image, - half_pad_length) + if iono_method == 'main_side_band': + side_image, _ = compute_phase_jump( + previous_side_with_pad, + side_image, + half_pad_length) + elif iono_method == 'main_diff_ms_band': + diff_ms_image, _ = compute_phase_jump( + previous_diff_ms_with_pad, + diff_ms_image, + half_pad_length) previous_main_with_pad = main_image previous_side_with_pad = side_image + previous_diff_ms_with_pad = diff_ms_image # Estimating phase unwrapping errors com_unw_err, diff_unw_err = \ diff --git a/tests/python/packages/isce3/atmosphere/ionosphere.py b/tests/python/packages/isce3/atmosphere/ionosphere.py index 842c45d1b..0bffeebb6 100644 --- a/tests/python/packages/isce3/atmosphere/ionosphere.py +++ b/tests/python/packages/isce3/atmosphere/ionosphere.py @@ -2,35 +2,38 @@ from osgeo import gdal import isce3 +from isce3.atmosphere.ionosphere_filter import IonosphereFilter, write_array +from isce3.atmosphere.main_band_estimation import ( + compute_unwrapp_error_main_diff_ms_band, + compute_unwrapp_error_main_side_band, + estimate_iono_main_side, + estimate_iono_main_diff) +from isce3.atmosphere.split_band_estimation import ( + compute_unwrapp_error_split_main_band, + estimate_iono_low_high, + estimate_iono_main_diff_low_high) from isce3.signal.interpolate_by_range import decimate_freq_a_array -from isce3.atmosphere.ionosphere_filter import IonosphereFilter, write_array -from isce3.atmosphere.main_band_estimation import (compute_unwrapp_error_main_diff_ms_band, - compute_unwrapp_error_main_side_band, - estimate_iono_main_side, - estimate_iono_main_diff) -from isce3.atmosphere.split_band_estimation import (compute_unwrapp_error_split_main_band, - estimate_iono_low_high, - estimate_iono_main_diff_low_high) def simulate_ifgrams(f, dr, dTEC): phi_non_dispersive = 4.0*np.pi*f*dr/isce3.core.speed_of_light K = 40.31 - phi_TEC = (-4.0*np.pi*K/(isce3.core.speed_of_light*f))*dTEC + phi_TEC = (-4.0 * np.pi * K / (isce3.core.speed_of_light * f)) * dTEC phi = phi_non_dispersive + phi_TEC - return phi, phi_non_dispersive , phi_TEC + return phi, phi_non_dispersive, phi_TEC + def test_ionosphere_methods(): f0 = 1.27e9 f1 = f0 + 60.0e6 - BW = 20.0e6 # bandwidth of the main band + BW = 20.0e6 # bandwidth of the main band - f0L = f0 - BW/3.0 - f0H = f0 + BW/3.0 + f0L = f0 - BW / 3.0 + f0H = f0 + BW / 3.0 dr = np.array([[0.2]]) - dTEC = np.array([[2.0*1e16]]) + dTEC = np.array([[2.0 * 1e16]]) phi0, phi0_non, phi0_iono = simulate_ifgrams(f0, dr, dTEC) phi1, j0, j1 = simulate_ifgrams(f1, dr, dTEC) phi0L, j0, j1 = simulate_ifgrams(f0L, dr, dTEC) @@ -44,7 +47,7 @@ def test_ionosphere_methods(): freq_high=f0H, phi0_low=phi0L, phi0_high=phi0H) - + phi_iono_main_LH, phi_n_main_LH = estimate_iono_main_diff_low_high( f0=f0, freq_low=f0L, @@ -70,10 +73,11 @@ def test_ionosphere_methods(): assert difference_ref_ms_abs < 1e-5 assert difference_ref_md_abs < 1e-5 + def test_unwrap_error_methods(): f0 = 1.27e9 f1 = f0 + 60.0e6 - BW = 20.0e6 # bandwidt + BW = 20.0e6 # bandwidth f0L = f0 - BW/3.0 f0H = f0 + BW/3.0 @@ -98,7 +102,12 @@ def test_unwrap_error_methods(): diff_ref_err[50:100, 50:100] = diff_ref_err[50:100, 50:100] + 4 * np.pi # test for split_main_band - phi_iono_lh, phi_n_lh = estimate_iono_low_high(f0, f0L, f0H, phaseL, phaseH) + phi_iono_lh, phi_n_lh = estimate_iono_low_high( + f0, + f0L, + f0H, + phaseL, + phaseH) # add unwrapping errors to subbands phaseL_unwErr = phaseL.copy() @@ -110,18 +119,28 @@ def test_unwrap_error_methods(): # assume that ionosphere phase is correctly estimated through filtering com_unw_err, diff_unw_err = compute_unwrapp_error_split_main_band( - f0=f0,freq_low=f0L, freq_high=f0H, - disp_array=phi_iono_lh, nondisp_array=phi_n_lh, - low_sub_runw=phaseL_unwErr, high_sub_runw=phaseH_unwErr) - - difference_comref_ls_abs = np.sum(np.abs(com_unw_err * 2*np.pi - common_ref_err)) - difference_diffref_ls_abs = np.sum(np.abs(diff_unw_err* 2*np.pi - diff_ref_err)) + f0=f0, + freq_low=f0L, + freq_high=f0H, + disp_array=phi_iono_lh, + nondisp_array=phi_n_lh, + low_sub_runw=phaseL_unwErr, + high_sub_runw=phaseH_unwErr) + + difference_comref_ls_abs = np.sum( + np.abs(com_unw_err * 2*np.pi - common_ref_err)) + difference_diffref_ls_abs = np.sum( + np.abs(diff_unw_err * 2*np.pi - diff_ref_err)) assert difference_comref_ls_abs < 1e-5 assert difference_diffref_ls_abs < 1e-5 # test for main_side_band - phi_iono_ms, phi_n_ms = estimate_iono_main_side(f0, f1, phase0, phaseSideBand) + phi_iono_ms, phi_n_ms = estimate_iono_main_side( + f0, + f1, + phase0, + phaseSideBand) # add unwrapping errors to subbands phaseSideBand_unwErr = phaseSideBand.copy() @@ -137,8 +156,10 @@ def test_unwrap_error_methods(): disp_array=phi_iono_ms, nondisp_array=phi_n_ms, main_runw=phase0_unwErr, side_runw=phaseSideBand_unwErr) - difference_comref_ms_abs = np.sum(np.abs(com_unw_err * 2*np.pi - common_ref_err)) - difference_diffref_ms_abs = np.sum(np.abs(diff_unw_err* 2*np.pi - diff_ref_err)) + difference_comref_ms_abs = np.sum( + np.abs(com_unw_err * 2*np.pi - common_ref_err)) + difference_diffref_ms_abs = np.sum( + np.abs(diff_unw_err * 2*np.pi - diff_ref_err)) assert difference_comref_ms_abs < 1e-5 assert difference_diffref_ms_abs < 1e-5 @@ -159,12 +180,15 @@ def test_unwrap_error_methods(): disp_array=phi_iono_md, nondisp_array=phi_n_md, main_runw=phase0_unwErr, diff_ms_runw=phase_diffband_unwerr) - difference_comref_md_abs = np.sum(np.abs(com_unw_err * 2*np.pi - common_ref_err)) - difference_diffref_md_abs = np.sum(np.abs(diff_unw_err* 2*np.pi - diff_ref_err)) + difference_comref_md_abs = np.sum(np.abs( + com_unw_err * 2 * np.pi - common_ref_err)) + difference_diffref_md_abs = np.sum(np.abs( + diff_unw_err * 2 * np.pi - diff_ref_err)) assert difference_comref_md_abs < 1e-5 assert difference_diffref_md_abs < 1e-5 + def test_ionosphere_filter(): ''' Check if filtered dispersive changes depending on iteration number @@ -201,15 +225,20 @@ def test_ionosphere_filter(): maskraster[50:70, 50:70] = 0 # write data into files - write_array(simul_disp_path, - phase0, data_shape=np.shape(phase0)) - write_array(simul_disp_sig_path, - phase_sig, data_shape=np.shape(phase_sig)) - write_array(mask_path, - maskraster, data_shape=np.shape(maskraster)) - - #initialize ionosphere filter object - + write_array( + simul_disp_path, + phase0, + data_shape=np.shape(phase0)) + write_array( + simul_disp_sig_path, + phase_sig, + data_shape=np.shape(phase_sig)) + write_array( + mask_path, + maskraster, + data_shape=np.shape(maskraster)) + + # initialize ionosphere filter object iono_filter_obj = IonosphereFilter( x_kernel=kernel_range_size, y_kernel=kernel_azimuth_size, @@ -258,22 +287,27 @@ def test_ionosphere_filter(): del filt_iter_3_gdal # only compare the regions which is not affected by invalid region difference = np.sum(np.abs( - filt_iter_1[int(kernel_azimuth_size/2):30, int(kernel_range_size/2):30] - - filt_iter_3[int(kernel_azimuth_size/2):30, int(kernel_range_size/2):30])) + filt_iter_1[int(kernel_azimuth_size/2):30, + int(kernel_range_size/2):30] - + filt_iter_3[int(kernel_azimuth_size/2):30, + int(kernel_range_size/2):30])) assert difference < 1e-5 + def test_decimate_runw(): main_slant = np.arange(500, 1000, 2) side_slant = np.arange(500, 1000, 4) - main_runw = np.reshape(np.arange(500, 1000, 2) ,[1, -1]) + main_runw = np.reshape(np.arange(500, 1000, 2), [1, -1]) - decimate_test = decimate_freq_a_array(main_slant, + decimate_test = decimate_freq_a_array( + main_slant, side_slant, main_runw) difference = np.sum(np.abs(decimate_test - side_slant)) assert difference < 1e-5 + if __name__ == '__main__': test_ionosphere_methods() test_unwrap_error_methods()