Get corrections at a corner reflector (CR) location using GridGeococoding
[1]:
%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
from s1etad.geometry import GridGeocoding
Searching bursts in which the CR is present
Load the S1-ETAD product.
[2]:
filename = (
"data/"
"S1B_IW_ETA__AXDV_20200124T221416_20200124T221444_019964_025C43_0A63.SAFE"
)
[3]:
eta = Sentinel1Etad(filename)
The CR position has been chosen to be in the overlap region of different bursts and swaths.
[4]:
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.
[5]:
selection = eta.query_burst(geometry=cr)
selection
[5]:
| 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.
[6]:
burst = next(eta.iter_bursts(selection))
Get the grid of geodetic coordinates.
[7]:
lats, lons, heights = burst.get_lat_lon_height()
Get the range and azimuth time axes.
[8]:
azimuth_time, range_time = burst.get_burst_grid()
Initialize the Grid Geocoding object.
NOTE: if one uses time axis then consistent time coordinates shall be provided in all back-geocoding requests.
[9]:
ebg = GridGeocoding(lats, lons, heights, xaxis=range_time, yaxis=azimuth_time)
Now it is possible to perform the back-geocoding i.e. computation of RADAR coordinates (tau, t) starting form geodetic coordinates (lat, lon, h):
(lat, lon, h) -> (tau, t)
[10]:
tau, t = ebg.backward_geocode(lat0, lon0, h0)
tau, t
[10]:
(array([0.00032243]), array([15.04354041]))
Of course it is also possible to make the inverse conversion:
[11]:
lat1, lon1, h1 = ebg.forward_geocode(tau, t)
print(f"Initial coordinates: (lat0, lon0, h0) = ({lat0}, {lon0}, {h0})")
print(
f"Foeward geocoding output: (lat1, lon1, h1) = ({lat1.item()}, {lon1.item()}, {h1.item()})"
)
Initial coordinates: (lat0, lon0, h0) = (-0.38, -59.37, 223.57)
Foeward geocoding output: (lat1, lon1, h1) = (-0.3799999998146794, -59.36999999969689, 223.56783954882505)
Using image coordinates (lines, samples)
It is also possible to initialize the GridGeocoding without providing time axes information.
In this case the geocoder will work using image coordinates (lines and samples) instead of range/azimuth times.
[12]:
ebg = GridGeocoding(lats, lons, heights)
It is possible to perform back-geocoding, i.e. (lat, lon, h) -> (sample, line):
[13]:
sample, line = ebg.backward_geocode(lat0, lon0, h0)
sample, line
[13]:
(array([14.51639907]), array([104.98472783]))
and also to perform the forward conversion: (sample, line) -> (lat, lon, h)
[14]:
lat1, lon1, h1 = ebg.forward_geocode(sample, line)
print(f"Initial coordinates: (lat0, lon0, h0) = ({lat0}, {lon0}, {h0})")
print(
f"Foeward geocoding output: (lat1, lon1, h1) = ({lat1.item()}, {lon1.item()}, {h1.item()})"
)
Initial coordinates: (lat0, lon0, h0) = (-0.38, -59.37, 223.57)
Foeward geocoding output: (lat1, lon1, h1) = (-0.37999999999661643, -59.369999999994484, 223.5678411001362)
Putting all together
[15]:
from matplotlib.patches import Rectangle
fig, ax = plt.subplots(nrows=len(selection), ncols=1, figsize=[13, 8])
for loop, burst in enumerate(eta.iter_bursts(selection)):
# coordinate grids
lats, lons, heights = burst.get_lat_lon_height()
# backward geocoding with image coordinates
ebg = GridGeocoding(lats, lons, heights)
x0, y0 = ebg.backward_geocode(lat0, lon0, h0)
print(
"xy ",
burst.swath_id,
burst.burst_index,
x0,
y0,
ebg.latitude(x0, y0),
ebg.longitude(x0, y0),
ebg.height(x0, y0),
)
# get the range and azimuth times
azimuth_time, range_time = burst.get_burst_grid()
# backward geocoding with time coordinates
ebg = GridGeocoding(
lats, lons, heights, xaxis=range_time, yaxis=azimuth_time
)
tau0, t0 = ebg.backward_geocode(lat0, lon0, h0)
print(
"time",
burst.swath_id,
burst.burst_index,
tau0,
t0,
ebg.latitude(tau0, t0),
ebg.longitude(tau0, t0),
ebg.height(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(
(x0 - rec_half_size, y0 - rec_half_size),
width=rec_half_size * 2 + 1,
height=rec_half_size * 2 + 1,
color="red",
fill=True,
)
ax[loop].add_patch(p)
# 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((x0, y0))} or time {f_t((tau0, t0))} should be the same"
)
print(
f"The total correction at lat/lon {lat0, lon0} is {f_ij((x0, y0))} m in range"
)
print()
xy IW2 14 [14.51639907] [104.98472783] [-0.38] [-59.37] [223.5678411]
time IW2 14 [0.00032243] [15.04354041] [-0.38] [-59.37] [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 [396.51924056] [42.98307295] [-0.38] [-59.37] [223.48899793]
time IW1 16 [0.00032244] [15.04349188] [-0.38] [-59.37] [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 [14.51925805] [10.98391415] [-0.38] [-59.37] [223.51050799]
time IW2 17 [0.00032244] [15.04351655] [-0.38] [-59.37] [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
[16]:
burst.radar_to_geodetic(tau0, t0)
[16]:
(array([-0.38]), array([-59.37]), array([223.51050639]))
[17]:
cor
[17]:
{'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'}