Merge pull request #142 from ACCESS-Cloud-Based-InSAR/iono_adding_bur…

…st_jumps

Iono adding burst jumps

CHANGELOG.md

@@ -20,6 +20,7 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 ### Added
 * localize_data within __main__.py added option to use/not use water mask for ionosphere processing
 * Added option to estimate burst phase jumps in ionosphere computation
+* Added additional attributes for ionosphere computation into GUNW ionosphere layer metadata: processing_steps, water_mask, mask_connected_component_zero (flag) , do_phase_bridging (flag), swath_mode (flag), swath_ramp_removal (flag), swath_mode_description, multilook_az_rg1, multilook_az_rg2, iono_height
 * Added packing of additional attributes for ionosphere computation into GUNW
 * fsspec is now required in environment due to burst processing.
 

isce2_topsapp/__main__.py

@@ -200,28 +200,23 @@ def gunw_slc():
     )
 
     # Run ionospheric correction
-    # MG: correct burst jumps, e.g. needed for Arabian
-    #   processing. TODO: We need a trigger function for
-    #   this option (it adds 10min to iono for frame processing)
-    #   example: look at filt_toposphase.unw or .flat
-    #   and analyze if there are any burst jumps,
-    #   if yes, set this option True, or we could run it always
-    #   and store both iono-long wavelength and burst jumps(az_shifts)
     if args.estimate_ionosphere_delay:
