Use case 4: corner reflector

This section explains how to identify bursts in which the corner reflector (CR), having known geodetic coordinates, is visible, and how to get corrections at CR location.

Notebook setup

%matplotlib inline
%load_ext autoreload
%autoreload 2

import sys

sys.path.append("../..")

import numpy as np
from scipy.interpolate import RegularGridInterpolator
from matplotlib import pyplot as plt

import s1etad
from s1etad import Sentinel1Etad

Searching bursts in which the CR is present

Load the S1-ETAD product.

filename = (
    "data/"
    "S1B_IW_ETA__AXDV_20200124T221416_20200124T221444_019964_025C43_0A63.SAFE"
)

eta = Sentinel1Etad(filename)

The CR position has been chosen to be in the overlap region of different bursts and swaths.

from shapely.geometry import Point

lat0 = -0.38  #   0°22'40.16"S
lon0 = -59.37  # 59°23'20.70"W
h0 = 223.57

cr = Point(lon0, lat0)

Query for burst covering the CR.

selection = eta.query_burst(geometry=cr)
selection

	bIndex	pIndex	sIndex	productID	swathID	azimuthTimeMin	azimuthTimeMax
14	14	1	2	S1B_IW_SLC__1ADV_20200124T221416_20200124T2214...	IW2	2020-01-24 22:14:28.445747384	2020-01-24 22:14:31.612902809
5	16	1	1	S1B_IW_SLC__1ADV_20200124T221416_20200124T2214...	IW1	2020-01-24 22:14:30.263929202	2020-01-24 22:14:33.401759114
15	17	1	2	S1B_IW_SLC__1ADV_20200124T221416_20200124T2214...	IW2	2020-01-24 22:14:31.202345624	2020-01-24 22:14:34.369501049

Get corrections at CR location

Retrieve the first burst.

burst = next(eta.iter_bursts(selection))

Get RADAR coordinates (range and azimuth time) of the CR.

tau0, t0 = burst.geodetic_to_radar(lat0, lon0, h0)
tau0, t0

(array([0.00032243]), array([15.04354041]))

Get image coordinates (line and sample) of the CR.

NOTE: image coordinates area floating point numbers including pixel fractions.

line, sample = burst.radar_to_image(t0, tau0)
print(f"Line {line} of {burst.lines}, sample {sample} of {burst.samples}.")

Line [104.98472789] of 109, sample [14.51639925] of 476.

Get the correction.

correction = burst.get_correction(s1etad.ECorrectionType.SUM, meter=True)
data = correction["x"]  # summ of corrections in the range direction

Interpolate the correction at the specified RADAR coordinates (tau0, t0).

azimuth_time, range_time = burst.get_burst_grid()

interpolator = RegularGridInterpolator((range_time, azimuth_time), data.T)

value_at_radar_coordinates = interpolator((tau0, t0))

print(
    f"The correction value at RADAR coordinates (tau0, t0) = ({tau0}, {t0}) is {value_at_radar_coordinates}"
)

The correction value at RADAR coordinates (tau0, t0) = ([0.00032243], [15.04354041]) is [3.17136665]

Interpolate the correction at the specified image coordinates (line, sample).

xaxis = np.arange(burst.samples)
yaxis = np.arange(burst.lines)

interpolator = RegularGridInterpolator((xaxis, yaxis), data.T)

value_at_image_coordinates = interpolator((sample, line))

print(
    f"The correction value at image coordinates (sample, line) = ({sample}, {line}) is {value_at_image_coordinates}"
)

The correction value at image coordinates (sample, line) = ([14.51639925], [104.98472789]) is [3.17136665]

Putting all together

from matplotlib.patches import Rectangle

fig, ax = plt.subplots(nrows=len(selection), ncols=1, figsize=[13, 8])

if len(selection) < 2:
    ax = [ax]

