Merge pull request #883 from jungkyoj/split_main_band

Split-spectrum (InSAR ionosphere phase correction)

Author: jungkyoJung

python/packages/isce3/splitspectrum/splitspectrum.py

@@ -295,8 +295,11 @@ def bandpass_spectrum(self,
             )
 
         # remove the windowing effect from the spectrum 
-        spectrum_target = fft(slc_raster, n=fft_size) / \
-                              window_target
+        spectrum_target = fft(slc_raster, n=fft_size)
+        spectrum_target = np.divide(spectrum_target,
+                                    window_target,
+                                    out=np.zeros_like(spectrum_target),
+                                    where=window_target!=0)
         # apply new bandpass window to spectrum 
         slc_bandpassed = ifft(spectrum_target 
                       * window_bandpass

python/packages/nisar/workflows/insar.py

@@ -2,10 +2,10 @@
 import time
 
 import journal
-from nisar.workflows import (crossmul, dense_offsets, geo2rdr,
+from nisar.workflows import (bandpass_insar, crossmul, dense_offsets, geo2rdr,
                              geocode_insar, h5_prep, filter_interferogram,
-                             rdr2geo, resample_slc, rubbersheet, unwrap,
-                             bandpass_insar)
+                             rdr2geo, resample_slc, rubbersheet,
+                             split_spectrum, unwrap)
 from nisar.workflows.insar_runconfig import InsarRunConfig
 from nisar.workflows.persistence import Persistence
 from nisar.workflows.yaml_argparse import YamlArgparse
@@ -64,6 +64,10 @@ def run(cfg: dict, out_paths: dict, run_steps: dict):
     if run_steps['unwrap'] and 'RUNW' in out_paths:
         unwrap.run(cfg, out_paths['RIFG'], out_paths['RUNW'])
 
+    if run_steps['ionosphere'] and \
+            cfg['processing']['ionosphere_phase_correction']['enabled']:
+        split_spectrum.run(cfg)
+
     if run_steps['geocode'] and 'GUNW' in out_paths:
         geocode_insar.run(cfg, out_paths['RUNW'], out_paths['GUNW'])
 
@@ -86,4 +90,4 @@ def run(cfg: dict, out_paths: dict, run_steps: dict):
     if persist.run:
         _, out_paths = h5_prep.get_products_and_paths(insar_runcfg.cfg)
 
-        run(insar_runcfg.cfg, out_paths, persist.run_steps)
+        run(insar_runcfg.cfg, out_paths, persist.run_steps)

python/packages/nisar/workflows/insar_runconfig.py

@@ -1,6 +1,10 @@
+from nisar.products.readers import SLC
 from nisar.workflows.geo2rdr_runconfig import Geo2rdrRunConfig
+import nisar.workflows.helpers as helpers
 import journal
 import os
+import h5py
+import numpy as np
 import warnings
 
 class InsarRunConfig(Geo2rdrRunConfig):
@@ -150,6 +154,185 @@ def yaml_check(self):
         if algorithm not in self.cfg['processing']['phase_unwrap']:
             self.cfg['processing']['phase_unwrap'][algorithm]={}
 
+        iono_cfg = self.cfg['processing']['ionosphere_phase_correction']
+
+        # If ionosphere phase correction is enabled, check defaults
+        if iono_cfg['enabled']:
+            # Extract split-spectrum dictionary
+            split_cfg = iono_cfg['split_range_spectrum']
+            iono_freq_pol = iono_cfg['list_of_frequencies']
+            iono_method = iono_cfg['spectral_diversity']
+
+            # Extract main range bandwidth from reference SLC
+            ref_slc_path = self.cfg['input_file_group']['input_file_path']
+            sec_slc_path = self.cfg['input_file_group']['secondary_file_path']
+        
+            ref_slc = SLC(hdf5file=ref_slc_path)
+            sec_slc = SLC(hdf5file=sec_slc_path)
+
+            # extract the polarizations from reference and secondary hdf5
+            with h5py.File(ref_slc_path, 'r', libver='latest', 
+                swmr=True) as ref_h5, \
+                h5py.File(sec_slc_path, 'r', libver='latest', 
+                swmr=True) as sec_h5:
+                
+                ref_pol_path = os.path.join(
+                    ref_slc.SwathPath, 'frequencyA', 'listOfPolarizations')
+                ref_pols_freqA = list(
+                    np.array(ref_h5[ref_pol_path][()], dtype=str))
+                
+                sec_pol_path = os.path.join(
+                    sec_slc.SwathPath, 'frequencyA', 'listOfPolarizations')
+                sec_pols_freqA = list(
+                    np.array(sec_h5[sec_pol_path][()], dtype=str))
+
+            rg_main_bandwidth = ref_slc.getSwathMetadata(
+                'A').processed_range_bandwidth
+
+            # get common polarzations of freqA from reference and secondary
+            common_pol_refsec_freqA = set.intersection(
+                set(ref_pols_freqA), set(sec_pols_freqA))
+
+            # If no common polarizations found between reference and secondary, 
+            # then throw errors. 
+            if not common_pol_refsec_freqA:
+                err_str = "No common polarization between frequency A rasters"
+                error_channel.log(err_str)
+                raise FileNotFoundError(err_str)
+
+            # If polarizations are given, then check if HDF5 has them. 
+            # If not, then throw error. 
+            if iono_freq_pol['A']:
+                for iono_pol in iono_freq_pol['A']:
+                    if (iono_pol not in ref_pols_freqA) or \
+                       (iono_pol not in sec_pols_freqA):
+                        err_str = f"polarzations {iono_pol} for ionosphere estimation are requested, but not found"
+                        error_channel.log(err_str)
+                        raise FileNotFoundError(err_str)
+
+            # If common polarization found, but input polarizations are not given, 
+            # then assign the common polarization for split_main_band
+            if (common_pol_refsec_freqA) and (not iono_freq_pol['A']):
+                # Co-polarizations are found, split_main_band will be used for co-pols
+                common_copol_ref_sec = [pol for pol in common_pol_refsec_freqA 
+                    if pol in ['VV', 'HH']]
+                iono_freq_pol['A'] = common_copol_ref_sec
+                    
+                # If common co-pols not found, cross-pol will be alternatively used.
+                if not common_copol_ref_sec:
+                    iono_freq_pol['A'] = common_pol_refsec_freqA
+
+                warning_str = f"{iono_freq_pol} will be used for {iono_method}"
+                warning_channel.log(warning_str)
+                self.cfg['processing'][
+                    'ionosphere_phase_correction'][
+                    'list_of_frequencies'] = iono_freq_pol
+
+            # Depending on how the user has selected "spectral_diversity" check if
+            # "low_bandwidth" and "high_bandwidth" are assigned. Otherwise, use default
+            if iono_method == 'split_main_band':
+                # If "low_bandwidth" or 'high_bandwidth" is not allocated, split the main range bandwidth
+                # into two 1/3 sub-bands.
+                if split_cfg['low_band_bandwidth'] is None:
+                    split_cfg['low_band_bandwidth'] = rg_main_bandwidth / 3.0
+                    info_str = "low band width for low sub-bands are not given;"\
+                        "It is automatically set by 1/3 of range bandwidth of frequencyA"
+                    warning_channel.log(info_str)
+
+                if split_cfg['high_band_bandwidth'] is None:
+                    split_cfg['high_band_bandwidth'] = rg_main_bandwidth / 3.0
+                    info_str = "high band width for high sub-band are not given;"\
+                        "It is automatically set by 1/3 of range bandwidth of frequencyA"
+                    warning_channel.log(info_str)
+
+                # If polarzations for frequency B are requested 
+                # for split_main_band method, then throw error
+                if iono_freq_pol['B']:
+                    err_str = f"Incorrect polarzations {iono_freq_pol['B']} for frequency B are requested. "\
+                        f"{iono_method} should not have polarizations in frequency B."
+                    error_channel.log(err_str)
+                    raise FileNotFoundError(err_str)
+
+            # methods that use side band
+            if iono_method in ['main_side_band', 'main_diff_main_side_band']:
+                # extract the polarizations from reference and secondary hdf5
+                with h5py.File(ref_slc_path, 'r', libver='latest', 
+                    swmr=True) as ref_h5, \
+                    h5py.File(sec_slc_path, 'r', libver='latest', 
+                    swmr=True) as sec_h5:
+                    
+                    ref_pol_path = os.path.join(
+                        ref_slc.SwathPath, 'frequencyB', 'listOfPolarizations')
+                    ref_pols_freqB = list(
+                        np.array(ref_h5[ref_pol_path][()], dtype=str))
+                    
+                    sec_pol_path = os.path.join(
+                        sec_slc.SwathPath, 'frequencyB', 'listOfPolarizations')
+                    sec_pols_freqB = list(
+                        np.array(sec_h5[sec_pol_path][()], dtype=str))
+
+                # find common polarizations for freq B between ref and sec HDF5
+                common_pol_refsec_freqB = set.intersection(
+                    set(ref_pols_freqB), set(sec_pols_freqB))
+
+                # when common polarzations are not found, throw error. 
+                if not common_pol_refsec_freqB:
+                    err_str = "No common polarization between frequencyB rasters"
+                    error_channel.log(err_str)
+                    raise FileNotFoundError(err_str)
+
+                common_pol_ref_freq_a_b = set.intersection(
+                    set(ref_pols_freqA), set(ref_pols_freqB))
+                if not common_pol_ref_freq_a_b:
+                    err_str = "No common polarization between frequency A and B rasters"
+                    error_channel.log(err_str)
+                    raise FileNotFoundError(err_str)
+
+                # If polarizations are given, then check if HDF5 has them. 
+                # If not, then throw error. 
+                if iono_freq_pol['B']:
+                    for iono_pol in iono_freq_pol['B']:
+                        if (iono_pol not in ref_pols_freqB) or \
+                            (iono_pol not in sec_pols_freqB):
+                            err_str = f"polarzations {iono_pol} for ionosphere"\
+                                "estimation are requested, but not found"
+                            error_channel.log(err_str)
+                            raise FileNotFoundError(err_str)
+                            
+                # Co-polarizations are found and input pol for freq B is not given
+                else:
+                    common_pol_refsec_freq_ab = set.intersection(
+                    set(common_pol_refsec_freqB), set(common_pol_ref_freq_a_b))
+
+                    # if pol of freq A is given, this pol is used for freq B. 
+                    if  iono_freq_pol['A']:
+                        common_pol_refsec_freq_ab = set.intersection(
+                        set(iono_freq_pol['A']), set(common_pol_refsec_freq_ab))
+            
+                    common_copol_ref_sec = [pol for pol in common_pol_refsec_freq_ab 
+                        if pol in ['VV', 'HH']]
+
+                    if common_copol_ref_sec:
+                        info_str = f"{common_copol_ref_sec} will be "\
+                            f"used for {iono_method}"
+                        warning_channel.log(info_str)
+                        self.cfg['processing'][
+                            'ionosphere_phase_correction'][
+                            'list_of_frequencies']['A'] = common_copol_ref_sec
+                        self.cfg['processing'][
+                            'ionosphere_phase_correction'][
+                            'list_of_frequencies']['A'] = common_copol_ref_sec
+
+                    # If common co-pols not found, cross-pol will be alternatively used.
+                    else:
+                        info_str = f"{common_pol_refsec_freq_ab} will be used for split_main_band"
+                        warning_channel.log(info_str)
+                        self.cfg['processing'][
+                            'ionosphere_phase_correction'][
+                            'list_of_frequencies']['A'] = common_pol_refsec_freq_ab
+                        self.cfg['processing'][
+                            'ionosphere_phase_correction'][
+                            'list_of_frequencies']['B'] = common_pol_refsec_freq_ab
 
         if 'interp_method' not in self.cfg['processing']['geocode']:
             self.cfg['processing']['geocode']['interp_method'] = 'BILINEAR'

python/packages/nisar/workflows/persistence.py

@@ -7,9 +7,9 @@ class Persistence():
     basic class that determines InSAR persistence
     '''
     # init InSAR steps in reverse chronological run order
-    insar_steps = ['geocode', 'unwrap', 'filter_interferogram', 'crossmul', 'fine_resample', 
-                   'rubbersheet', 'dense_offsets', 'coarse_resample', 'geo2rdr', 'rdr2geo', 
-                   'h5_prep', 'bandpass_insar']
+    insar_steps = ['geocode', 'ionosphere', 'unwrap', 'filter_interferogram', 'crossmul',
+                   'fine_resample', 'rubbersheet', 'dense_offsets', 'coarse_resample',
+                   'geo2rdr', 'rdr2geo', 'h5_prep', 'bandpass_insar']
 
     def __init__(self, restart=False):
         # bool flag that determines if insar.run is called
@@ -50,18 +50,25 @@ def read_log(self):
                 if 'successfully ran INSAR' in log_line:
                     break
 
+                # success message found
+                success_msg_found = False
                 # check for message indicating successful run of step
-                if 'Successfully ran' in log_line or 'successfully ran' in log_line:
+                if 'uccessfully ran' in log_line:
                     # iterate thru reverse chronological steps
                     for insar_step in self.insar_steps:
                         # any step not found in line will be step to run
                         if insar_step not in log_line:
                             # set step name found to True
                             self.run_steps[insar_step] = True
+                            success_msg_found = True
                         else:
                             # check if any steps need to be run
                             if any(self.run_steps.values()):
                                 self.run = True
 
-                            # assume all previous steps successfully run and stop
+                            # all previous steps successfully run and stop
                             break
+
+            # check if any steps need to be run or success msg not found
+            if not success_msg_found:
+                self.__init__(True)