Skip to content

Commit d8f9e80

Browse files
authored
Fix track number/burst id calculation (#77)
* add more small testing data scripts helped make very small version, orbits manually shrunk * make a failing test for anx crossing * create methods for burst/track, fill tests * fix erroneous burst id in existing test * add a sample burst db to compare * add script for remaking the burst sample in case we add more tests * add a geometry check for the esa database * perform test without pandas * codacy items * add the geometry comparison to the other test cases * add two more test cases * refactor tests for new cases * redo logic for strange track175 case * update burst csv * fix first test problem, cut bursts * get tests working again for track175 case * fix esa db csv * use nsmap instead of long manual urls * remove testing script * working version on full orbit cycle * fix tests to check all subswaths * try recursive include for circleci fail, codacy * retry manifest * add better docstrings, add print missing tiff again, comment for start/end track * case sensitivity clarification * revert tests folder back to current version * fix unit test for correct burst id * adjust anx crossing to use mod, add comment on lxml nsmap * split out burst id into class * fix burst load filtering for new class * use the class vars in initialization * formatting * just use __str__ instead of as_str, adjust as_dict * address comments for clarity
1 parent 647dd25 commit d8f9e80

File tree

6 files changed

+237
-70
lines changed

6 files changed

+237
-70
lines changed

MANIFEST.in

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
include src/s1reader/data/sentinel1_track_burst_id.txt
1+
recursive-include src/s1reader/data/ *

src/s1reader/s1_burst_id.py

Lines changed: 140 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,140 @@
1+
import datetime
2+
from dataclasses import dataclass
3+
from typing import ClassVar
4+
5+
import numpy as np
6+
7+
8+
@dataclass(frozen=True)
9+
class S1BurstId:
10+
# Constants in Table 9-7 of Sentinel-1 SLC Detailed Algorithm Definition
11+
T_beam: ClassVar[float] = 2.758273 # interval of one burst [s]
12+
T_pre: ClassVar[float] = 2.299849 # Preamble time interval [s]
13+
T_orb: ClassVar[float] = 12 * 24 * 3600 / 175 # Nominal orbit period [s]
14+
track_number: int
15+
esa_burst_id: int
16+
subswath: str
17+
18+
@classmethod
19+
def from_burst_params(
20+
cls,
21+
sensing_time: datetime.datetime,
22+
ascending_node_dt: datetime.datetime,
23+
start_track: int,
24+
end_track: int,
25+
subswath: str,
26+
):
27+
"""Calculate the unique burst ID (track, ESA burstId, swath) of a burst.
28+
29+
Accounts for equator crossing frames, where the current track number of
30+
a burst may change mid-frame. Uses the ESA convention defined in the
31+
Sentinel-1 Level 1 Detailed Algorithm Definition.
32+
33+
Parameters
34+
----------
35+
sensing_time : datetime
36+
Sensing time of the first input line of this burst [UTC]
37+
The XML tag is sensingTime in the annotation file.
38+
ascending_node_dt : datetime
39+
Time of the ascending node prior to the start of the scene.
40+
start_track : int
41+
Relative orbit number at the start of the acquisition, from 1-175.
42+
end_track : int
43+
Relative orbit number at the end of the acquisition.
44+
subswath : str, {'IW1', 'IW2', 'IW3'}
45+
Name of the subswath of the burst (not case sensitive).
46+
47+
Returns
48+
-------
49+
S1BurstId
50+
The burst ID object containing track number + ESA's burstId number + swath ID.
51+
52+
Notes
53+
-----
54+
The `start_track` and `end_track` parameters are used to determine if the
55+
scene crosses the equator. They are the same if the frame does not cross
56+
the equator.
57+
58+
References
59+
----------
60+
ESA Sentinel-1 Level 1 Detailed Algorithm Definition
61+
https://sentinels.copernicus.eu/documents/247904/1877131/S1-TN-MDA-52-7445_Sentinel-1+Level+1+Detailed+Algorithm+Definition_v2-4.pdf/83624863-6429-cfb8-2371-5c5ca82907b8
62+
"""
63+
swath_num = int(subswath[-1])
64+
# Since we only have access to the current subswath, we need to use the
65+
# burst-to-burst times to figure out
66+
# 1. if IW1 crossed the equator, and
67+
# 2. The mid-burst sensing time for IW2
68+
# IW1 -> IW2 takes ~0.83220 seconds
69+
# IW2 -> IW3 takes ~1.07803 seconds
70+
# IW3 -> IW1 takes ~0.84803 seconds
71+
burst_times = np.array([0.832, 1.078, 0.848])
72+
iw1_start_offsets = [
73+
0,
74+
-burst_times[0],
75+
-burst_times[0] - burst_times[1],
76+
]
77+
offset = iw1_start_offsets[swath_num - 1]
78+
start_iw1 = sensing_time + datetime.timedelta(seconds=offset)
79+
80+
start_iw1_to_mid_iw2 = burst_times[0] + burst_times[1] / 2
81+
mid_iw2 = start_iw1 + datetime.timedelta(seconds=start_iw1_to_mid_iw2)
82+
83+
has_anx_crossing = end_track == (start_track + 1) % 175
84+
85+
time_since_anx_iw1 = (start_iw1 - ascending_node_dt).total_seconds()
86+
time_since_anx = (mid_iw2 - ascending_node_dt).total_seconds()
87+
88+
if (time_since_anx_iw1 - cls.T_orb) < 0:
89+
# Less than a full orbit has passed
90+
track_number = start_track
91+
else:
92+
track_number = end_track
93+
# Additional check for scenes which have a given ascending node
94+
# that's more than 1 orbit in the past
95+
if not has_anx_crossing:
96+
time_since_anx = time_since_anx - cls.T_orb
97+
98+
# Eq. 9-89: ∆tb = tb − t_anx + (r - 1)T_orb
99+
# tb: mid-burst sensing time (sensing_time)
100+
# t_anx: ascending node time (ascending_node_dt)
101+
# r: relative orbit number (relative_orbit_start)
102+
dt_b = time_since_anx + (start_track - 1) * cls.T_orb
103+
104+
# Eq. 9-91 : 1 + floor((∆tb − T_pre) / T_beam )
105+
esa_burst_id = 1 + int(np.floor((dt_b - cls.T_pre) / cls.T_beam))
106+
107+
return cls(track_number, esa_burst_id, subswath)
108+
109+
@classmethod
110+
def from_str(cls, burst_id_str: str):
111+
"""Parse a S1BurstId object from a string.
112+
113+
Parameters
114+
----------
115+
burst_id_str : str
116+
The burst ID string, e.g. "t123_456_iw1"
117+
118+
Returns
119+
-------
120+
S1BurstId
121+
The burst ID object containing track number + ESA's burstId number + swath ID.
122+
"""
123+
track_number, esa_burst_id, subswath = burst_id_str.split("_")
124+
track_number = int(track_number[1:])
125+
esa_burst_id = int(esa_burst_id)
126+
return cls(track_number, esa_burst_id, subswath.lower())
127+
128+
def __str__(self):
129+
# Form the unique JPL ID by combining track/burst/swath
130+
return f"t{self.track_number:03d}_{self.esa_burst_id:06d}_{self.subswath.lower()}"
131+
132+
def __eq__(self, other) -> bool:
133+
# Allows for comparison with strings, as well as S1BurstId objects
134+
# e.g., you can filter down burst IDs with:
135+
# burst_ids = ["t012_024518_iw3", "t012_024519_iw3"]
136+
# bursts = [b for b in bursts if b.burst_id in burst_ids]
137+
if isinstance(other, str):
138+
return str(self) == other
139+
else:
140+
return super().__eq__(other)

src/s1reader/s1_burst_slc.py

Lines changed: 12 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,8 @@
1010
from osgeo import gdal
1111

1212
from s1reader import s1_annotation
13+
from .s1_burst_id import S1BurstId
14+
1315

1416
# Other functionalities
1517
def polyfit(xin, yin, zin, azimuth_order, range_order,
@@ -144,7 +146,7 @@ class Sentinel1BurstSlc:
144146
doppler: Doppler
145147
range_bandwidth: float
146148
polarization: str # {VV, VH, HH, HV}
147-
burst_id: str # t{track_number}_{burst_index}_iw{1,2,3}
149+
burst_id: S1BurstId
148150
platform_id: str # S1{A,B}
149151
safe_filename: str # SAFE file name
150152
center: tuple # {center lon, center lat} in degrees
@@ -363,6 +365,8 @@ def as_dict(self):
363365
temp['std'] = val.std
364366
temp['coeffs'] = val.coeffs
365367
val = temp
368+
elif key == 'burst_id':
369+
val = str(val)
366370
elif key == 'border':
367371
val = self.border[0].wkt
368372
elif key == 'doppler':
@@ -664,7 +668,7 @@ def width(self):
664668
@property
665669
def swath_name(self):
666670
'''Swath name in iw1, iw2, iw3.'''
667-
return self.burst_id.split('_')[2]
671+
return self.burst_id.swath_name
668672

669673
@property
670674
def thermal_noise_lut(self):
@@ -690,8 +694,9 @@ def eap_compensation_lut(self):
690694
f' IPF version = {self.ipf_version}')
691695

692696
return self.burst_eap.compute_eap_compensation_lut(self.width)
693-
def bbox(self):
694-
'''Returns the (west, south, east, north) bounding box of the burst.'''
695-
# Uses https://shapely.readthedocs.io/en/stable/manual.html#object.bounds
696-
# Returns a tuple of 4 floats representing (west, south, east, north) in degrees.
697-
return self.border[0].bounds
697+
698+
@property
699+
def relative_orbit_number(self):
700+
'''Returns the relative orbit number of the burst.'''
701+
orbit_number_offset = 73 if self.platform_id == 'S1A' else 202
702+
return (self.abs_orbit_number - orbit_number_offset) % 175 + 1

0 commit comments

Comments
 (0)