-        iono_attributes = iono_processing(
+        iono_attr = iono_processing(
             mask_filename=loc_data["water_mask"],
-            correct_burst_jumps=False,
+            correct_burst_ramps=True,
         )
 
     ref_properties = loc_data["reference_properties"]
     sec_properties = loc_data["secondary_properties"]
     extent = loc_data["extent"]
 
     additional_2d_layers = []
-    aditional_attributes = []
+    additional_attributes = []
     if args.estimate_ionosphere_delay:
         additional_2d_layers.append("ionosphere")
-        aditional_attributes.append(iono_attributes)
+        additional_attributes.append(iono_attr['ionosphere'])
+        additional_2d_layers.append("ionosphereBurstRamps")
+        additional_attributes.append(iono_attr['ionosphereBurstRamps'])
 
     additional_2d_layers = additional_2d_layers or None
     nc_path = package_gunw_product(
@@ -230,7 +225,7 @@ def gunw_slc():
         secondary_properties=sec_properties,
         extent=extent,
         additional_2d_layers=additional_2d_layers,
-        additional_attributes=aditional_attributes
+        additional_attributes=additional_attributes
     )
 
     if args.compute_solid_earth_tide:

isce2_topsapp/iono_proc.py

@@ -37,7 +37,7 @@ def iono_processing(
     *,
     topsapp_xml_filename: str = 'topsApp.xml',
     mask_filename: str = '',
-    correct_burst_jumps: bool = False,
+    correct_burst_ramps: bool = True,
         ) -> None:
 
     '''
@@ -119,8 +119,10 @@ def iono_processing(
                   coh_threshold=coh_threshold, sigma_rule=sigma_rule)
 
     # Step when using azimuth shift correction
-    # between bursts
-    if correct_burst_jumps:
+    # between bursts, misregistration due to high iono content
+    # based on long wavelength ionosphere delay estimation
+    # NOTE: wrong long-wavelength iono estimates can effect this steps
+    if correct_burst_ramps:
         # ionosphere shift
         runIon.ionosphere_shift(topsapp, ionParam)
 
@@ -131,10 +133,16 @@ def iono_processing(
         # esd
         runIon.esd(topsapp, ionParam)
 
-        # Create merged/topophase.ion file
+        # Create merged/topophase.ion.az_shift file
         # using ion/ion_burst/ files
         merge_multilook_bursts(topsapp, input_dir='ion/ion_burst',
-                               output_filename='topophase.ion')
+                               output_filename='topophase.ion.az_shift')
+
+        GEOCODE_LIST_ION.append('merged/topophase.ion.az_shift')
+
+        # Iono long wavelength
+        merge_bursts(input_file="ion/ion_cal/filt.ion",
+                     output_filename="topophase.ion")
 
     else:
         # Create merged/topophase.ion file
@@ -146,28 +154,17 @@ def iono_processing(
                        topsapp.geocode_bbox)
 
     # Create attribute iono processing dict
-    if correct_burst_jumps:
-        steps = ionParam.allSteps
-    else:
-        steps = ionParam.allSteps[:-3]
+    steps = ionParam.allSteps
     # as we skip grd2ion projection to ionospheric thin shell
     del steps[2]
-    steps.append('geocode')
 
     # attributes dict items must be in the following
     # format: str, Number, ndarray, number, list, tuple
     iono_dict = dict(
-        processing_steps=steps,
-        mask=str(mask_filename),
+        processing_steps=steps[0:3] + ['geocode'],
+        water_mask=str(mask_filename),
         mask_connected_component_zero=str(conncomp_flag),
         do_phase_bridging=str(True),
-        estimate_burst_jumps=str(ionParam.considerBurstProperties),
-        burst_jumps_description=('calculation of burst jumps'
-                                 ' (scalloping effect) due'
-                                 ' to misaligment in coregistration'
-                                 ' caused be large ionosphere'
-                                 ' content/delay, typically low or'
-                                 ' high latitudes'),
         swath_mode=str(False) if ionParam.calIonWithMerged else str(True),
         swath_ramp_removal=str(np.bool_(ionParam.rampRemovel)),
         swath_mode_description=('raw_ion ionosphere calculation'
@@ -184,7 +181,14 @@ def iono_processing(
         iono_height=ionParam.ionHeight
         )
 
-    return iono_dict
+    burst_ramps_dict = dict(
+        processing_steps=steps + ['geocode'],
+    )
+
+    attr = dict(ionosphere=iono_dict,
+                ionosphereBurstRamps=burst_ramps_dict)
+
+    return attr
 
 
 def mask_iono_ifg_bursts(tops_dir: Path,

isce2_topsapp/packaging.py

@@ -11,13 +11,13 @@
 
 import isce2_topsapp
 from isce2_topsapp.packaging_utils.additional_layers import add_2d_layer
-from isce2_topsapp.packaging_utils.ionosphere import format_ionosphere_for_gunw
+from isce2_topsapp.packaging_utils.ionosphere import format_iono_burst_ramps, format_ionosphere_for_gunw
 from isce2_topsapp.templates import read_netcdf_packaging_template
 
 DATASET_VERSION = '3.0.0'
 
 
-PERMISSIBLE_2D_LAYERS = ['ionosphere']
+PERMISSIBLE_2D_LAYERS = ['ionosphere', 'ionosphereBurstRamps']
 
 
 """Warning: the packaging scripts were written as command line scripts and
@@ -240,6 +240,9 @@ def package_additional_layers_into_gunw(gunw_path: Path,
     if 'ionosphere' in additional_2d_layers:
         # current working directory is ISCE directory
         _ = format_ionosphere_for_gunw(isce_data_directory, gunw_path)
+    if 'ionosphereBurstRamps' in additional_2d_layers:
+        # current working directory is ISCE directory
+        _ = format_iono_burst_ramps(isce_data_directory, gunw_path)
 
     # Assumes ionosphere raster is written to specific path
     additional_dataset = zip(additional_2d_layers, additional_attributes)

isce2_topsapp/packaging_utils/additional_layers.json

@@ -6,6 +6,15 @@
                                        "long_name": "ionospherePhaseCorrection",
                                        "description": "Estimated ionosphere phase correction"
                                        }
-                             }
+                             },
+ "ionosphereBurstRamps": {"dst_group": "/science/grids/corrections/derived/ionosphereBurstRamps",
+                    "dst_variable": "ionosphereBurstRamps",
+                    "input_relative_path": "merged/ionosphereBurstRamps_for_gunw.geo",
+                    "attrs": {
+                                "standard_name": "ionosphereBurstRamps",
+                                "long_name": "ionosphereBurstRamps",
+                                "description": "Estimated burst ramps correction due to ionosphere"
+                                }
+              } 
 
 }

isce2_topsapp/packaging_utils/additional_layers.py

@@ -18,6 +18,7 @@ def add_2d_layer(layer_name: str, gunw_netcdf_path: Path, additional_attrs: dict
     """
 
     layer_data = ADDITIONAL_LAYERS[layer_name]
+
     dst_group = layer_data['dst_group']
     dst_variable = layer_data['dst_variable']
     if additional_attrs:

isce2_topsapp/packaging_utils/ionosphere.py

@@ -45,3 +45,28 @@ def format_ionosphere_for_gunw(isce_directory: Path,
         ds.write(X_ion_low_res, 1)
 
     return out_path
+
+
+def format_iono_burst_ramps(isce_directory: Path,
+                            gunw_netcdf_path: Path) -> Path:
+    # open long-wavelength ionosphere layer
+    with rasterio.open(isce_directory / "merged/topophase.ion.geo") as ds:
+        X_ion = ds.read(1)
+        p_ion = ds.profile
+
+    X_ion[X_ion == 0] = np.nan
+    p_ion['nodata'] = np.nan
+
+    # open burstRamps layer
+    with rasterio.open(isce_directory / "merged/topophase.ion.az_shift.geo") as ds:
+        X_ramps = ds.read(1)
+    X_ramps[X_ramps == 0] = np.nan
+
+    # Get burst ramps without long-wavength ionospheric delay
+    X_ramps -= X_ion
+
+    out_path = isce_directory / 'merged/ionosphereBurstRamps_for_gunw.geo'
+    with rasterio.open(out_path, 'w', **p_ion) as ds:
+        ds.write(X_ramps, 1)
+
+    return out_path