Quick guide to correct a single S1 SLC burst with the S1-ETAD product

This guide should provide the basic steps which need to be performed to correct a single Sentinel-1 SLC burst with the timings provided in the S1-ETAD product.

For detailed explanations of the product format and grid definition, please consult the “Product Format Specification” Document (ETAD-DLR-PS-0014) v1.2.

All required information for timing correction is contained in the NetCDF product. The XML annotation can be consulted for informative purposes.

Notebook setup

[1]:
%matplotlib inline
[2]:
%load_ext autoreload
%autoreload 2
[3]:
import sys
sys.path.append('../..')
[4]:
import os
import datetime

from lxml import etree

import numpy as np
from scipy.constants import c
from matplotlib import pyplot as plt

1. Select a burst of a S1-SLC product and extract relevant parameters

The following parameters have been extracted by the S1 product:

  • S1A_IW_SLC__1SDH_20191216T194511_20191216T194536_030378_0379CF_9F82.SAFE

In particular it has been used the annotation file:

  • annotations/s1a-iw1-slc-hh-20191216t194512-20191216t194534-030378-0379cf-001.xml

[5]:
swath_name = 'IW1'
burst_t0 = datetime.datetime.fromisoformat('2019-12-16T19:45:20.475893')
burst_tau0 = 5.372502580223076e-03
burst_dt = 2.055556299999998e-03
burst_range_pixel_spacing = 2.329562e+00
burst_dtau = burst_range_pixel_spacing / (c / 2)
burst_lines = 1503
burst_samples = 20701

burst_duration = burst_lines * burst_dt
burst_center = burst_t0 + datetime.timedelta(seconds=0.5 * burst_lines * burst_dt)
t_ref = burst_center
[6]:
print('Burst parameters')
print('t0:            ', burst_t0)
print('tau0:          ', burst_tau0)
print('dt:            ', burst_dt)
print('dtau:          ', burst_dtau)
print('lines:         ', burst_lines)
print('sammples:      ', burst_samples)
print('burst duration:', burst_duration)
print('referenc time: ', t_ref, '(burst center, arbitrary choice)')
Burst parameters
t0:             2019-12-16 19:45:20.475893
tau0:           0.005372502580223076
dt:             0.002055556299999998
dtau:           1.554116481475995e-08
lines:          1503
sammples:       20701
burst duration: 3.089501118899997
referenc time:  2019-12-16 19:45:22.020644 (burst center, arbitrary choice)

2. Extract the relevant timing information from the ETAD product

The NetCDF file is structured in hierarchical groups. The first level contains groups for each sub-swath (i.e. IW1-IW3 for IW mode). Each swath group is further divided into groups for the individual Bursts which are named and ordered according to their acquisition start time.

The overall organization of the product is reflected by the Python API provided by the s1etad package, but internal details of XML and NetCDF files are hidden.

Load S1-ETAD product

Please refer to s1etad quick start guide for more information about the installation procedure of the s1etad Python package.

[7]:
import s1etad
[8]:
filename = 'S1A_IW_ETA__AXDH_20191216T194148_20191216T194536_030378_0379CF_2705.SAFE'
[9]:
eta = s1etad.Sentinel1Etad(filename)
[10]:
eta
[10]:
Sentinel1Etad("S1A_IW_ETA__AXDH_20191216T194148_20191216T194536_030378_0379CF_2705.SAFE")  # 0x7f92565708d0
Number of Sentinel-1 slices: 9
Sentinel-1 products list:
  S1A_IW_SLC__1SDH_20191216T194148_20191216T194217_030378_0379CF_C6E4.SAFE
  S1A_IW_SLC__1SDH_20191216T194215_20191216T194243_030378_0379CF_D303.SAFE
  S1A_IW_SLC__1SDH_20191216T194241_20191216T194308_030378_0379CF_030E.SAFE
  S1A_IW_SLC__1SDH_20191216T194306_20191216T194333_030378_0379CF_0890.SAFE
  S1A_IW_SLC__1SDH_20191216T194331_20191216T194358_030378_0379CF_E79A.SAFE
  S1A_IW_SLC__1SDH_20191216T194355_20191216T194422_030378_0379CF_FD83.SAFE
  S1A_IW_SLC__1SDH_20191216T194420_20191216T194448_030378_0379CF_E077.SAFE
  S1A_IW_SLC__1SDH_20191216T194446_20191216T194513_030378_0379CF_4680.SAFE
  S1A_IW_SLC__1SDH_20191216T194511_20191216T194536_030378_0379CF_9F82.SAFE
Number of swaths: 3
Swath list: IW1, IW2, IW3
Azimuth time:
  min: 2019-12-16 19:41:48.058815
  max: 2019-12-16 19:45:36.583231
Range time:
  min: 0.005371694439612913
  max: 0.006416620248553707
