Doppler shift (#34)

* range delay caused by Doppler shift

* update docs

* address multiple styling comments

* move the new functions next to bistatic_dely

* address codacy, remove space

* consolidate azimuth carrier methods
make carrier method class member

* make azimuth carrier a class

* docstring nits

* fix syntax error; add missing defs

* fix undefined orbit

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

Author: hfattahi

src/s1reader/s1_burst_slc.py

@@ -10,57 +10,6 @@
 
 
 # Other functionalities
-def compute_az_carrier(burst, orbit, offset, position):
-    '''
-    Estimate azimuth carrier and store in numpy arrary
-
-    Parameters
-    ----------
-    burst: Sentinel1BurstSlc
-       Sentinel1 burst object
-    orbit: isce3.core.Orbit
-       Sentinel1 orbit ephemerides
-    offset: float
-       Offset between reference and secondary burst
-    position: tuple
-       Tuple of locations along y and x directions
-
-    Returns
-    -------
-    carr: np.ndarray
-       Azimuth carrier
-    '''
-
-    # Get burst sensing mid relative to orbit reference epoch
-    fmt = "%Y-%m-%dT%H:%M:%S.%f"
-    orbit_ref_epoch = datetime.datetime.strptime(orbit.reference_epoch.__str__()[:-3], fmt)
-
-    t_mid = burst.sensing_mid - orbit_ref_epoch
-    _, v = orbit.interpolate(t_mid.total_seconds())
-    vs = np.linalg.norm(v)
-    ks = 2 * vs * burst.azimuth_steer_rate / burst.wavelength
-
-    y, x = position
-
-    n_lines, _ = burst.shape
-    eta = (y - (n_lines // 2) + offset) * burst.azimuth_time_interval
-    rng = burst.starting_range + x * burst.range_pixel_spacing
-
-    f_etac = np.array(
-        burst.doppler.poly1d.eval(rng.flatten().tolist())).reshape(rng.shape)
-    ka = np.array(
-        burst.azimuth_fm_rate.eval(rng.flatten().tolist())).reshape(rng.shape)
-
-    eta_ref = (burst.doppler.poly1d.eval(
-        burst.starting_range) / burst.azimuth_fm_rate.eval(
-        burst.starting_range)) - (f_etac / ka)
-    kt = ks / (1.0 - ks / ka)
-
-    carr = np.pi * kt * ((eta - eta_ref) ** 2)
-
-    return carr
-
-
 def polyfit(xin, yin, zin, azimuth_order, range_order,
             sig=None, snr=None, cond=1.0e-12,
             max_order=True):
@@ -155,6 +104,19 @@ class represents a polynomial function of range
     poly = isce3.core.Poly2d(coeffs, xmin, ymin, xnorm, ynorm)
     return poly
 
+@dataclass
+class AzimuthCarrierComponents:
+    kt: np.ndarray
+    eta: float
+    eta_ref: float
+
+    @property
+    def antenna_steering_doppler(self):
+        return self.kt * (self.eta - self.eta_ref)
+
+    @property
+    def carrier(self):
+        return np.pi * self.kt * ((self.eta - self.eta_ref) ** 2)
 
 @dataclass(frozen=True)
 class Doppler:
@@ -198,6 +160,7 @@ class Sentinel1BurstSlc:
     range_window_coefficient: float
     rank: int # The number of PRI between transmitted pulse and return echo.
     prf_raw_data: float  # Pulse repetition frequency (PRF) of the raw data [Hz]
+    range_chirp_rate: float # Range chirp rate [Hz]
 
     def as_isce3_radargrid(self):
         '''Init and return isce3.product.RadarGridParameters.
@@ -339,14 +302,14 @@ class represents a polynomial function of range
         x_mesh, y_mesh = np.meshgrid(x, y)
 
         # Estimate azimuth carrier
-        az_carrier = compute_az_carrier(self, self.orbit,
+        az_carr_comp = self.az_carrier_components(
                                         offset=offset,
                                         position=(y_mesh, x_mesh))
 
         # Fit azimuth carrier polynomial with x/y or range/azimuth
         if index_as_coord:
             az_carrier_poly = polyfit(x_mesh.flatten()+1, y_mesh.flatten()+1,
-                                      az_carrier.flatten(), az_order,
+                                      az_carr_comp.carrier.flatten(), az_order,
                                       rg_order)
         else:
             # Convert x/y to range/azimuth
@@ -356,7 +319,7 @@ class represents a polynomial function of range
 
             # Estimate azimuth carrier polynomials
             az_carrier_poly = polyfit(rg_mesh.flatten(), az_mesh.flatten(),
-                                  az_carrier.flatten(), az_order,
+                                  az_carr_comp.carrier.flatten(), az_order,
                                   rg_order)
 
         return az_carrier_poly
@@ -478,6 +441,126 @@ def bistatic_delay(self, xstep=1, ystep=1):
 
         return isce3.core.LUT2d(x, y, bistatic_correction)
 
+    def geometrical_and_steering_doppler(self, xstep=500, ystep=50):
+        """
+        Compute total Doppler which is the sum of two components:
+        (1) the geometrical Doppler induced by the relative movement
+        of the sensor and target
+        (2) the TOPS specicifc Doppler caused by the electric steering
+        of the beam along the azimuth direction resulting in Doppler varying
+        with azimuth time.
+        Parameters
+        ----------
+        xstep: int
+            Spacing along x direction [pixels]
+        ystep: int
+            Spacing along y direction [pixels]
+
+        Returns
+        -------
+        x : int
+           The index of samples in range direction as an 1D array
+        y : int
+           The index of samples in azimuth direction as an 1D array
+        total_doppler : float
+           Total Doppler which is the sum of the geometrical Doppler and
+           beam steering induced Doppler [Hz] as a 2D array
+        """
+
+        x = np.arange(0, self.width, xstep, dtype=int)
+        y = np.arange(0, self.length, ystep, dtype=int)
+        x_mesh, y_mesh = np.meshgrid(x, y)
+        az_carr_comp = self.az_carrier_components(
+                                        offset=0.0,
+                                        position=(y_mesh, x_mesh))
+
+        slant_range = self.starting_range + x * self.range_pixel_spacing
+        geometrical_doppler = self.doppler.poly1d.eval(slant_range)
+
+        total_doppler = az_carr_comp.antenna_steering_doppler + geometrical_doppler
+
+        return x, y, total_doppler
+
+    def doppler_induced_range_shift(self, xstep=500, ystep=50):
+        """
+        Computes the range delay caused by the Doppler shift as described
+        by Gisinger et al 2021
+
+        Parameters
+        ----------
+        xstep: int
+            Spacing along x direction [pixels]
+        ystep: int
+            Spacing along y direction [pixels]
+
+        Returns
+        -------
+        isce3.core.LUT2d:
+           LUT2D object of range delay correction [seconds] as a function
+           of the x and y indices.
+
+        """
+
+        x, y, doppler_shift = self.geometrical_and_steering_doppler(
+                                                    xstep=xstep, ystep=ystep)
+        tau_corr = doppler_shift / self.range_chirp_rate
+
+        return isce3.core.LUT2d(x, y, tau_corr)
+
+    def az_carrier_components(self, offset, position):
+        '''
+        Estimate azimuth carrier and store in numpy arrary. Also return
+        contributing components.
+
+        Parameters
+        ----------
+        offset: float
+           Offset between reference and secondary burst
+        position: tuple
+           Tuple of locations along y and x directions
+
+        Returns
+        -------
+        eta: float
+            zero-Doppler azimuth time centered in the middle of the burst
+        eta_ref: float
+            refernce time
+        kt: np.ndarray
+            Doppler centroid rate in the focused TOPS SLC data [Hz/s]
+        carr: np.ndarray
+           Azimuth carrier
+
+        Reference
+        ---------
+        https://sentinels.copernicus.eu/documents/247904/0/Sentinel-1-TOPS-SLC_Deramping/b041f20f-e820-46b7-a3ed-af36b8eb7fa0
+        '''
+        # Get self.sensing mid relative to orbit reference epoch
+        fmt = "%Y-%m-%dT%H:%M:%S.%f"
+        orbit_ref_epoch = datetime.datetime.strptime(self.orbit.reference_epoch.__str__()[:-3], fmt)
+
+        t_mid = self.sensing_mid - orbit_ref_epoch
+        _, v = self.orbit.interpolate(t_mid.total_seconds())
+        vs = np.linalg.norm(v)
+        ks = 2 * vs * self.azimuth_steer_rate / self.wavelength
+
+        y, x = position
+
+        n_lines, _ = self.shape
+        eta = (y - (n_lines // 2) + offset) * self.azimuth_time_interval
+        rng = self.starting_range + x * self.range_pixel_spacing
+
+        f_etac = np.array(
+            self.doppler.poly1d.eval(rng.flatten().tolist())).reshape(rng.shape)
+        ka = np.array(
+            self.azimuth_fm_rate.eval(rng.flatten().tolist())).reshape(rng.shape)
+
+        eta_ref = (self.doppler.poly1d.eval(
+            self.starting_range) / self.azimuth_fm_rate.eval(
+            self.starting_range)) - (f_etac / ka)
+        kt = ks / (1.0 - ks / ka)
+
+
+        return AzimuthCarrierComponents(kt, eta, eta_ref)
 
     @property
     def sensing_mid(self):

src/s1reader/s1_reader.py

@@ -299,6 +299,7 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str,
         downlink_element =  tree.find('generalAnnotation/downlinkInformationList/downlinkInformation')
         prf_raw_data = float(downlink_element.find('prf').text)
         rank = int(downlink_element.find('downlinkValues/rank').text)
+        range_chirp_ramp_rate = float(downlink_element.find('downlinkValues/txPulseRampRate').text)
 
         n_lines = int(tree.find('swathTiming/linesPerBurst').text)
         n_samples = int(tree.find('swathTiming/samplesPerBurst').text)
@@ -399,8 +400,7 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str,
                                       tiff_path, i, first_valid_sample,
                                       last_sample, first_valid_line, last_line,
                                       range_window_type, range_window_coeff,
-                                      rank, prf_raw_data)
-
+                                      rank, prf_raw_data, range_chirp_ramp_rate)
     return bursts
 
 def _is_zip_annotation_xml(path: str, id_str: str) -> bool: