diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 412f064..69d05c0 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -201,114 +201,6 @@ def from_et(cls, et_in: ET, path_annotation: str): return cls -@dataclass -class NoiseAnnotation(AnnotationBase): - ''' - Reader for Noise Annotation Data Set (NADS) for IW SLC - Based on ESA documentation: "Thermal Denoising of Products Generated by the S-1 IPF" - ''' - - basename_annotation: str - rg_list_azimuth_time: np.ndarray - rg_list_line: list - rg_list_pixel: list - rg_list_noise_range_lut: list - az_first_azimuth_line: int - az_first_range_sample: int - az_last_azimuth_line: int - az_last_range_sample: int - az_line: np.ndarray - az_noise_azimuth_lut: np.ndarray - - @classmethod - def from_et(cls,et_in: ET, ipf_version: version.Version, path_annotation: str): - ''' - Extracts list of noise information from etree - - Parameter - ---------- - et_in : xml.etree.ElementTree - Parsed NADS annotation .xml - - Return - ------- - cls: NoiseAnnotation - Parsed NADS from et_in - ''' - - cls.xml_et = et_in - cls.basename_annotation = os.path.basename(path_annotation) - - if ipf_version < min_ipf_version_az_noise_vector: # legacy SAFE data - cls.rg_list_azimuth_time = \ - cls._parse_vectorlist('noiseVectorList', - 'azimuthTime', - 'datetime') - cls.rg_list_line = \ - cls._parse_vectorlist('noiseVectorList', - 'line', - 'scalar_int') - cls.rg_list_pixel = \ - cls._parse_vectorlist('noiseVectorList', - 'pixel', - 'vector_int') - cls.rg_list_noise_range_lut = \ - cls._parse_vectorlist('noiseVectorList', - 'noiseLut', - 'vector_float') - - cls.az_first_azimuth_line = None - cls.az_first_range_sample = None - cls.az_last_azimuth_line = None - cls.az_last_range_sample = None - cls.az_line = None - cls.az_noise_azimuth_lut = None - - else: - cls.rg_list_azimuth_time = \ - cls._parse_vectorlist('noiseRangeVectorList', - 'azimuthTime', - 'datetime') - cls.rg_list_line = \ - cls._parse_vectorlist('noiseRangeVectorList', - 'line', - 'scalar_int') - cls.rg_list_pixel = \ - cls._parse_vectorlist('noiseRangeVectorList', - 'pixel', - 'vector_int') - cls.rg_list_noise_range_lut = \ - cls._parse_vectorlist('noiseRangeVectorList', - 'noiseRangeLut', - 'vector_float') - cls.az_first_azimuth_line = \ - cls._parse_vectorlist('noiseAzimuthVectorList', - 'firstAzimuthLine', - 'scalar_int')[0] - cls.az_first_range_sample = \ - cls._parse_vectorlist('noiseAzimuthVectorList', - 'firstRangeSample', - 'scalar_int')[0] - cls.az_last_azimuth_line = \ - cls._parse_vectorlist('noiseAzimuthVectorList', - 'lastAzimuthLine', - 'scalar_int')[0] - cls.az_last_range_sample = \ - cls._parse_vectorlist('noiseAzimuthVectorList', - 'lastRangeSample', - 'scalar_int')[0] - cls.az_line = \ - cls._parse_vectorlist('noiseAzimuthVectorList', - 'line', - 'vector_int')[0] - cls.az_noise_azimuth_lut = \ - cls._parse_vectorlist('noiseAzimuthVectorList', - 'noiseAzimuthLut', - 'vector_float')[0] - - return cls - - @dataclass class ProductAnnotation(AnnotationBase): ''' @@ -535,135 +427,6 @@ def closest_block_to_azimuth_time(vector_azimuth_time: np.ndarray, return np.argmin(np.abs(vector_azimuth_time-azimuth_time_burst)) -@dataclass -class BurstNoise: - '''Noise correction information for Sentinel-1 burst''' - basename_nads: str - range_azimith_time: datetime.datetime - range_line: float - range_pixel: np.ndarray - range_lut: np.ndarray - - azimuth_first_azimuth_line: int - azimuth_first_range_sample: int - azimuth_last_azimuth_line: int - azimuth_last_range_sample: int - azimuth_line: np.ndarray - azimuth_lut: np.ndarray - - line_from: int - line_to: int - - @classmethod - def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, - azimuth_time: datetime.datetime, - line_from: int, - line_to: int, - ipf_version: version.Version): - ''' - Extracts the noise correction information for individual burst from NoiseAnnotation - - Parameters - ---------- - noise_annotation: NoiseAnnotation - Subswath-wide noise annotation information - azimuth_time : datetime.datetime - Azimuth time of the burst - line_from: int - First line of the burst in the subswath - line_to: int - Last line of the burst in the subswath - ipf_version: float - IPF version of the SAFE data - - Returns - ------- - cls: BurstNoise - Instance of BurstNoise initialized by the input parameters - - ''' - - basename_nads = noise_annotation.basename_annotation - id_closest = closest_block_to_azimuth_time(noise_annotation.rg_list_azimuth_time, - azimuth_time) - - range_azimith_time = noise_annotation.rg_list_azimuth_time[id_closest] - range_line = noise_annotation.rg_list_line[id_closest] - range_pixel = noise_annotation.rg_list_pixel[id_closest] - range_lut = noise_annotation.rg_list_noise_range_lut[id_closest] - - azimuth_first_azimuth_line = noise_annotation.az_first_azimuth_line - azimuth_first_range_sample = noise_annotation.az_first_range_sample - azimuth_last_azimuth_line = noise_annotation.az_last_azimuth_line - azimuth_last_range_sample = noise_annotation.az_last_range_sample - - if ipf_version >= min_ipf_version_az_noise_vector: - # Azimuth noise LUT exists - crop to the extent of the burst - id_top = np.argmin(np.abs(noise_annotation.az_line-line_from)) - id_bottom = np.argmin(np.abs(noise_annotation.az_line-line_to)) - - # put some margin when possible - if id_top > 0: - id_top -= 1 - if id_bottom < len(noise_annotation.az_line)-1: - id_bottom += 1 - azimuth_line = noise_annotation.az_line[id_top:id_bottom + 1] - azimuth_lut = noise_annotation.az_noise_azimuth_lut[id_top:id_bottom + 1] - - else: - azimuth_line = None - azimuth_lut = None - - return cls(basename_nads, range_azimith_time, range_line, range_pixel, range_lut, - azimuth_first_azimuth_line, azimuth_first_range_sample, - azimuth_last_azimuth_line, azimuth_last_range_sample, - azimuth_line, azimuth_lut, - line_from, line_to) - - - def compute_thermal_noise_lut(self, shape_lut): - ''' - Calculate thermal noise LUT whose shape is `shape_lut` - - Parameter: - ---------- - shape_lut: tuple or list - Shape of the output LUT - - Returns - ------- - arr_lut_total: np.ndarray - 2d array containing thermal noise correction look up table values - ''' - - nrows, ncols = shape_lut - - # Interpolate the range noise vector - rg_lut_interp_obj = InterpolatedUnivariateSpline(self.range_pixel, - self.range_lut, - k=1) - if self.azimuth_last_range_sample is not None: - vec_rg = np.arange(self.azimuth_last_range_sample + 1) - else: - vec_rg = np.arange(ncols) - rg_lut_interpolated = rg_lut_interp_obj(vec_rg) - - # Interpolate the azimuth noise vector - if (self.azimuth_line is None) or (self.azimuth_lut is None): - az_lut_interpolated = np.ones(nrows) - else: # IPF >= 2.90 - az_lut_interp_obj = InterpolatedUnivariateSpline(self.azimuth_line, - self.azimuth_lut, - k=1) - vec_az = np.arange(self.line_from, self.line_to + 1) - az_lut_interpolated = az_lut_interp_obj(vec_az) - - arr_lut_total = np.matmul(az_lut_interpolated[..., np.newaxis], - rg_lut_interpolated[np.newaxis, ...]) - - return arr_lut_total - - @dataclass class BurstCalibration: '''Calibration information for Sentinel-1 IW SLC burst diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 27a3c20..bc209b4 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -638,6 +638,16 @@ def thermal_noise_lut(self): return self.burst_noise.compute_thermal_noise_lut(self.shape) + @property + def thermal_noise_lut(self): + ''' + Returns the LUT for thermal noise correction for the burst + ''' + if self.burst_noise is None: + raise ValueError('burst_noise is not defined for this burst.') + + return self.burst_noise.compute_thermal_noise_lut(self.shape) + @property def eap_compensation_lut(self): '''Returns LUT for EAP compensation. diff --git a/src/s1reader/s1_noise.py b/src/s1reader/s1_noise.py new file mode 100644 index 0000000..75bf389 --- /dev/null +++ b/src/s1reader/s1_noise.py @@ -0,0 +1,156 @@ +from dataclasses import dataclass +import datetime +import xml.etree.ElementTree as ET + +import numpy as np +from packaging import version +from scipy.interpolate import InterpolatedUnivariateSpline + +from s1reader.utils.utils import as_datetime, min_ipf_version_az_noise_vector + +@dataclass +class BurstNoise: #For thermal noise correction + '''Noise correction information for Sentinel-1 burst''' + range_azimith_time: datetime.datetime + range_line: float # TODO is this type correct? + range_pixel: np.ndarray + range_lut: np.ndarray + + azimuth_first_azimuth_line: int + azimuth_first_range_sample: int + azimuth_last_azimuth_line: int + azimuth_last_range_sample: int + azimuth_line: np.ndarray + azimuth_lut: np.ndarray + + line_from: int + line_to: int + + + def compute_thermal_noise_lut(self, shape_lut): + ''' + Calculate thermal noise LUT whose shape is `shape_lut` + + Parameter: + ---------- + shape_lut: tuple or list + Shape of the output LUT + + Returns + ------- + arr_lut_total: np.ndarray + 2d array containing thermal noise correction look up table values + ''' + + nrows, ncols = shape_lut + + # Interpolate the range noise vector + rg_lut_interp_obj = InterpolatedUnivariateSpline(self.range_pixel, + self.range_lut, + k=1) + if self.azimuth_last_range_sample is not None: + vec_rg = np.arange(self.azimuth_last_range_sample + 1) + else: + vec_rg = np.arange(ncols) + rg_lut_interpolated = rg_lut_interp_obj(vec_rg) + + # Interpolate the azimuth noise vector + if (self.azimuth_line is None) or (self.azimuth_lut is None): + az_lut_interpolated = np.ones(nrows) + else: # IPF >= 2.90 + az_lut_interp_obj = InterpolatedUnivariateSpline(self.azimuth_line, + self.azimuth_lut, + k=1) + vec_az = np.arange(self.line_from, self.line_to + 1) + az_lut_interpolated = az_lut_interp_obj(vec_az) + + arr_lut_total = np.matmul(az_lut_interpolated[..., np.newaxis], + rg_lut_interpolated[np.newaxis, ...]) + + return arr_lut_total + + +@dataclass +class BurstNoiseLoader: + ''' document me plz + ''' + azimuth_first_azimuth_line: int + azimuth_first_range_sample: int + azimuth_last_azimuth_line: int + azimuth_last_range_sample: int + azimuth_line: np.ndarray + azimuth_lut: np.ndarray + + noise_rg_vec_list: ET + noise_rg_key: str + + @classmethod + def from_file(cls, noise_path: str, ipf_version: version.Version, + open_method=open): + ''' document me plz + ''' + # TODO comments to explain WTF is going on + with open_method(noise_path, 'r') as f_noise: + noise_tree = ET.parse(f_noise) + + #legacy SAFE data + if ipf_version < min_ipf_version_az_noise_vector: + az_first_azimuth_line = None + az_first_range_sample = None + az_last_azimuth_line = None + az_last_range_sample = None + az_line = None + az_lut = None + + noise_rg_vec_list = noise_tree.find('noiseVectorList') + noise_rg_key = 'noiseLut' + else: + az_noise_tree = noise_tree.find('noiseAzimuthVector') + az_first_azimuth_line = int(az_noise_tree.find('firstAzimuthLine').text) + az_first_range_sample = int(az_noise_tree.find('firstRangeSample').text) + az_last_azimuth_line = int(az_noise_tree.find('lastAzimuthLine').text) + az_last_range_sample = int(az_noise_tree.find('lastRangeSample').text) + az_line = np.array([int(x) + for x in az_noise_tree.find('line').text.split()]) + az_lut = np.array([float(x) + for x in az_noise_tree.find('noiseAzimuthLut').text.split()]) + + noise_rg_vec_list = noise_tree.find('noiseRangeVectorList') + noise_rg_key = 'noiseRangeLut' + + return cls(az_first_azimuth_line, az_first_range_sample, + az_last_azimuth_line, az_last_range_sample, + az_line, az_lut, noise_rg_vec_list, noise_rg_key) + + def get_nearest_noise(self, burst_az_time, line_from, line_to): + ''' document me plz + ''' + # find closest az time + nearest_noise_rg_vec = None + nearest_rg_az_time = None + min_dt = 365 * 24 * 2600 + for noise_rg_vec in self.noise_rg_vec_list: + rg_az_time = as_datetime(noise_rg_vec.find('azimuthTime').text) + dt = abs(rg_az_time - burst_az_time) + if dt < min_dt: + nearest_noise_rg_vec = noise_rg_vec + nearest_rg_az_time = rg_az_time + min_dt = dt + continue + if dt > min_dt: + break + + rg_line = int(nearest_noise_rg_vec.find('line').text) + rg_pixels = np.array([int(x) + for x in nearest_noise_rg_vec.find('pixel').text.split()]) + + rg_lut = np.array([float(x) + for x in nearest_noise_rg_vec.find(self.noise_rg_key).text.split()]) + + return BurstNoise(nearest_rg_az_time, rg_line, rg_pixels, rg_lut, + self.azimuth_first_azimuth_line, + self.azimuth_first_range_sample, + self.azimuth_last_azimuth_line, + self.azimuth_last_range_sample, + self.azimuth_line, self.azimuth_lut, + line_from, line_to) diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index df641bd..4a13f78 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -15,34 +15,14 @@ from nisar.workflows.stage_dem import check_dateline from s1reader import s1_annotation # to access __file__ -from s1reader.s1_annotation import ProductAnnotation, NoiseAnnotation,\ - CalibrationAnnotation, AuxCal, \ - BurstCalibration, BurstEAP, BurstNoise - +from s1reader.s1_annotation import (ProductAnnotation, CalibrationAnnotation, + AuxCal, BurstCalibration, BurstEAP) from s1reader.s1_burst_slc import Doppler, Sentinel1BurstSlc - +from s1reader.s1_noise import BurstNoise, BurstNoiseLoader +from s1reader.utils.utils import as_datetime esa_track_burst_id_file = f"{os.path.dirname(os.path.realpath(__file__))}/data/sentinel1_track_burst_id.txt" -# TODO evaluate if it make sense to combine below into a class -def as_datetime(t_str, fmt = "%Y-%m-%dT%H:%M:%S.%f"): - '''Parse given time string to datetime.datetime object. - - Parameters: - ---------- - 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: - ------ - _ : datetime.datetime - datetime.datetime object parsed from given time string. - ''' - return datetime.datetime.strptime(t_str, fmt) - def parse_polynomial_element(elem, poly_name): '''Parse azimuth FM (Frequency Modulation) rate element to reference time and poly1d tuples. @@ -376,7 +356,6 @@ def get_path_aux_cal(directory_aux_cal: str, str_annotation: str): return list_aux_cal[id_match] - def is_eap_correction_necessary(ipf_version: version.Version) -> SimpleNamespace : ''' Examines if what level of elevation antenna pattern (EAP) correction is necessary. @@ -495,10 +474,8 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, # load the Noise annotation noise_annotation_path = annotation_path.replace('annotation/', 'annotation/calibration/noise-') - with open_method(noise_annotation_path, 'r') as f_nads: - tree_nads = ET.parse(f_nads) - noise_annotation = NoiseAnnotation.from_et(tree_nads, ipf_version, - noise_annotation_path) + noise_loader = BurstNoiseLoader.from_file(noise_annotation_path, + ipf_version) # load AUX_CAL annotation eap_necessity = is_eap_correction_necessary(ipf_version) @@ -638,8 +615,11 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, # Extract burst-wise information for Calibration, Noise, and EAP correction burst_calibration = BurstCalibration.from_calibration_annotation(calibration_annotation, sensing_start) - burst_noise = BurstNoise.from_noise_annotation(noise_annotation, sensing_start, - i*n_lines, (i+1)*n_lines-1, ipf_version) + + burst_noise = noise_loader.get_nearest_noise(sensing_start, + i * n_lines, + (i + 1) * n_lines - 1) + if aux_cal_subswath is None: # Not applying EAP correction; (IPF high enough or user turned that off) # No need to fill in `burst_aux_cal` diff --git a/src/s1reader/utils/__init__.py b/src/s1reader/utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/s1reader/utils/utils.py b/src/s1reader/utils/utils.py new file mode 100644 index 0000000..63160c6 --- /dev/null +++ b/src/s1reader/utils/utils.py @@ -0,0 +1,25 @@ +import datetime + +from packaging import version + +# Minimum IPF version from which the S1 product's Noise Annotation +# Data Set (NADS) includes azimuth noise vector annotation +min_ipf_version_az_noise_vector = version.parse('2.90') + +def as_datetime(t_str, fmt = "%Y-%m-%dT%H:%M:%S.%f"): + '''Parse given time string to datetime.datetime object. + + Parameters: + ---------- + 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: + ------ + _ : datetime.datetime + datetime.datetime object parsed from given time string. + ''' + return datetime.datetime.strptime(t_str, fmt)