Grid sampling:
  x: 8.081406101630269e-07
  y: 0.028777788199999974
  unit: s
Grid spacing:
  x: 200.0
  y: 200.0
  unit: m
Processing settings:
  troposphericDelayCorrection: True
  ionosphericDelayCorrection: True
  solidEarthTideCorrection: True
  bistaticAzimuthCorrection: True
  dopplerShiftRangeCorrection: True
  FMMismatchAzimuthCorrection: True

Query ETAD product

[11]:
margin = 0.5
t0_query = burst_center - datetime.timedelta(seconds=(burst_duration + margin) / 2)
t1_query = burst_center + datetime.timedelta(seconds=(burst_duration + margin) / 2)
[12]:
selection = eta.query_burst(swath='IW1', first_time=t0_query, last_time=t1_query)
[13]:
selection
[13]:
bIndex pIndex sIndex productID swathID azimuthTimeMin azimuthTimeMax
77 232 9 1 S1A_IW_SLC__1SDH_20191216T194511_20191216T1945... IW1 2019-12-16 19:45:20.438891915 2019-12-16 19:45:23.604448617

Get the burst object

[14]:
swath = eta[swath_name]
burst = swath[selection.bIndex.values[0]]
[15]:
burst
[15]:
Sentinel1EtadBurst("/IW1/Burst0232")  0x7f923de715d0
Swaths ID: IW1
Burst index: 232
Shape: (111, 402)
Sampling start:
  x: 0.0
  y: 212.3800769159272
  units: s
Sampling:
  x: 8.081406101630269e-07
  y: 0.028777788199999974
  units: s

Compute ETAD time grids for the selected burst

Start time computation

The burst group contains range and azimuth start times which need to be calculated:

Az = rootGroup:azimuthTimeMin + burst:gridStartAzimuthTime
Rg = rootGroup:rangeTimeMin + burst:gridStartRangeTime

The above formula corresponds to the following expression when one uses the s1etad Python API:

[16]:
eta_az_start = (eta.min_azimuth_time - t_ref).total_seconds() + burst.sampling_start['y']
eta_rg_start = eta.min_range_time                             + burst.sampling_start['x']

NOTE: please note that the t_ref is used to be able to express eta_az_ax_start as a float.

[17]:
print('eta_az_start:', eta_az_start, 'sec (relative to "t_ref")')
print('eta_rg_start:', eta_rg_start, 'sec')
eta_az_start: -1.5817520840728037 sec (relative to "t_ref")
eta_rg_start: 0.005371694439612913 sec

Grid axis computation

[18]:
eta_az_ax = eta_az_start + np.arange(burst.lines)   * burst.sampling['y']
eta_rg_ax = eta_rg_start + np.arange(burst.samples) * burst.sampling['x']

A simpler alternative is to use the burst.get_burst_grid method:

[19]:
eta_az_ax_rel, eta_rg_ax_rel = burst.get_burst_grid()
[20]:
eta_az_ax2 = (eta.min_azimuth_time - t_ref).total_seconds() + eta_az_ax_rel
eta_rg_ax2 = eta.min_range_time                             + eta_rg_ax_rel

NOTE: in this case the burst.sampling_start is not used because it is already taken into account in burst.get_burst_grid.

Extraction of correction data

The relevant timing corrections are named (in the NetCDF file):

  • sumOfCorrectionsRg

  • sumOfCorrectionsAz

They can be accessed with the Python API as follows:

[21]:
corrections = burst.get_correction(s1etad.ECorrectionType.SUM)
rg_correction = corrections['x']
az_correction = corrections['y']
[22]:
fig, ax = plt.subplots(1, 2, figsize=[13, 5], sharex=True, sharey=True)

img = ax[0].imshow(rg_correction * 1e9, aspect='auto', origin='lower')
ax[0].set_title('Sum of Range Corrections')
ax[0].set_xlabel('Range [samples]')
ax[0].set_ylabel('Azimuth [lines]')
fig.colorbar(img, ax=ax[0]).set_label(r'$[ns]$')

img = ax[1].imshow(az_correction * 1e3, aspect='auto', origin='lower')
ax[1].set_title('Sum of Azimuth Corrections')
ax[1].set_xlabel('Range [samples]')
# ax[1].set_ylabel('Azimuth [lines]')
fig.colorbar(img, ax=ax[1]).set_label(r'$[ms]$')
../_images/notebooks_S1-ETAD-Product_Quick_User_Guide_44_0.png

3. Resample the regularly sampled ETAD burst grids onto the regularly sampled S1-SLC grid with a bi-linear interpolation

Azimuth times \(t\) and range times \(\tau\) of the burst grid are defined by the following NetCDF parameters:

Az = rootGroup:azimuthTimeMin + burst:gridStartAzimuthTime + burst:azimuth
Rg = rootGroup:rangeTimeMin   + burst:gridStartRangeTime   + burst:range