for loop, burst in enumerate(eta.iter_bursts(selection)):
    # get RADAR and image coordinates
    tau0, t0 = burst.geodetic_to_radar(lat0, lon0)
    line0, sample0 = burst.radar_to_image(t0, tau0)
    print(
        "xy  ",
        burst.swath_id,
        burst.burst_index,
        line0,
        sample0,
        burst.radar_to_geodetic(tau0, t0),
    )
    print(
        "time",
        burst.swath_id,
        burst.burst_index,
        tau0,
        t0,
        burst.radar_to_geodetic(tau0, t0),
    )

    # correction
    cor = burst.get_correction(s1etad.ECorrectionType.SUM, meter="True")

    ax[loop].imshow(cor["x"], aspect="auto")
    ax[loop].set_ylabel(f"{burst.swath_id}/{burst.burst_id}")
    rec_half_size = 1
    p = Rectangle(
        (sample0 - rec_half_size, line0 - rec_half_size),
        width=rec_half_size * 2 + 1,
        height=rec_half_size * 2 + 1,
        color="red",
        fill=True,
    )
    ax[loop].add_patch(p)

    # get the range and azimuth time axes
    azimuth_time, range_time = burst.get_burst_grid()

    # interpolate at the desired working (RADAR) coordinates
    f_t = RegularGridInterpolator((range_time, azimuth_time), cor["x"].T)

    # get the image (lines and samples) axes
    yaxis = np.arange(azimuth_time.size)
    xaxis = np.arange(range_time.size)

    # interpolate at the desired working (image) coordinates
    f_ij = RegularGridInterpolator((xaxis, yaxis), cor["x"].T)

    print(
        f"Interpolation by array coordinate {f_ij((sample0, line0))} or time {f_t((tau0, t0))} should be the same"
    )
    print(
        f"The total correction at lat/lon {lat0, lon0} is  {f_ij((sample0, line0))} m in range"
    )
    print()

xy   IW2 14 [104.98472789] [14.51639925] (array([-0.38]), array([-59.37]), array([223.56783955]))
time IW2 14 [0.00032243] [15.04354041] (array([-0.38]), array([-59.37]), array([223.56783955]))
Interpoaltion by array coordinate [3.17136665] or time [3.17136665] should be the same
The total correction at lat/lon (-0.38, -59.37) is  [3.17136665] m in range

xy   IW1 16 [42.98307299] [396.51924071] (array([-0.38]), array([-59.37]), array([223.48899634]))
time IW1 16 [0.00032244] [15.04349188] (array([-0.38]), array([-59.37]), array([223.48899634]))
Interpoaltion by array coordinate [3.6754839] or time [3.6754839] should be the same
The total correction at lat/lon (-0.38, -59.37) is  [3.6754839] m in range

xy   IW2 17 [10.98391419] [14.51925819] (array([-0.38]), array([-59.37]), array([223.51050639]))
time IW2 17 [0.00032244] [15.04351655] (array([-0.38]), array([-59.37]), array([223.51050639]))
Interpoaltion by array coordinate [3.96727542] or time [3.96727542] should be the same
The total correction at lat/lon (-0.38, -59.37) is  [3.96727542] m in range

burst.radar_to_geodetic(tau0, t0)

(array([-0.38]), array([-59.37]), array([223.51050639]))

cor

{'x': array([[4.04582886, 4.05006726, 4.06043874, ..., 4.28537005, 4.27479399,
         4.27751953],
        [4.05717985, 4.05634625, 4.05804824, ..., 4.27876634, 4.26313737,
         4.26012205],
        [4.03887094, 4.04289338, 4.04641388, ..., 4.27190359, 4.27784639,
         4.27295177],
        ...,
        [3.17816942, 3.18652249, 3.19618868, ..., 3.40269193, 3.41371017,
         3.40664538],
        [3.16365578, 3.16691524, 3.17739579, ..., 3.39818295, 3.40377122,
         3.40371275],
        [3.16692545, 3.17197457, 3.17454612, ..., 3.38868677, 3.39083797,
         3.39346216]]),
 'y': array([[-2.01295019, -2.00427015, -1.98119158, ..., -3.05862508,
         -3.08279639, -3.0770417 ],
        [-1.96329126, -1.96693968, -1.96464479, ..., -3.05755989,
         -3.09247042, -3.0996791 ],
        [-1.98382678, -1.97638487, -1.97001614, ..., -3.05710019,
         -3.04477515, -3.0558878 ],
        ...,
        [-1.82117034, -1.8436438 , -1.86872783, ..., -3.27037334,
         -3.29933656, -3.28893836],
        [-1.80620571, -1.81735789, -1.84509837, ..., -3.2799184 ,
         -3.29739384, -3.30231877],
        [-1.83240352, -1.84772252, -1.85734387, ..., -3.27850086,
         -3.28845637, -3.29956304]]),
 'unit': 'm',
 'name': 'sum'}