From d8f9e80364c02ef338062525065b1122599075d6 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Tue, 22 Nov 2022 20:39:28 -0500 Subject: [PATCH] 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 --- MANIFEST.in | 2 +- src/s1reader/s1_burst_id.py | 140 ++++++++++++++++++++++++++++++++++ src/s1reader/s1_burst_slc.py | 19 +++-- src/s1reader/s1_reader.py | 143 ++++++++++++++++++++--------------- src/s1reader/version.py | 1 - tests/test_bursts.py | 2 +- 6 files changed, 237 insertions(+), 70 deletions(-) create mode 100644 src/s1reader/s1_burst_id.py diff --git a/MANIFEST.in b/MANIFEST.in index c2b49866..5ff768ec 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1 +1 @@ -include src/s1reader/data/sentinel1_track_burst_id.txt +recursive-include src/s1reader/data/ * diff --git a/src/s1reader/s1_burst_id.py b/src/s1reader/s1_burst_id.py new file mode 100644 index 00000000..1b8b6709 --- /dev/null +++ b/src/s1reader/s1_burst_id.py @@ -0,0 +1,140 @@ +import datetime +from dataclasses import dataclass +from typing import ClassVar + +import numpy as np + + +@dataclass(frozen=True) +class S1BurstId: + # Constants in Table 9-7 of Sentinel-1 SLC Detailed Algorithm Definition + T_beam: ClassVar[float] = 2.758273 # interval of one burst [s] + T_pre: ClassVar[float] = 2.299849 # Preamble time interval [s] + T_orb: ClassVar[float] = 12 * 24 * 3600 / 175 # Nominal orbit period [s] + track_number: int + esa_burst_id: int + subswath: str + + @classmethod + def from_burst_params( + cls, + sensing_time: datetime.datetime, + ascending_node_dt: datetime.datetime, + start_track: int, + end_track: int, + subswath: str, + ): + """Calculate the unique burst ID (track, ESA burstId, swath) of a burst. + + Accounts for equator crossing frames, where the current track number of + a burst may change mid-frame. Uses the ESA convention defined in the + Sentinel-1 Level 1 Detailed Algorithm Definition. + + Parameters + ---------- + sensing_time : datetime + Sensing time of the first input line of this burst [UTC] + The XML tag is sensingTime in the annotation file. + ascending_node_dt : datetime + Time of the ascending node prior to the start of the scene. + start_track : int + Relative orbit number at the start of the acquisition, from 1-175. + end_track : int + Relative orbit number at the end of the acquisition. + subswath : str, {'IW1', 'IW2', 'IW3'} + Name of the subswath of the burst (not case sensitive). + + Returns + ------- + S1BurstId + The burst ID object containing track number + ESA's burstId number + swath ID. + + Notes + ----- + The `start_track` and `end_track` parameters are used to determine if the + scene crosses the equator. They are the same if the frame does not cross + the equator. + + References + ---------- + ESA Sentinel-1 Level 1 Detailed Algorithm Definition + 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 + """ + swath_num = int(subswath[-1]) + # Since we only have access to the current subswath, we need to use the + # burst-to-burst times to figure out + # 1. if IW1 crossed the equator, and + # 2. The mid-burst sensing time for IW2 + # IW1 -> IW2 takes ~0.83220 seconds + # IW2 -> IW3 takes ~1.07803 seconds + # IW3 -> IW1 takes ~0.84803 seconds + burst_times = np.array([0.832, 1.078, 0.848]) + iw1_start_offsets = [ + 0, + -burst_times[0], + -burst_times[0] - burst_times[1], + ] + offset = iw1_start_offsets[swath_num - 1] + start_iw1 = sensing_time + datetime.timedelta(seconds=offset) + + start_iw1_to_mid_iw2 = burst_times[0] + burst_times[1] / 2 + mid_iw2 = start_iw1 + datetime.timedelta(seconds=start_iw1_to_mid_iw2) + + has_anx_crossing = end_track == (start_track + 1) % 175 + + time_since_anx_iw1 = (start_iw1 - ascending_node_dt).total_seconds() + time_since_anx = (mid_iw2 - ascending_node_dt).total_seconds() + + if (time_since_anx_iw1 - cls.T_orb) < 0: + # Less than a full orbit has passed + track_number = start_track + else: + track_number = end_track + # Additional check for scenes which have a given ascending node + # that's more than 1 orbit in the past + if not has_anx_crossing: + time_since_anx = time_since_anx - cls.T_orb + + # Eq. 9-89: ∆tb = tb − t_anx + (r - 1)T_orb + # tb: mid-burst sensing time (sensing_time) + # t_anx: ascending node time (ascending_node_dt) + # r: relative orbit number (relative_orbit_start) + dt_b = time_since_anx + (start_track - 1) * cls.T_orb + + # Eq. 9-91 : 1 + floor((∆tb − T_pre) / T_beam ) + esa_burst_id = 1 + int(np.floor((dt_b - cls.T_pre) / cls.T_beam)) + + return cls(track_number, esa_burst_id, subswath) + + @classmethod + def from_str(cls, burst_id_str: str): + """Parse a S1BurstId object from a string. + + Parameters + ---------- + burst_id_str : str + The burst ID string, e.g. "t123_456_iw1" + + Returns + ------- + S1BurstId + The burst ID object containing track number + ESA's burstId number + swath ID. + """ + track_number, esa_burst_id, subswath = burst_id_str.split("_") + track_number = int(track_number[1:]) + esa_burst_id = int(esa_burst_id) + return cls(track_number, esa_burst_id, subswath.lower()) + + def __str__(self): + # Form the unique JPL ID by combining track/burst/swath + return f"t{self.track_number:03d}_{self.esa_burst_id:06d}_{self.subswath.lower()}" + + def __eq__(self, other) -> bool: + # Allows for comparison with strings, as well as S1BurstId objects + # e.g., you can filter down burst IDs with: + # burst_ids = ["t012_024518_iw3", "t012_024519_iw3"] + # bursts = [b for b in bursts if b.burst_id in burst_ids] + if isinstance(other, str): + return str(self) == other + else: + return super().__eq__(other) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 35a25d1a..2c68330a 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -10,6 +10,8 @@ from osgeo import gdal from s1reader import s1_annotation +from .s1_burst_id import S1BurstId + # Other functionalities def polyfit(xin, yin, zin, azimuth_order, range_order, @@ -144,7 +146,7 @@ class Sentinel1BurstSlc: doppler: Doppler range_bandwidth: float polarization: str # {VV, VH, HH, HV} - burst_id: str # t{track_number}_{burst_index}_iw{1,2,3} + burst_id: S1BurstId platform_id: str # S1{A,B} safe_filename: str # SAFE file name center: tuple # {center lon, center lat} in degrees @@ -363,6 +365,8 @@ def as_dict(self): temp['std'] = val.std temp['coeffs'] = val.coeffs val = temp + elif key == 'burst_id': + val = str(val) elif key == 'border': val = self.border[0].wkt elif key == 'doppler': @@ -664,7 +668,7 @@ def width(self): @property def swath_name(self): '''Swath name in iw1, iw2, iw3.''' - return self.burst_id.split('_')[2] + return self.burst_id.swath_name @property def thermal_noise_lut(self): @@ -690,8 +694,9 @@ def eap_compensation_lut(self): f' IPF version = {self.ipf_version}') return self.burst_eap.compute_eap_compensation_lut(self.width) - def bbox(self): - '''Returns the (west, south, east, north) bounding box of the burst.''' - # Uses https://shapely.readthedocs.io/en/stable/manual.html#object.bounds - # Returns a tuple of 4 floats representing (west, south, east, north) in degrees. - return self.border[0].bounds + + @property + def relative_orbit_number(self): + '''Returns the relative orbit number of the burst.''' + orbit_number_offset = 73 if self.platform_id == 'S1A' else 202 + return (self.abs_orbit_number - orbit_number_offset) % 175 + 1 diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index d6874f98..8e98de0d 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -3,6 +3,7 @@ import os import warnings import lxml.etree as ET +from typing import Union import zipfile from types import SimpleNamespace @@ -20,6 +21,7 @@ BurstCalibration, BurstEAP, BurstNoise from s1reader.s1_burst_slc import Doppler, Sentinel1BurstSlc +from s1reader.s1_burst_id import S1BurstId esa_track_burst_id_file = f"{os.path.dirname(os.path.realpath(__file__))}/data/sentinel1_track_burst_id.txt" @@ -32,9 +34,6 @@ def as_datetime(t_str): ---------- t_str : string Time string to be parsed. (e.g., "2021-12-10T12:00:0.0") - fmt : string - Format of string provided. Defaults to az time format found in annotation XML. - (e.g., "%Y-%m-%dT%H:%M:%S.%f"). Returns: ------ @@ -268,23 +267,62 @@ def get_burst_centers_and_boundaries(tree): return center_pts, boundary_pts def get_ipf_version(tree: ET): - '''Extract the IPF version from the ET of manifest.safe - ''' - # path to xmlData in manifest - xml_meta_path = 'metadataSection/metadataObject/metadataWrap/xmlData' + '''Extract the IPF version from the ET of manifest.safe. - # piecemeal build path to software path to access version attrib - esa_http = '{http://www.esa.int/safe/sentinel-1.0}' - processing = xml_meta_path + f'/{esa_http}processing' - facility = processing + f'/{esa_http}facility' - software = facility + f'/{esa_http}software' + Parameters + ---------- + tree : xml.etree.ElementTree + ElementTree containing the parsed 'manifest.safe' file. + Returns + ------- + ipf_version : version.Version + IPF version of the burst. + ''' # get version from software element - software_elem = tree.find(software) + search_term, nsmap = _get_manifest_pattern(tree, ['processing', 'facility', 'software']) + software_elem = tree.find(search_term, nsmap) ipf_version = version.parse(software_elem.attrib['version']) return ipf_version +def get_start_end_track(tree: ET): + '''Extract the start/end relative orbits from manifest.safe file''' + search_term, nsmap = _get_manifest_pattern(tree, ['orbitReference', 'relativeOrbitNumber']) + elem_start, elem_end = tree.findall(search_term, nsmap) + return int(elem_start.text), int(elem_end.text) + + +def _get_manifest_pattern(tree: ET, keys: list): + '''Get the search path to extract data from the ET of manifest.safe. + + Parameters + ---------- + tree : xml.etree.ElementTree + ElementTree containing the parsed 'manifest.safe' file. + keys : list + List of keys to search for in the manifest file. + + Returns + ------- + str + Search path to extract from the ET of the manifest.safe XML. + + ''' + # https://lxml.de/tutorial.html#namespaces + # Get the namespace from the root element to avoid full urls + # This is a dictionary with a short name containing the labels in the + # XML tree, and the fully qualified URL as the value. + try: + nsmap = tree.nsmap + except AttributeError: + nsmap = tree.getroot().nsmap + # path to xmlData in manifest + xml_meta_path = 'metadataSection/metadataObject/metadataWrap/xmlData' + safe_terms = "/".join([f'safe:{key}' for key in keys]) + return f'{xml_meta_path}/{safe_terms}', nsmap + + def get_path_aux_cal(directory_aux_cal: str, str_annotation: str): ''' Decide which aux_cal to load @@ -323,8 +361,7 @@ def get_path_aux_cal(directory_aux_cal: str, str_annotation: str): list_aux_cal = glob.glob(f'{directory_aux_cal}/{str_platform}_AUX_CAL_V*.SAFE.zip') if len(list_aux_cal) == 0: - raise ValueError( 'Cannot find AUX_CAL files from directory: ' - f'{directory_aux_cal}') + raise ValueError(f'Cannot find AUX_CAL files from {directory_aux_cal} .') format_datetime = '%Y%m%dT%H%M%S' @@ -457,27 +494,18 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, bursts : list List of Sentinel1BurstSlc objects found in annotation XML. ''' - - # a dict where the key is the track number and the value is a list of - # two integers for the start and stop burst number - track_burst_num = get_track_burst_num() - _, tail = os.path.split(annotation_path) - platform_id, swath_name, _, pol = [x.upper() for x in tail.split('-')[:4]] + platform_id, subswath, _, pol = [x.upper() for x in tail.split('-')[:4]] safe_filename = os.path.basename(annotation_path.split('.SAFE')[0]) - # For IW mode, one burst has a duration of ~2.75 seconds and a burst - # overlap of approximately ~0.4 seconds. - # https://sentinels.copernicus.eu/web/sentinel/user-guides/sentinel-1-sar/product-types-processing-levels/level-1 - # Additional precision calculated from averaging the differences between - # burst sensing starts in prototyping test data - burst_interval = 2.758277 - # parse manifest.safe to retrieve IPF version manifest_path = os.path.dirname(annotation_path).replace('annotation','') + 'manifest.safe' with open_method(manifest_path, 'r') as f_manifest: - tree_manfest = ET.parse(f_manifest) - ipf_version = get_ipf_version(tree_manfest) + tree_manifest = ET.parse(f_manifest) + ipf_version = get_ipf_version(tree_manifest) + # Parse out the start/end track to determine if we have an + # equator crossing (for the burst_id calculation). + start_track, end_track = get_start_end_track(tree_manifest) # Load the Product annotation - for EAP calibration with open_method(annotation_path, 'r') as f_lads: @@ -543,7 +571,7 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, slant_range_time = float(image_info_element.find('slantRangeTime').text) ascending_node_time = as_datetime(image_info_element.find('ascendingNodeTime').text) - downlink_element = tree.find('generalAnnotation/downlinkInformationList/downlinkInformation') + 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) @@ -565,8 +593,6 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, range_window_coeff = float(rng_processing_element.find('windowCoefficient').text) orbit_number = int(tree.find('adsHeader/absoluteOrbitNumber').text) - orbit_number_offset = 73 if platform_id == 'S1A' else 202 - track_number = (orbit_number - orbit_number_offset) % 175 + 1 center_pts, boundary_pts = get_burst_centers_and_boundaries(tree) @@ -590,31 +616,28 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, osv_list = [] # load individual burst + half_burst_in_seconds = 0.5 * (n_lines - 1) * azimuth_time_interval burst_list_elements = tree.find('swathTiming/burstList') n_bursts = int(burst_list_elements.attrib['count']) bursts = [[]] * n_bursts - sensing_starts = [[]] * n_bursts - sensing_times = [[]] * n_bursts for i, burst_list_element in enumerate(burst_list_elements): - # get burst timing + # Zero Doppler azimuth time of the first line of this burst sensing_start = as_datetime(burst_list_element.find('azimuthTime').text) - sensing_starts[i] = sensing_start - sensing_times[i] = as_datetime(burst_list_element.find('sensingTime').text) - dt = sensing_times[i] - ascending_node_time - # local burst_num within one track, starting from 0 - burst_num = int((dt.seconds + dt.microseconds / 1e6) // burst_interval) + azimuth_time_mid = sensing_start + datetime.timedelta(seconds=half_burst_in_seconds) - # convert the local burst_num to the global burst_num, starting from 1 - burst_num += track_burst_num[track_number][0] + # Sensing time of the first input line of this burst [UTC] + sensing_time = as_datetime(burst_list_element.find('sensingTime').text) + # Create the burst ID to match the ESA ID scheme + burst_id = S1BurstId.from_burst_params( + sensing_time, ascending_node_time, start_track, end_track, subswath + ) # choose nearest azimuth FM rate - d_seconds = 0.5 * (n_lines - 1) * azimuth_time_interval - sensing_mid = sensing_start + datetime.timedelta(seconds=d_seconds) - az_fm_rate = get_nearest_polynomial(sensing_mid, az_fm_rate_list) + az_fm_rate = get_nearest_polynomial(azimuth_time_mid, az_fm_rate_list) # choose nearest doppler - poly1d = get_nearest_polynomial(sensing_mid, doppler_list) + poly1d = get_nearest_polynomial(azimuth_time_mid, doppler_list) lut2d = doppler_poly1d_to_lut2d(poly1d, starting_range, range_pxl_spacing, (n_lines, n_samples), azimuth_time_interval) @@ -643,9 +666,6 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, last_sample = min(last_valid_samples[first_valid_line], last_valid_samples[last_line]) - - burst_id = f't{track_number:03d}_{burst_num}_{swath_name.lower()}' - # Extract burst-wise information for Calibration, Noise, and EAP correction if calibration_annotation is None: burst_calibration = None @@ -706,7 +726,7 @@ def _is_zip_annotation_xml(path: str, id_str: str) -> bool: def load_bursts(path: str, orbit_path: str, swath_num: int, pol: str = 'vv', - burst_ids: list[str] = None, + burst_ids: list[Union[str, S1BurstId]] = None, flag_apply_eap: bool = True): '''Find bursts in a Sentinel-1 zip file or a SAFE structured directory. @@ -720,11 +740,11 @@ def load_bursts(path: str, orbit_path: str, swath_num: int, pol: str = 'vv', Integer of subswath of desired burst. {1, 2, 3} pol : str Polarization of desired burst. {hh, vv, hv, vh} - burst_ids : list[str] + burst_ids : list[str] or list[S1BurstId] List of burst IDs for which their Sentinel1BurstSlc objects will be - returned. Default of None returns all bursts. Empty list returned if - none of the burst IDs are found. If not all burst IDs are found, a list - containing found bursts will be returned. + returned. Default of None returns all bursts. + If not all burst IDs are found, a list containing found bursts will be + returned (empty list if none are found). Returns: -------- @@ -738,7 +758,7 @@ def load_bursts(path: str, orbit_path: str, swath_num: int, pol: str = 'vv', burst_ids = [] # ensure burst IDs is a list - if not isinstance(burst_ids, list): + if isinstance(burst_ids, (str, S1BurstId)): burst_ids = [burst_ids] # lower case polarity to be consistent with file naming convention @@ -761,10 +781,12 @@ def load_bursts(path: str, orbit_path: str, swath_num: int, pol: str = 'vv', if burst_ids: bursts = [b for b in bursts if b.burst_id in burst_ids] - burst_ids_found = set([b.burst_id for b in bursts]) + # convert S1BurstId to string for consistent object comparison later + burst_ids_found = set([str(b.burst_id) for b in bursts]) warnings.simplefilter("always") - set_burst_ids = set(burst_ids) + # burst IDs as strings for consistent object comparison + set_burst_ids = set([str(b) for b in burst_ids]) if not burst_ids_found: warnings.warn("None of provided burst IDs found in sub-swath {swath_num}") elif burst_ids_found != set_burst_ids: @@ -785,7 +807,7 @@ def _burst_from_zip(zip_path: str, id_str: str, orbit_path: str, flag_apply_eap: path : str Path to zip file. id_str: str - Identifcation of desired burst. Format: iw[swath_num]-slc-[pol] + Identification of desired burst. Format: iw[swath_num]-slc-[pol] orbit_path : str Path the orbit file. @@ -819,6 +841,7 @@ def _burst_from_zip(zip_path: str, id_str: str, orbit_path: str, flag_apply_eap: flag_apply_eap=flag_apply_eap) return bursts + def _burst_from_safe_dir(safe_dir_path: str, id_str: str, orbit_path: str, flag_apply_eap: bool): '''Find bursts in a Sentinel-1 SAFE structured directory. @@ -827,7 +850,7 @@ def _burst_from_safe_dir(safe_dir_path: str, id_str: str, orbit_path: str, flag_ path : str Path to SAFE directory. id_str: str - Identifcation of desired burst. Format: iw[swath_num]-slc-[pol] + Identification of desired burst. Format: iw[swath_num]-slc-[pol] orbit_path : str Path the orbit file. diff --git a/src/s1reader/version.py b/src/s1reader/version.py index 31d1d6f8..812e3028 100644 --- a/src/s1reader/version.py +++ b/src/s1reader/version.py @@ -14,4 +14,3 @@ # latest release version number and date release_version = release_history[0].version release_date = release_history[0].date - diff --git a/tests/test_bursts.py b/tests/test_bursts.py index 11204126..29c95d75 100644 --- a/tests/test_bursts.py +++ b/tests/test_bursts.py @@ -28,7 +28,7 @@ def test_burst(bursts): [-2056.701472691132, 353389.9614836443, -54143009.57327797]] for i, burst in enumerate(bursts): - expected_burst_id = f't071_{151200 + i}_iw3' + expected_burst_id = f't071_{151199 + i}_iw3' assert burst.burst_id == expected_burst_id assert burst.i_burst == i assert burst.abs_orbit_number == 32518