The relevant ETAD NetCDF corrections which need to be applied are:

  • sumOfCorrectionsRg: \(\Delta_{\tau}(t, \tau)\)

  • sumOfCorrectionsAz: \(\Delta_{t}(t, \tau)\)

Computation of the regularly sampled S1-SLC grid

[23]:
burst_t0_sec = (burst_t0 - t_ref).total_seconds()
slc_az_ax = burst_t0_sec + np.arange(burst_lines)   * burst_dt
slc_rg_ax = burst_tau0   + np.arange(burst_samples) * burst_dtau

slc_rg, slc_az = np.meshgrid(slc_rg_ax, slc_az_ax)
[24]:
eta_box = np.asarray([
    (eta_rg_ax[0],  eta_az_ax[0]),
    (eta_rg_ax[0],  eta_az_ax[-1]),
    (eta_rg_ax[-1], eta_az_ax[-1]),
    (eta_rg_ax[-1], eta_az_ax[0]),
    (eta_rg_ax[0],  eta_az_ax[0]),
])
slc_box = np.asarray([
    (slc_rg_ax[0],  slc_az_ax[0]),
    (slc_rg_ax[0],  slc_az_ax[-1]),
    (slc_rg_ax[-1], slc_az_ax[-1]),
    (slc_rg_ax[-1], slc_az_ax[0]),
    (slc_rg_ax[0],  slc_az_ax[0]),
])

fix, ax = plt.subplots()
ax.plot(eta_box[:, 0] * 1e3, eta_box[:, 1], 'r', label='eta')
ax.plot(slc_box[:, 0] * 1e3, slc_box[:, 1], 'b', label='slc')
ax.set_xlabel('Range [ms]')
ax.set_ylabel('Azimuth [s]')
ax.legend()
ax.grid()
../_images/notebooks_S1-ETAD-Product_Quick_User_Guide_49_0.png

Grid resampling

[25]:
from scipy.interpolate import RectBivariateSpline

rg_corr_interpolator = RectBivariateSpline(eta_az_ax, eta_rg_ax, rg_correction, kx=1, ky=1)  # bi-linear
az_corr_interpolator = RectBivariateSpline(eta_az_ax, eta_rg_ax, az_correction, kx=1, ky=1)  # bi-linear

delta_tau = rg_corr_interpolator(slc_az_ax, slc_rg_ax)
delta_t   = az_corr_interpolator(slc_az_ax, slc_rg_ax)

4. Subtract the resampled range corrections \(\Delta_{\tau}(t, \tau)\) from the S1-SLC range times

[26]:
slc_rg -= delta_tau

5. Subtract the resampled azimuth corrections \(\Delta_{t}(t, \tau)\) from the S1-SLC azimuth times

[27]:
slc_az -= delta_t

6. The result is an irregularly spaced burst grid with corrected timings for the S1 SLC data

NOTE: dummy synthetic (float) data are used just to demonstrate the procedure.

[28]:
burst_slc = np.linspace(0, 1, burst_lines).reshape(burst_lines, 1) + np.linspace(0, 1, burst_samples)

NOTE: in this example it is used a standard scipy interpolator. In a real use case a proper interpolator for complex interferometric data shall be used.

NOTE: the interpolation step is very slow. To reduce the computation time data are sub-sampled in the range direction.

[29]:
from scipy.interpolate import RectBivariateSpline

subsamp = 10
slc_interpolator = RectBivariateSpline(slc_az_ax, slc_rg_ax[::subsamp], burst_slc[:, ::subsamp])
[30]:
interpolated_slc = slc_interpolator(slc_az[:, ::subsamp], slc_rg[:, ::subsamp], grid=False)
[31]:
extent = (0, burst_samples, burst_lines, 0)

fig, ax = plt.subplots(1, 2, figsize=[13, 5], sharex=True, sharey=True)

img = ax[0].imshow(burst_slc[:,::subsamp], aspect='auto', origin='lower', extent=extent)
ax[0].set_title('Oridinal SLC (dummy) data')
ax[0].set_xlabel('Range [samles]')
ax[0].set_ylabel('Azimuth [lines]')
ax[0].grid()
fig.colorbar(img, ax=ax[0]).set_label('Amplitude')

img = ax[1].imshow(burst_slc[:,::subsamp] - interpolated_slc, aspect='auto', origin='lower', extent=extent)
ax[1].set_title(r'$\Delta$ = data - corrected_data')
ax[1].set_xlabel('Range [samles]')
# ax[1].set_ylabel('Azimuth [lines]')
ax[1].grid()
cb = fig.colorbar(img, ax=ax[1])
cb.set_label(r'$\Delta$')
cb.formatter.set_scientific(True)
cb.formatter.set_powerlimits((0, 0))
../_images/notebooks_S1-ETAD-Product_Quick_User_Guide_63_0.png