Patch for "no polynomial" case for azimuth FM rate mitigation & handl…

…ing the burst without burst polygon (#134)

* Take care of the case when the az. time of the polynomial list does not cover sensing start / stop

* Let `InterpolatedUnivariateSpline` to return boundary values when the input parameter is out of the interpolator's range

* extrapolate the burst center point the border polygon when the burst polygon data is missing in the metadata

* Interface to turn on / off the correction

* ditch out burst border polygon extrapolation, and fill that with `None`

* Place empty geometry instead of `None`

* whitespace removal

* Apply suggestions from code review

Co-authored-by: Liang Yu <[email protected]>

* add description about "no polynomials for burst" case

* get rid of unnecessary module in `s1-reader.py`

* bump the version

* more clarification about the missing polynomial case

---------

Co-authored-by: Seongsu Jeong <[email protected]>
Co-authored-by: Liang Yu <[email protected]>

src/s1reader/s1_annotation.py

@@ -1180,13 +1180,15 @@ def from_polynomial_lists(cls,
          fm_rate_coeff_burst_arr,
          fm_rate_tau0_burst_vec) = cls.extract_polynomial_sequence(az_fm_rate_list,
                                                                    sensing_start,
-                                                                   sensing_end)
+                                                                   sensing_end,
+                                                                   handle_out_of_range=True)
 
         (dc_aztime_burst_vec,
          dc_coeff_burst_arr,
          dc_tau0_burst_vec) = cls.extract_polynomial_sequence(doppler_centroid_list,
                                                               sensing_start,
-                                                              sensing_end)
+                                                              sensing_end,
+                                                              handle_out_of_range=True)
 
         return cls(fm_rate_aztime_burst_vec, fm_rate_coeff_burst_arr, fm_rate_tau0_burst_vec,
                    dc_aztime_burst_vec, dc_coeff_burst_arr, dc_tau0_burst_vec)
@@ -1195,7 +1197,8 @@ def from_polynomial_lists(cls,
     @classmethod
     def extract_polynomial_sequence(cls, polynomial_list: list,
                                     datetime_start: datetime.datetime,
-                                    datetime_end: datetime.datetime):
+                                    datetime_end: datetime.datetime,
+                                    handle_out_of_range=True):
         '''
         Scan `vec_azimuth_time` end find indices of the vector
         that covers the period defined with
@@ -1246,6 +1249,20 @@ def extract_polynomial_sequence(cls, polynomial_list: list,
 
         # Scale factor to convert range (in meters) to seconds (tau)
         range_to_tau = 2.0 / speed_of_light
+
+        # Take care of the case when the az. time of the polynomial list does not cover
+        # the sensing start/stop
+        if (index_end == index_start) and handle_out_of_range:
+            #      0--1--2--3--4--5 <- az. time of polynomial list (index shown on the left)
+            #|--|                   <- sensing start / stop
+            if index_start == 0:
+                index_end += 1
+
+            # 0--1--2--3--4--5      <- az. time of polynomial list (index shown on the left)
+            #                  |--| <- sensing start / stop
+            else:
+                index_start -= 1
+
         for poly in polynomial_list[index_start:index_end+1]:
             vec_aztime_sequence.append(poly[0])
             arr_coeff_sequence.append(poly[1].coeffs)

src/s1reader/s1_burst_slc.py

@@ -911,15 +911,29 @@ def az_fm_rate_mismatch_from_llh(self,
                              for datetime_vec in self.extended_coeffs.dc_aztime_vec]
 
         # calculate splined interpolation of the coeffs. and tau_0s
+        if len(fm_rate_aztime_sec_vec) <= 1:
+            # Interpolator object cannot be created with only one set of polynomials.
+            # Such case happens when there is no polygon that falls between a burst's start / stop
+            # which was found from very few Sentinel-1 IW SLCs processed by IPF 002.36
+            # Return an empty LUT2d in that case
+            return isce3.core.LUT2d()
         interpolator_tau0_ka = InterpolatedUnivariateSpline(
                                         fm_rate_aztime_sec_vec,
                                         self.extended_coeffs.fm_rate_tau0_vec,
+                                        ext=3,
                                         k=1)
         tau0_ka_interp = interpolator_tau0_ka(vec_t)[..., np.newaxis]
 
+        if len(dc_aztime_sec_vec) <=1:
+            # Interpolator object cannot be created with only one set of polynomials.
+            # Such case happens when there is no polygon that falls between a burst's start / stop
+            # which was found from very few Sentinel-1 IW SLCs processed by IPF 002.36
+            # Return an empty LUT2d in that case
+            return isce3.core.LUT2d()
         interpolator_tau0_fdc_interp = InterpolatedUnivariateSpline(
                                         dc_aztime_sec_vec,
                                         self.extended_coeffs.dc_tau0_vec,
+                                        ext=3,
                                         k=1)
         tau0_fdc_interp = interpolator_tau0_fdc_interp(vec_t)[..., np.newaxis]
 

src/s1reader/s1_reader.py

@@ -19,7 +19,7 @@
 from nisar.workflows.stage_dem import check_dateline
 
 from s1reader import s1_annotation  # to access __file__
-from s1reader.s1_annotation import (RFI_INFO_AVAILABLE_FROM, 
+from s1reader.s1_annotation import (RFI_INFO_AVAILABLE_FROM,
                                     CalibrationAnnotation,
                                     AuxCal, BurstCalibration,
                                     BurstEAP, BurstNoise, BurstExtendedCoeffs,
@@ -221,16 +221,20 @@ def calculate_centroid(lons, lats):
 
     return shapely.geometry.Point(llh_centroid[:2])
 
-def get_burst_centers_and_boundaries(tree):
+def get_burst_centers_and_boundaries(tree, num_bursts=None):
     '''Parse grid points list and calculate burst center lat and lon
 
-    Parameters:
-    -----------
+    Parameters
+    ----------
     tree : Element
         Element containing geolocation grid points.
-
-    Returns:
-    --------
+    num_bursts: int or None
+        Expected number of bursts in the subswath.
+        To check if the # of polygon is the same as the parsed polygons.
+        When it's None, the number of burst is the same as # polygons found in this function.
+    
+    Returns
+    -------
     center_pts : list
         List of burst centroids ass shapely Points
     boundary_pts : list
@@ -252,6 +256,9 @@ def get_burst_centers_and_boundaries(tree):
         lons[i] = float(grid_pt[5].text)
 
     unique_line_indices = np.unique(lines)
+
+    if num_bursts is None:
+        num_bursts = len(unique_line_indices) - 1
     n_bursts = len(unique_line_indices) - 1
     center_pts = [[]] * n_bursts
     boundary_pts = [[]] * n_bursts
@@ -273,6 +280,14 @@ def get_burst_centers_and_boundaries(tree):
         poly = shapely.geometry.Polygon(zip(burst_lons, burst_lats))
         boundary_pts[i] = check_dateline(poly)
 
+    num_border_polygon = len(unique_line_indices) - 1
+    if num_bursts > num_border_polygon:
+        warnings.warn('Inconsistency between # bursts in subswath, and the # polygons. ')
+        num_missing_polygons = num_bursts - num_border_polygon
+
+        center_pts += [shapely.Point()] * num_missing_polygons
+        boundary_pts += [[shapely.Polygon()]] * num_missing_polygons
+
     return center_pts, boundary_pts
 
 
@@ -853,8 +868,6 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str,
 
         orbit_number = int(tree.find('adsHeader/absoluteOrbitNumber').text)
 
-        center_pts, boundary_pts = get_burst_centers_and_boundaries(tree)
-
     wavelength = isce3.core.speed_of_light / radar_freq
     starting_range = slant_range_time * isce3.core.speed_of_light / 2
     range_pxl_spacing = isce3.core.speed_of_light / (2 * range_sampling_rate)
@@ -912,6 +925,8 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str,
     n_bursts = int(burst_list_elements.attrib['count'])
     bursts = [[]] * n_bursts
 
+    center_pts, boundary_pts = get_burst_centers_and_boundaries(tree, num_bursts=n_bursts)
+
     for i, burst_list_element in enumerate(burst_list_elements):
         # Zero Doppler azimuth time of the first line of this burst
         sensing_start = as_datetime(burst_list_element.find('azimuthTime').text)

src/s1reader/version.py

@@ -4,6 +4,7 @@
 # release history
 Tag = collections.namedtuple('Tag', 'version date')
 release_history = (
+    Tag('0.2.3', '2023-09-21'),
     Tag('0.2.2', '2023-09-08'),
     Tag('0.2.1', '2023-08-23'),
     Tag('0.2.0', '2023-07-25'),