From dffe813f7eb6b7d8cfa71ca84016b7b95da705e5 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Tue, 9 Aug 2022 22:57:00 -0700 Subject: [PATCH 01/49] Replacing the loaders in Burst* class into class methods, with further implementation for thermal and EAP correction --- src/s1reader/s1_annotation.py | 228 ++++++++++++++++++++++++++-------- src/s1reader/s1_reader.py | 5 +- 2 files changed, 176 insertions(+), 57 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 5d2f7d1e..0f78ed1a 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -10,9 +10,9 @@ import numpy as np from packaging import version -from scipy.interpolate import InterpolatedUnivariateSpline +from scipy.interpolate import InterpolatedUnivariateSpline, interp1d -#some thresholds +#A IPF version from which has azimuth noise vector version_threshold_azimuth_noise_vector=version.parse('2.90') @dataclass @@ -152,8 +152,6 @@ def from_et(cls, et_in=None): return cls - - @dataclass class NoiseAnnotation(AnnotationBase): '''Reader for Noise Annotation Data Set (NADS) for IW SLC''' @@ -337,7 +335,7 @@ def closest_block_to_azimuth_time(vector_azimuth_time: np.ndarray, Returns ------- - int + _: int Index of vector_azimuth_time that is the closest to azimuth_burst_time ''' @@ -363,8 +361,8 @@ class BurstNoise: #For thermal noise correction line_from: int = None line_to: int = None - - def from_noise_annotation(self, noise_annotation: NoiseAnnotation, azimuth_time: datetime.datetime, + @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 @@ -381,21 +379,27 @@ def from_noise_annotation(self, noise_annotation: NoiseAnnotation, azimuth_time: ipf_version: float IPF version of the SAFE data + Returns + ------- + cls: BurstNoise + An instance from BurstNoise initialized by the input parameters + ''' id_closest = closest_block_to_azimuth_time(noise_annotation.rg_list_azimuth_time, azimuth_time) - self.range_azimith_time = noise_annotation.rg_list_azimuth_time[id_closest] - self.range_line = noise_annotation.rg_list_line[id_closest] - self.range_pixel = noise_annotation.rg_list_pixel[id_closest] - self.range_lut = noise_annotation.rg_list_noise_range_lut[id_closest] - self.azimuth_first_azimuth_line = noise_annotation.az_first_azimuth_line - self.azimuth_first_range_sample = noise_annotation.az_first_range_sample - self.azimuth_last_azimuth_line = noise_annotation.az_last_azimuth_line - self.azimuth_last_range_sample = noise_annotation.az_last_range_sample + 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] - self.line_from = line_from - self.line_to = line_to + 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 + + line_from = line_from + line_to = line_to if ipf_version >= version_threshold_azimuth_noise_vector: #Azimuth noise LUT exists - crop to the extent of the burst @@ -406,28 +410,45 @@ def from_noise_annotation(self, noise_annotation: NoiseAnnotation, azimuth_time: id_top -= 1 if id_bottom < len(noise_annotation.az_line)-1: id_bottom += 1 - self.azimuth_line = noise_annotation.az_line[id_top:id_bottom + 1] - self.azimuth_lut = noise_annotation.az_noise_azimuth_lut[id_top: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(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 export_lut(self): - '''Gives out the LUT table whose size is the same as the burst SLC''' - ncols = self.azimuth_last_range_sample-self.azimuth_first_range_sample+1 - nrows = self.line_to-self.line_from+1 + '''Returns the LUT table for thermal noise correction whose size is the same as the burst SLC''' + ncols = self.azimuth_last_range_sample - self.azimuth_first_range_sample + 1 + nrows = self.line_to - self.line_from + 1 + #interpolator for range noise vector intp_rg_lut = InterpolatedUnivariateSpline(self.range_pixel, self.range_lut, k=1) - intp_az_lut = InterpolatedUnivariateSpline(self.azimuth_line, self.azimuth_lut, k=1) + grid_rg = np.arange(self.azimuth_last_range_sample + 1) + rg_lut_interp = intp_rg_lut(grid_rg).reshape((1, ncols)) - grid_rg = np.arange(self.azimuth_last_range_sample+1) - grid_az = np.arange(self.line_from, self.line_to+1) + #interpolator for azimuth noise vector - take IPF version into consideration + if (self.azimuth_line is None) or (self.azimuth_lut is None): # IPF <2.90 + az_lut_interp = np.ones(nrows).reshape((nrows, 1)) - rg_lut_interp = intp_rg_lut(grid_rg).reshape((1,ncols)) - az_lut_interp = intp_az_lut(grid_az).reshape((nrows,1)) + else: #IPF >= 2.90 + intp_az_lut = InterpolatedUnivariateSpline(self.azimuth_line, self.azimuth_lut, k=1) + grid_az = np.arange(self.line_from,self.line_to + 1) + az_lut_interp = intp_az_lut(grid_az).reshape((nrows, 1)) arr_lut_total = np.matmul(az_lut_interp, rg_lut_interp) + return arr_lut_total + @dataclass class BurstCalibration: '''Calibration information for Sentinel-1 IW SLC burst @@ -458,36 +479,38 @@ def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotati EAP correction info for the burst ''' id_closest = closest_block_to_azimuth_time(calibration_annotation.list_azimuth_time, azimuth_time) - cls.azimuth_time = calibration_annotation.list_azimuth_time[id_closest] - cls.line = calibration_annotation.list_line[id_closest] - cls.pixel = calibration_annotation.list_pixel[id_closest] - cls.sigma_naught = calibration_annotation.list_sigma_nought[id_closest] - cls.beta_naught = calibration_annotation.list_beta_nought[id_closest] - cls.gamma = calibration_annotation.list_gamma[id_closest] - cls.dn = calibration_annotation.list_dn[id_closest] + azimuth_time = calibration_annotation.list_azimuth_time[id_closest] + line = calibration_annotation.list_line[id_closest] + pixel = calibration_annotation.list_pixel[id_closest] + sigma_naught = calibration_annotation.list_sigma_nought[id_closest] + beta_naught = calibration_annotation.list_beta_nought[id_closest] + gamma = calibration_annotation.list_gamma[id_closest] + dn = calibration_annotation.list_dn[id_closest] - return cls + return cls(azimuth_time, line, pixel, + sigma_naught, beta_naught, gamma, dn) @dataclass class BurstEAP: - '''Elevation antenna pattern (EAP) correction information for Sentinel-1 IW SLC burst + '''EAP correction information for Sentinel-1 IW SLC burst ''' #from LADS - Ns: int #number of samples - fs: float #range sampling rate + Ns: int # number of samples + fs: float # range sampling rate eta_start: datetime.datetime - tau_0: float #imageInformation/slantRangeTime - tau_sub: np.ndarray #antennaPattern/slantRangeTime - theta_sub: np.ndarray #antennaPattern/elevationAngle + tau_0: float # imageInformation/slantRangeTime + tau_sub: np.ndarray # antennaPattern/slantRangeTime + theta_sub: np.ndarray # antennaPattern/elevationAngle azimuth_time: datetime.datetime + ascending_node_time: datetime.datetime #from AUX_CAL - G_eap: np.ndarray #elevationAntennaPattern - delta_theta:float #elavationAngleIncrement + G_eap: np.ndarray # elevationAntennaPattern + delta_theta:float # elavationAngleIncrement @classmethod - def from_product_annotation_and_aux_cal(cls, product_annotation: ProductAnnotation, aux_cal: AuxCal, azimuth_time: datetime.datetime): + def from_product_annotation_and_aux_cal(cls,product_annotation: ProductAnnotation, aux_cal: AuxCal, azimuth_time: datetime.datetime): ''' A class method that extracts the EAP correction info for the IW SLC burst @@ -495,25 +518,122 @@ def from_product_annotation_and_aux_cal(cls, product_annotation: ProductAnnotati ---------- product_annotation: ProductAnnotation A swath-wide product annotation class + aux_cal: AuxCal AUX_CAL information that corresponds to the sensing time + azimuth_time: datetime.datetime Azimuth time of the burst Returns ------- cls: BurstEAP - EAP correction info for the burst + A burst-wide information for EAP correction + ''' id_closest = closest_block_to_azimuth_time(product_annotation.antenna_pattern_azimuth_time, azimuth_time) - cls.Ns = product_annotation.number_of_samples - cls.fs = product_annotation.range_sampling_rate - cls.eta_start = azimuth_time - cls.tau_0 = product_annotation.antenna_pattern_slant_range_time[id_closest] - cls.tau_sub = product_annotation.antenna_pattern_slant_range_time[id_closest] - cls.theta_am = product_annotation.antenna_pattern_elevation_angle + Ns = product_annotation.number_of_samples + fs = product_annotation.range_sampling_rate + eta_start = azimuth_time + tau_0 = product_annotation.antenna_pattern_slant_range_time[id_closest] + tau_sub = product_annotation.antenna_pattern_slant_range_time[id_closest] + theta_sub = product_annotation.antenna_pattern_elevation_pattern[id_closest] + #self.theta_am = product_annotation.antenna_pattern_elevation_angle + G_eap = aux_cal.elevation_antenna_pattern + delta_theta = aux_cal.elevation_angle_increment + + ascending_node_time = product_annotation.ascending_node_time + + return cls(Ns, fs, eta_start, tau_0, tau_sub, theta_sub, + azimuth_time, ascending_node_time, + G_eap, delta_theta) + + + def export_lut(self): + '''Returns LUT for EAP correction. Parts of the codes were borrowed from ISCE2. + Based on ESA docuemnt "Impact of the Elevation Antenna Pattern Phase Compensation on the Interferometric Phase Preservation"''' + + #Step 1. Retrieve two-way complex EAP term + n_elt = len(self.G_eap) + + theta_am = (np.arange(n_elt)-(n_elt-1)/2)* self.delta_theta + + #Step 2. finding the roll steering angle + #2.1. get ascending node time - DONE + delta_anx = self.eta_start-self.ascending_node_time + theta_offnadir = self._anx2roll(delta_anx) + + #2.2. get state vector to calculate satellite height - Is it necessary? + + + #Step 3. Computing the elevtion angle in the geometry of view + theta_eap = theta_am+theta_offnadir + + #Step 4. re-interpolating the 2-way complex EAP + tau = self.tau_0+np.arange(self.Ns)/self.fs + + #4.1. set up interpolator + theta = np.interp(tau, self.tau_sub, self.theta_sub) + + interpolator_G = interp1d(theta_eap,self.G_eap) + G_eap_interpolated = interpolator_G(theta) + phi_EAP = np.angle(G_eap_interpolated) + cJ = np.complex64(1.0j) + G_EAP = np.exp(cJ * phi_EAP) + return G_EAP + + def _anx2roll(self, delta_anx): + ''' + Returns the Platform nominal roll as function of elapsed time from + ascending node crossing time (ANX). + Straight from S1A documentation. + ''' + + ####Estimate altitude based on time elapsed since ANX + altitude = self._anx2height(delta_anx) + + ####Reference altitude + href = 711.700 #;km + + ####Reference boresight at reference altitude + boresight_ref = 29.450 # ; deg + + ####Partial derivative of roll vs altitude + alpha_roll = 0.0566 # ;deg/km + + ####Estimate nominal roll + nominal_roll = boresight_ref - alpha_roll * (altitude/1000.0 - href) #Theta off nadir + + return nominal_roll + + + def _anx2height(self, delta_anx): + ''' + Returns the platform nominal height as function of elapse time from + ascending node crossing time (ANX). + Straight from S1A documention. + ''' + + ###Average height + h0 = 707714.8 #;m + + ####Perturbation amplitudes + h = np.array([8351.5, 8947.0, 23.32, 11.74]) #;m + + ####Perturbation phases + phi = np.array([3.1495, -1.5655 , -3.1297, 4.7222]) #;radians + + ###Orbital time period in seconds + Torb = (12*24*60*60)/175. + + ###Angular velocity + worb = 2*np.pi / Torb + + ####Evaluation of series + ht=h0 + for i in range(len(h)): + ht += h[i] * np.sin((i+1) * worb * delta_anx + phi[i]) + + return ht - cls.G_eap = aux_cal.elevation_antenna_pattern - cls.delta_theta = aux_cal.elevation_angle_increment - return cls diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index 4317ae6a..d43eb43c 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -492,9 +492,8 @@ 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 = s1_annotation.BurstCalibration.from_calibration_annotation(calibration_annotation, sensing_start) - bursts_noise=s1_annotation.BurstNoise() - bursts_noise.from_noise_annotation(noise_annotation,sensing_start,i*n_lines,(i+1)*n_lines-1,ipf_version) - + bursts_noise=s1_annotation.BurstNoise.from_noise_annotation(noise_annotation, sensing_start, i*n_lines, (i+1)*n_lines-1, ipf_version) + bursts[i] = Sentinel1BurstSlc(ipf_version, sensing_start, radar_freq, wavelength, azimuth_steer_rate, azimuth_time_interval, slant_range_time, starting_range, iw2_mid_range, From 0ec4fc4b78a688818ca8ef66fca6317881ad6a28 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 10 Aug 2022 15:30:10 -0700 Subject: [PATCH 02/49] fix on determining beta_naught; addressing PEP8 issues --- src/s1reader/s1_annotation.py | 128 ++++++++++++++++++---------------- src/s1reader/s1_reader.py | 14 ++-- 2 files changed, 78 insertions(+), 64 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 0f78ed1a..90a81054 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -31,8 +31,10 @@ def _parse_scalar(cls, path_field: str, str_type: str): ---------- path_field : str Field in the xml_et to parse - str_type : str {'datetime', 'scalar_int', 'scalar_float', 'vector_int', 'vector_float', 'str'} - Specify how the texts in the field will be parsed + str_type : str + Specify how the texts in the field will be parsed. + accepted values: + {'datetime', 'scalar_int', 'scalar_float', 'vector_int', 'vector_float', 'str'} Returns ------- @@ -75,9 +77,12 @@ def _parse_vectorlist(cls, name_vector_list: str, name_vector: str, str_type: st name_vector_list : str List Field in the xml_et to parse name_vector : str - Name of the field in each elements of the VectorList (e.g. 'noiseLut' in 'noiseVectorList') - str_type : str {'datetime', 'scalar_int', 'scalar_float', 'vector_int', 'vector_float', 'str'} - Specify how the texts in the field will be parsed: + Name of the field in each elements of the VectorList + (e.g. 'noiseLut' in 'noiseVectorList') + str_type : str + Specify how the texts in the field will be parsed + accepted values: + {'datetime', 'scalar_int', 'scalar_float', 'vector_int', 'vector_float', 'str'} Returns ------- @@ -138,8 +143,17 @@ class CalibrationAnnotation(AnnotationBase): list_dn: list @classmethod - def from_et(cls, et_in=None): - '''Extracts the list of calibration informaton from etree from CADS''' + def from_et(cls, et_in: ET): + '''Extracts the list of calibration informaton from etree from CADS + Parameters: + ----------- + et_in: ElementTree From CADS .xml file + + Returns: + cls: CalibrationAnnotation + instance from CalibrationAnnotation initialized by the input parameter + ''' + cls.xml_et = et_in cls.list_azimuth_time = cls._parse_vectorlist('calibrationVectorList', 'azimuthTime', 'datetime') cls.list_line = cls._parse_vectorlist('calibrationVectorList', 'line', 'scalar_int') @@ -154,13 +168,10 @@ def from_et(cls, et_in=None): @dataclass class NoiseAnnotation(AnnotationBase): - '''Reader for Noise Annotation Data Set (NADS) for IW SLC''' - # NOTE Schema of the NADS is slightly different before/after IPF version 2.90. Needs to be adaptive in accordance with the version. - # The issue above was fixed in further implement of thermal noise correction. - # A separete PR regarding this will be submitted upon the acceptance of this code. - # in ISCE2 code: if float(self.IPFversion) < 2.90: + '''Reader for Noise Annotation Data Set (NADS) for IW SLC # REF: .../isce2/components/isceobj/Sensor/GRD/Sentinel1.py - + # : ESA Documentation "Thermal Denoising of Products Generated by the S-1 IPF" + ''' rg_list_azimuth_time: np.ndarray rg_list_line: list rg_list_pixel: list @@ -308,14 +319,13 @@ def from_et(cls,et_in: ET, pol: str, str_swath: str): if n_val == len(arr_eap_val): #Provided in real numbers: In case of AUX_CAL for old IPFs. cls.azimuth_antenna_element_pattern=arr_eap_val elif n_val*2 == len(arr_eap_val): #Provided in complex numbers: In case of recent IPFs e.g. 3.10 - cls.azimuth_antenna_element_pattern=arr_eap_val[0::2]+arr_eap_val[1::2]*1.0j + cls.azimuth_antenna_element_pattern = arr_eap_val[0::2]+arr_eap_val[1::2]*1.0j else: raise ValueError(f'The number of values does not match. n_val={n_val}, #values in elevationAntennaPattern/values={len(arr_eap_val)}') cls.azimuth_angle_increment = float(calibration_params.find('azimuthAntennaPattern/azimuthAngleIncrement').text) - cls.azimuth_antenna_pattern = np.array( - [float(token_val) for token_val in calibration_params.find('azimuthAntennaPattern/values').text.split()] - ) + cls.azimuth_antenna_pattern = \ + np.array([float(token_val) for token_val in calibration_params.find('azimuthAntennaPattern/values').text.split()]) cls.absolute_calibration_constant = float(calibration_params.find('absoluteCalibrationConstant').text) cls.noise_calibration_factor = float(calibration_params.find('noiseCalibrationFactor').text) @@ -324,7 +334,8 @@ def from_et(cls,et_in: ET, pol: str, str_swath: str): def closest_block_to_azimuth_time(vector_azimuth_time: np.ndarray, azimuth_time_burst: datetime.datetime) -> int: - '''Find the id of the closest data block in annotation. To be used when populating BurstNoise, BurstCalibration, and BurstEAP. + '''Find the id of the closest data block in annotation. + To be used when populating BurstNoise, BurstCalibration, and BurstEAP. Parameters ---------- @@ -346,20 +357,20 @@ def closest_block_to_azimuth_time(vector_azimuth_time: np.ndarray, @dataclass class BurstNoise: #For thermal noise correction '''Noise correction information for Sentinel-1 burst''' - range_azimith_time: datetime.datetime = None - range_line: float = None - range_pixel: np.ndarray = None - range_lut: np.ndarray = None + range_azimith_time: datetime.datetime + range_line: float + range_pixel: np.ndarray + range_lut: np.ndarray - azimuth_first_azimuth_line: int = None - azimuth_first_range_sample: int = None - azimuth_last_azimuth_line: int = None - azimuth_last_range_sample: int = None - azimuth_line: np.ndarray = None - azimuth_lut: np.ndarray = None + 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 = None - line_to: int = None + line_from: int + line_to: int @classmethod def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: datetime.datetime, @@ -398,9 +409,6 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: azimuth_last_azimuth_line = noise_annotation.az_last_azimuth_line azimuth_last_range_sample = noise_annotation.az_last_range_sample - line_from = line_from - line_to = line_to - if ipf_version >= version_threshold_azimuth_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)) @@ -425,7 +433,7 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: def export_lut(self): - '''Returns the LUT table for thermal noise correction whose size is the same as the burst SLC''' + '''Returns the burst-sized LUT for thermal noise correction''' ncols = self.azimuth_last_range_sample - self.azimuth_first_range_sample + 1 nrows = self.line_to - self.line_from + 1 @@ -434,11 +442,10 @@ def export_lut(self): grid_rg = np.arange(self.azimuth_last_range_sample + 1) rg_lut_interp = intp_rg_lut(grid_rg).reshape((1, ncols)) - #interpolator for azimuth noise vector - take IPF version into consideration - if (self.azimuth_line is None) or (self.azimuth_lut is None): # IPF <2.90 + #interpolator for azimuth noise vector + if (self.azimuth_line is None) or (self.azimuth_lut is None): # IPF <2.90 az_lut_interp = np.ones(nrows).reshape((nrows, 1)) - - else: #IPF >= 2.90 + else: # IPF >= 2.90 intp_az_lut = InterpolatedUnivariateSpline(self.azimuth_line, self.azimuth_lut, k=1) grid_az = np.arange(self.line_from,self.line_to + 1) az_lut_interp = intp_az_lut(grid_az).reshape((nrows, 1)) @@ -462,7 +469,8 @@ class BurstCalibration: dn: np.ndarray = None @classmethod - def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotation, azimuth_time: datetime.datetime): + def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotation, + azimuth_time: datetime.datetime): ''' A class method that extracts the calibration info for the burst @@ -478,15 +486,24 @@ def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotati cls: BurstCalibration EAP correction info for the burst ''' + id_closest = closest_block_to_azimuth_time(calibration_annotation.list_azimuth_time, azimuth_time) azimuth_time = calibration_annotation.list_azimuth_time[id_closest] line = calibration_annotation.list_line[id_closest] pixel = calibration_annotation.list_pixel[id_closest] sigma_naught = calibration_annotation.list_sigma_nought[id_closest] - beta_naught = calibration_annotation.list_beta_nought[id_closest] + #beta_naught = calibration_annotation.list_beta_nought[id_closest] gamma = calibration_annotation.list_gamma[id_closest] dn = calibration_annotation.list_dn[id_closest] + #Check if all values in the list of beta_naught LUTs are the same + matrix_beta_naught = np.array(calibration_annotation.list_beta_nought) + if matrix_beta_naught.min() == matrix_beta_naught.max(): + beta_naught = np.min(matrix_beta_naught) + else: + #TODO Switch to LUT-based method when there is significant changes in the array + beta_naught = np.mean(matrix_beta_naught) + return cls(azimuth_time, line, pixel, sigma_naught, beta_naught, gamma, dn) @@ -510,7 +527,8 @@ class BurstEAP: delta_theta:float # elavationAngleIncrement @classmethod - def from_product_annotation_and_aux_cal(cls,product_annotation: ProductAnnotation, aux_cal: AuxCal, azimuth_time: datetime.datetime): + def from_product_annotation_and_aux_cal(cls, product_annotation: ProductAnnotation, + aux_cal: AuxCal, azimuth_time: datetime.datetime): ''' A class method that extracts the EAP correction info for the IW SLC burst @@ -550,29 +568,21 @@ def from_product_annotation_and_aux_cal(cls,product_annotation: ProductAnnotatio def export_lut(self): - '''Returns LUT for EAP correction. Parts of the codes were borrowed from ISCE2. - Based on ESA docuemnt "Impact of the Elevation Antenna Pattern Phase Compensation on the Interferometric Phase Preservation"''' + '''Returns LUT for EAP correction. + Parts of the codes and subfunctions were borrowed from ISCE2. + Based on ESA docuemnt : "Impact of the Elevation Antenna Pattern Phase Compensation on the Interferometric Phase Preservation"''' - #Step 1. Retrieve two-way complex EAP term n_elt = len(self.G_eap) theta_am = (np.arange(n_elt)-(n_elt-1)/2)* self.delta_theta - #Step 2. finding the roll steering angle - #2.1. get ascending node time - DONE delta_anx = self.eta_start-self.ascending_node_time theta_offnadir = self._anx2roll(delta_anx) - #2.2. get state vector to calculate satellite height - Is it necessary? - - - #Step 3. Computing the elevtion angle in the geometry of view theta_eap = theta_am+theta_offnadir - #Step 4. re-interpolating the 2-way complex EAP tau = self.tau_0+np.arange(self.Ns)/self.fs - #4.1. set up interpolator theta = np.interp(tau, self.tau_sub, self.theta_sub) interpolator_G = interp1d(theta_eap,self.G_eap) @@ -580,6 +590,7 @@ def export_lut(self): phi_EAP = np.angle(G_eap_interpolated) cJ = np.complex64(1.0j) G_EAP = np.exp(cJ * phi_EAP) + return G_EAP def _anx2roll(self, delta_anx): @@ -615,7 +626,7 @@ def _anx2height(self, delta_anx): ''' ###Average height - h0 = 707714.8 #;m + h_0 = 707714.8 #;m ####Perturbation amplitudes h = np.array([8351.5, 8947.0, 23.32, 11.74]) #;m @@ -630,10 +641,9 @@ def _anx2height(self, delta_anx): worb = 2*np.pi / Torb ####Evaluation of series - ht=h0 - for i in range(len(h)): - ht += h[i] * np.sin((i+1) * worb * delta_anx + phi[i]) - - return ht - + h_t = h_0 + #for i in range(len(h)): + for i, h_i in enumerate(h): + h_t += h_i * np.sin((i+1) * worb * delta_anx + phi[i]) + return h_t diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index d43eb43c..83dd8abd 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -4,8 +4,8 @@ import xml.etree.ElementTree as ET import zipfile -from packaging import version from types import SimpleNamespace +from packaging import version import isce3 import numpy as np @@ -361,7 +361,8 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, # product_annotation = s1_annotation.ProductAnnotation.from_et(tree_lads) #load the Calibraton annotation - calibration_annotation_path = annotation_path.replace('annotation/', 'annotation/calibration/calibration-') + calibration_annotation_path =\ + annotation_path.replace('annotation/', 'annotation/calibration/calibration-') with open_method(calibration_annotation_path, 'r') as f_cads: tree_cads = ET.parse(f_cads) calibration_annotation = s1_annotation.CalibrationAnnotation.from_et(tree_cads) @@ -491,9 +492,12 @@ 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 = s1_annotation.BurstCalibration.from_calibration_annotation(calibration_annotation, sensing_start) - bursts_noise=s1_annotation.BurstNoise.from_noise_annotation(noise_annotation, sensing_start, i*n_lines, (i+1)*n_lines-1, ipf_version) - + burst_calibration =\ + s1_annotation.BurstCalibration.from_calibration_annotation(calibration_annotation, sensing_start) + bursts_noise=\ + s1_annotation.BurstNoise.from_noise_annotation(noise_annotation, sensing_start, + i*n_lines, (i+1)*n_lines-1, ipf_version) + bursts[i] = Sentinel1BurstSlc(ipf_version, sensing_start, radar_freq, wavelength, azimuth_steer_rate, azimuth_time_interval, slant_range_time, starting_range, iw2_mid_range, From b4d3c0e9dcd0507d3f84a56c3dfe7634c24dc69b Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Fri, 19 Aug 2022 17:48:39 -0700 Subject: [PATCH 03/49] Bug fix and feature addition to BurstEAP; restructuring LUT exportation --- src/s1reader/s1_annotation.py | 108 ++++++++++++++-------------------- src/s1reader/s1_burst_slc.py | 60 +++++++++++++++++++ src/s1reader/s1_reader.py | 72 +++++++++++++++++++---- 3 files changed, 165 insertions(+), 75 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 90a81054..4286a619 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -13,7 +13,7 @@ from scipy.interpolate import InterpolatedUnivariateSpline, interp1d #A IPF version from which has azimuth noise vector -version_threshold_azimuth_noise_vector=version.parse('2.90') +version_threshold_azimuth_noise_vector = version.parse('2.90') @dataclass class AnnotationBase: @@ -147,10 +147,12 @@ def from_et(cls, et_in: ET): '''Extracts the list of calibration informaton from etree from CADS Parameters: ----------- - et_in: ElementTree From CADS .xml file + et_in: ET + ElementTree From CADS .xml file Returns: - cls: CalibrationAnnotation + -------- + cls: CalibrationAnnotation instance from CalibrationAnnotation initialized by the input parameter ''' @@ -232,6 +234,9 @@ class ProductAnnotation(AnnotationBase): '''For L1 SLC product annotation file. For EAP correction.''' image_information_slant_range_time: float + # Attributes to be used when determining what AUX_CAL to load + inst_config_id: int + #elevation_angle: antenna_pattern_azimuth_time: list antenna_pattern_slant_range_time: list @@ -242,6 +247,8 @@ class ProductAnnotation(AnnotationBase): number_of_samples: int range_sampling_rate: float + slant_range_time: float + @classmethod def from_et(cls, et_in: ET): '''Extracts list of product information from etree from L1 annotation data set (LADS) @@ -259,6 +266,7 @@ def from_et(cls, et_in: ET): cls.xml_et = et_in cls.image_information_slant_range_time = cls._parse_scalar('imageAnnotation/imageInformation/slantRangeTime', 'scalar_float') + cls.inst_config_id = cls._parse_scalar('generalAnnotation/downlinkInformationList/downlinkInformation/downlinkValues/instrumentConfigId', 'scalar_int') cls.antenna_pattern_azimuth_time = cls._parse_vectorlist('antennaPattern/antennaPatternList', 'azimuthTime', 'datetime') cls.antenna_pattern_slant_range_time = cls._parse_vectorlist('antennaPattern/antennaPatternList', 'slantRangeTime', 'vector_float') cls.antenna_pattern_elevation_angle = cls._parse_vectorlist('antennaPattern/antennaPatternList', 'elevationAngle', 'vector_float') @@ -268,6 +276,8 @@ def from_et(cls, et_in: ET): cls.number_of_samples = cls._parse_scalar('imageAnnotation/imageInformation/numberOfSamples', 'scalar_int') cls.range_sampling_rate = cls._parse_scalar('generalAnnotation/productInformation/rangeSamplingRate', 'scalar_float') + cls.slant_range_time = cls._parse_scalar('imageAnnotation/imageInformation/slantRangeTime','scalar_float') + return cls @@ -276,8 +286,10 @@ class AuxCal(AnnotationBase): '''AUX_CAL information for elevation antenna pattern (EAP) correction''' beam_nominal_near_range: float beam_nominal_far_range: float + elevation_angle_increment: float elevation_antenna_pattern: np.ndarray + azimuth_angle_increment: float azimuth_antenna_pattern: np.ndarray azimuth_antenna_element_pattern_increment: float @@ -285,8 +297,9 @@ class AuxCal(AnnotationBase): absolute_calibration_constant: float noise_calibration_factor: float + @classmethod - def from_et(cls,et_in: ET, pol: str, str_swath: str): + def from_et(cls, et_in: ET, pol: str, str_swath: str): '''A class method that Extracts list of information AUX_CAL from the input ET. Parameters @@ -317,17 +330,26 @@ def from_et(cls,et_in: ET, pol: str, str_swath: str): n_val = int(calibration_params.find('elevationAntennaPattern/values').attrib['count']) arr_eap_val = np.array([float(token_val) for token_val in calibration_params.find('elevationAntennaPattern/values').text.split()]) if n_val == len(arr_eap_val): #Provided in real numbers: In case of AUX_CAL for old IPFs. - cls.azimuth_antenna_element_pattern=arr_eap_val + cls.elevation_antenna_pattern = arr_eap_val elif n_val*2 == len(arr_eap_val): #Provided in complex numbers: In case of recent IPFs e.g. 3.10 - cls.azimuth_antenna_element_pattern = arr_eap_val[0::2]+arr_eap_val[1::2]*1.0j + cls.elevation_antenna_pattern = arr_eap_val[0::2]+arr_eap_val[1::2]*1.0j else: raise ValueError(f'The number of values does not match. n_val={n_val}, #values in elevationAntennaPattern/values={len(arr_eap_val)}') - cls.azimuth_angle_increment = float(calibration_params.find('azimuthAntennaPattern/azimuthAngleIncrement').text) + cls.azimuth_angle_increment = \ + float(calibration_params.find('azimuthAntennaPattern/azimuthAngleIncrement').text) cls.azimuth_antenna_pattern = \ np.array([float(token_val) for token_val in calibration_params.find('azimuthAntennaPattern/values').text.split()]) - cls.absolute_calibration_constant = float(calibration_params.find('absoluteCalibrationConstant').text) - cls.noise_calibration_factor = float(calibration_params.find('noiseCalibrationFactor').text) + + cls.azimuth_antenna_element_pattern_increment = \ + float(calibration_params.find('azimuthAntennaElementPattern/azimuthAngleIncrement').text) + cls.azimuth_antenna_element_pattern = \ + np.array([float(token_val) for token_val in calibration_params.find('azimuthAntennaElementPattern/values').text.split()]) + + cls.absolute_calibration_constant = \ + float(calibration_params.find('absoluteCalibrationConstant').text) + cls.noise_calibration_factor = \ + float(calibration_params.find('noiseCalibrationFactor').text) return cls @@ -432,30 +454,6 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: line_from, line_to) - def export_lut(self): - '''Returns the burst-sized LUT for thermal noise correction''' - ncols = self.azimuth_last_range_sample - self.azimuth_first_range_sample + 1 - nrows = self.line_to - self.line_from + 1 - - #interpolator for range noise vector - intp_rg_lut = InterpolatedUnivariateSpline(self.range_pixel, self.range_lut, k=1) - grid_rg = np.arange(self.azimuth_last_range_sample + 1) - rg_lut_interp = intp_rg_lut(grid_rg).reshape((1, ncols)) - - #interpolator for azimuth noise vector - if (self.azimuth_line is None) or (self.azimuth_lut is None): # IPF <2.90 - az_lut_interp = np.ones(nrows).reshape((nrows, 1)) - else: # IPF >= 2.90 - intp_az_lut = InterpolatedUnivariateSpline(self.azimuth_line, self.azimuth_lut, k=1) - grid_az = np.arange(self.line_from,self.line_to + 1) - az_lut_interp = intp_az_lut(grid_az).reshape((nrows, 1)) - - arr_lut_total = np.matmul(az_lut_interp, rg_lut_interp) - - return arr_lut_total - - - @dataclass class BurstCalibration: '''Calibration information for Sentinel-1 IW SLC burst @@ -487,7 +485,9 @@ def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotati EAP correction info for the burst ''' - id_closest = closest_block_to_azimuth_time(calibration_annotation.list_azimuth_time, azimuth_time) + id_closest = closest_block_to_azimuth_time(calibration_annotation.list_azimuth_time, + azimuth_time) + azimuth_time = calibration_annotation.list_azimuth_time[id_closest] line = calibration_annotation.list_line[id_closest] pixel = calibration_annotation.list_pixel[id_closest] @@ -501,7 +501,7 @@ def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotati if matrix_beta_naught.min() == matrix_beta_naught.max(): beta_naught = np.min(matrix_beta_naught) else: - #TODO Switch to LUT-based method when there is significant changes in the array + # TODO Switch to LUT-based method when there is significant changes in the array beta_naught = np.mean(matrix_beta_naught) return cls(azimuth_time, line, pixel, @@ -553,9 +553,10 @@ def from_product_annotation_and_aux_cal(cls, product_annotation: ProductAnnotati Ns = product_annotation.number_of_samples fs = product_annotation.range_sampling_rate eta_start = azimuth_time - tau_0 = product_annotation.antenna_pattern_slant_range_time[id_closest] + tau_0 = product_annotation.slant_range_time tau_sub = product_annotation.antenna_pattern_slant_range_time[id_closest] - theta_sub = product_annotation.antenna_pattern_elevation_pattern[id_closest] + #theta_sub = product_annotation.antenna_pattern_elevation_pattern[id_closest] + theta_sub = product_annotation.antenna_pattern_elevation_angle[id_closest] #self.theta_am = product_annotation.antenna_pattern_elevation_angle G_eap = aux_cal.elevation_antenna_pattern delta_theta = aux_cal.elevation_angle_increment @@ -567,34 +568,10 @@ def from_product_annotation_and_aux_cal(cls, product_annotation: ProductAnnotati G_eap, delta_theta) - def export_lut(self): - '''Returns LUT for EAP correction. - Parts of the codes and subfunctions were borrowed from ISCE2. - Based on ESA docuemnt : "Impact of the Elevation Antenna Pattern Phase Compensation on the Interferometric Phase Preservation"''' - - n_elt = len(self.G_eap) - - theta_am = (np.arange(n_elt)-(n_elt-1)/2)* self.delta_theta - - delta_anx = self.eta_start-self.ascending_node_time - theta_offnadir = self._anx2roll(delta_anx) - - theta_eap = theta_am+theta_offnadir - - tau = self.tau_0+np.arange(self.Ns)/self.fs - - theta = np.interp(tau, self.tau_sub, self.theta_sub) - - interpolator_G = interp1d(theta_eap,self.G_eap) - G_eap_interpolated = interpolator_G(theta) - phi_EAP = np.angle(G_eap_interpolated) - cJ = np.complex64(1.0j) - G_EAP = np.exp(cJ * phi_EAP) - - return G_EAP - def _anx2roll(self, delta_anx): ''' + Borrowed from ISCE2. The original docstring as below: + Returns the Platform nominal roll as function of elapsed time from ascending node crossing time (ANX). Straight from S1A documentation. @@ -620,9 +597,12 @@ def _anx2roll(self, delta_anx): def _anx2height(self, delta_anx): ''' + Borrowed from ISCE2. The original docstring as below: + Returns the platform nominal height as function of elapse time from ascending node crossing time (ANX). Straight from S1A documention. + ''' ###Average height @@ -635,7 +615,7 @@ def _anx2height(self, delta_anx): phi = np.array([3.1495, -1.5655 , -3.1297, 4.7222]) #;radians ###Orbital time period in seconds - Torb = (12*24*60*60)/175. + Torb = (12*24*60*60) / 175. ###Angular velocity worb = 2*np.pi / Torb diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 9059acc8..08f72e12 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -8,6 +8,7 @@ import isce3 import numpy as np from osgeo import gdal +from scipy.interpolate import InterpolatedUnivariateSpline, interp1d from s1reader import s1_annotation @@ -621,3 +622,62 @@ def swath_name(self): '''Swath name in iw1, iw2, iw3.''' return self.burst_id.split('_')[1] + @property + def thermal_noise_lut(self): + '''Returns the burst-sized LUT for thermal noise correction''' + #ncols = self.azimuth_last_range_sample - self.azimuth_first_range_sample + 1 + #nrows = self.line_to - self.line_from + 1 + ncols = self.shape[1] + nrows = self.shape[0] + + # Interpolate the range noise vector + intp_rg_lut = InterpolatedUnivariateSpline(self.burst_noise.range_pixel, + self.burst_noise.range_lut, + k=1) + if self.burst_noise.azimuth_last_range_sample is not None: + grid_rg = np.arange(self.burst_noise.azimuth_last_range_sample + 1) + else: + grid_rg = np.arange(ncols) + rg_lut_interp = intp_rg_lut(grid_rg).reshape((1, ncols)) + + # Interpolate the azimuth noise vector + if (self.burst_noise.azimuth_line is None) or (self.burst_noise.azimuth_lut is None): + az_lut_interp = np.ones(nrows).reshape((nrows, 1)) + else: # IPF >= 2.90 + intp_az_lut = InterpolatedUnivariateSpline(self.burst_noise.azimuth_line, + self.burst_noise.azimuth_lut, + k=1) + grid_az = np.arange(self.burst_noise.line_from, self.burst_noise.line_to + 1) + az_lut_interp = intp_az_lut(grid_az).reshape((nrows, 1)) + + arr_lut_total = np.matmul(az_lut_interp, rg_lut_interp) + + return arr_lut_total + + @property + def eap_compensation_lut(self): + '''Returns LUT for EAP compensation. + Based on ESA docuemnt : + "Impact of the Elevation Antenna Pattern Phase Compensation on the Interferometric Phase Preservation" + ''' + + n_elt = len(self.burst_eap.G_eap) + + theta_am = (np.arange(n_elt)-(n_elt-1)/2)* self.burst_eap.delta_theta + + delta_anx = self.burst_eap.eta_start-self.burst_eap.ascending_node_time + theta_offnadir = self.burst_eap._anx2roll(delta_anx.seconds + delta_anx.microseconds * 1.0e-6) + + theta_eap = theta_am + theta_offnadir + + tau = self.burst_eap.tau_0 + np.arange(self.burst_eap.Ns) / self.burst_eap.fs + + theta = np.interp(tau, self.burst_eap.tau_sub, self.burst_eap.theta_sub) + + interpolator_G = interp1d(theta_eap, self.burst_eap.G_eap) + G_eap_interpolated = interpolator_G(theta) + phi_EAP = np.angle(G_eap_interpolated) + cJ = np.complex64(1.0j) + G_EAP = np.exp(cJ * phi_EAP) + + return G_EAP diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index 83dd8abd..6062b619 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -1,4 +1,5 @@ import datetime +import glob import os import warnings import xml.etree.ElementTree as ET @@ -279,6 +280,30 @@ def get_ipf_version(tree: ET): return ipf_version +def get_path_aux_cal(directory_aux_cal: str, str_platform: str, inst_config_id: int): + ''' + Decide what aux_cal to load + + Parameters: + ----------- + diectory_aux_cal: str + Directory for the AUX_CAL .zip files + str_platform: str + 'S1A' or 'S1B' + insc_config_id: int + Instrument configuration ID + + Returns: + -------- + path_aux_cal: AUX_CAL file that corresponds to the criteria provided + + ''' + list_aux_cal = glob.glob(f'{directory_aux_cal}/{str_platform}_AUX_CAL_V*.SAFE.zip') + list_aux_cal.sort() + path_aux_cal = list_aux_cal[inst_config_id-1] + + return path_aux_cal + def is_eap_correction_necessary(ipf_version: version.Version) -> SimpleNamespace : '''Examines if what level of EAP correction is necessary, based on the IPF version. Based on the comment on PR: https://github.com/opera-adt/s1-reader/pull/48#discussion_r926138372 @@ -349,30 +374,46 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, # burst sensing starts in prototyping test data burst_interval = 2.758277 - #parse manifest.safe to retrieve IPF version + # 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) - #Load the Product annotation - for EAP calibration - #with open_method(annotation_path, 'r') as f_lads: - # tree_lads = ET.parse(f_lads) - # product_annotation = s1_annotation.ProductAnnotation.from_et(tree_lads) + # Load the Product annotation - for EAP calibration + with open_method(annotation_path, 'r') as f_lads: + tree_lads = ET.parse(f_lads) + product_annotation = s1_annotation.ProductAnnotation.from_et(tree_lads) - #load the Calibraton annotation + # load the Calibraton annotation calibration_annotation_path =\ annotation_path.replace('annotation/', 'annotation/calibration/calibration-') with open_method(calibration_annotation_path, 'r') as f_cads: tree_cads = ET.parse(f_cads) calibration_annotation = s1_annotation.CalibrationAnnotation.from_et(tree_cads) - #load the Noise annotation + # 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 = s1_annotation.NoiseAnnotation.from_et(tree_nads, ipf_version) + # load AUX_CAL annotation + flag_apply_eap = is_eap_correction_necessary(ipf_version) + if flag_apply_eap.phase_correction: + path_aux_cal = get_path_aux_cal(f'{os.path.dirname(s1_annotation.__file__)}/data/aux_cal', + platform_id, + product_annotation.inst_config_id) + + with zipfile.ZipFile(path_aux_cal,'r') as zipfile_aux_cal: + str_safe_aux_cal = os.path.basename(path_aux_cal).replace('.zip','') + with zipfile_aux_cal.open(f'{str_safe_aux_cal}/data/s1a-aux-cal.xml','r') as f_aux_cal: + tree_aux_cal = ET.parse(f_aux_cal) + aux_cal_subswath = s1_annotation.AuxCal.from_et(tree_aux_cal, pol, subswath_id) + else: + #No need to load aux_cal (not applying EAP correction) + aux_cal_subswath = None + # Nearly all metadata loaded here is common to all bursts in annotation XML with open_method(annotation_path, 'r') as f: tree = ET.parse(f) @@ -492,11 +533,20 @@ 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 =\ - s1_annotation.BurstCalibration.from_calibration_annotation(calibration_annotation, sensing_start) - bursts_noise=\ + burst_calibration = \ + s1_annotation.BurstCalibration.from_calibration_annotation(calibration_annotation, + sensing_start) + burst_noise = \ s1_annotation.BurstNoise.from_noise_annotation(noise_annotation, sensing_start, i*n_lines, (i+1)*n_lines-1, ipf_version) + if flag_apply_eap.phase_correction: + burst_aux_cal = \ + s1_annotation.BurstEAP.from_product_annotation_and_aux_cal(product_annotation, + aux_cal_subswath, + sensing_start) + else: + # Not applying EAP correction; No need to fill in `burst_aux_cal` + burst_aux_cal = None bursts[i] = Sentinel1BurstSlc(ipf_version, sensing_start, radar_freq, wavelength, azimuth_steer_rate, azimuth_time_interval, @@ -510,7 +560,7 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, last_sample, first_valid_line, last_line, range_window_type, range_window_coeff, rank, prf_raw_data, range_chirp_ramp_rate, - burst_calibration, bursts_noise, None) + burst_calibration, burst_noise, burst_aux_cal) return bursts From 6d11d3d5f3206e18af3aae3c6d9902c1d193359f Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 24 Aug 2022 15:54:48 -0700 Subject: [PATCH 04/49] Readibility improvement; removing unnecessary imports --- src/s1reader/s1_annotation.py | 6 +++--- src/s1reader/s1_burst_slc.py | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 4286a619..793ee561 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -10,7 +10,7 @@ import numpy as np from packaging import version -from scipy.interpolate import InterpolatedUnivariateSpline, interp1d +#from scipy.interpolate import InterpolatedUnivariateSpline, interp1d #A IPF version from which has azimuth noise vector version_threshold_azimuth_noise_vector = version.parse('2.90') @@ -377,7 +377,7 @@ def closest_block_to_azimuth_time(vector_azimuth_time: np.ndarray, @dataclass -class BurstNoise: #For thermal noise correction +class BurstNoise: '''Noise correction information for Sentinel-1 burst''' range_azimith_time: datetime.datetime range_line: float @@ -496,7 +496,7 @@ def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotati gamma = calibration_annotation.list_gamma[id_closest] dn = calibration_annotation.list_dn[id_closest] - #Check if all values in the list of beta_naught LUTs are the same + # Check if all values in the list of beta_naught LUTs are the same matrix_beta_naught = np.array(calibration_annotation.list_beta_nought) if matrix_beta_naught.min() == matrix_beta_naught.max(): beta_naught = np.min(matrix_beta_naught) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 08f72e12..ac63b60e 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -167,10 +167,10 @@ class Sentinel1BurstSlc: prf_raw_data: float # Pulse repetition frequency (PRF) of the raw data [Hz] range_chirp_rate: float # Range chirp rate [Hz] - #Correction information - burst_calibration: s1_annotation.BurstCalibration #Radiometric correction - burst_noise: s1_annotation.BurstNoise #Thermal noise correction - burst_eap: s1_annotation.BurstEAP #EAP correction + # Correction information + burst_calibration: s1_annotation.BurstCalibration # Radiometric correction + burst_noise: s1_annotation.BurstNoise # Thermal noise correction + burst_eap: s1_annotation.BurstEAP # EAP correction def as_isce3_radargrid(self): From 142ece4672f005ce2e1e7a9fd562dceed187e14e Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 24 Aug 2022 16:37:40 -0700 Subject: [PATCH 05/49] Format change on `burst_id`; keeping the absolute orbit number inside `Sentinel1BurstSlc` --- src/s1reader/s1_burst_slc.py | 1 + src/s1reader/s1_reader.py | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 9059acc8..3f6739ca 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -152,6 +152,7 @@ class Sentinel1BurstSlc: border: list # list of lon, lat coordinate tuples (in degrees) representing burst border orbit: isce3.core.Orbit orbit_direction: str + abs_orbit_number: int # Absolute orbit number # VRT params tiff_path: str # path to measurement tiff in SAFE/zip i_burst: int diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index 4317ae6a..506f9263 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -487,7 +487,7 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, last_valid_samples[last_line]) - burst_id = f't{track_number}_{id_burst}_{subswath_id.lower()}' + burst_id = f't{track_number:03d}_{id_burst}_{subswath_id.lower()}' #Extract burst-wise information for Calibration, Noise, and EAP correction @@ -502,7 +502,7 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, (n_lines, n_samples), az_fm_rate, doppler, rng_processing_bandwidth, pol, burst_id, platform_id, center_pts[i], - boundary_pts[i], orbit, orbit_direction, + boundary_pts[i], orbit, orbit_direction, orbit_number, tiff_path, i, first_valid_sample, last_sample, first_valid_line, last_line, range_window_type, range_window_coeff, From bae4227e8cef39ea593aa5a809aeb48b63f96cac Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 24 Aug 2022 18:03:04 -0700 Subject: [PATCH 06/49] updates on test_bursts.py --- tests/test_bursts.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_bursts.py b/tests/test_bursts.py index 06a85b40..11204126 100644 --- a/tests/test_bursts.py +++ b/tests/test_bursts.py @@ -28,9 +28,10 @@ def test_burst(bursts): [-2056.701472691132, 353389.9614836443, -54143009.57327797]] for i, burst in enumerate(bursts): - expected_burst_id = f't71_{151200 + i}_iw3' + expected_burst_id = f't071_{151200 + i}_iw3' assert burst.burst_id == expected_burst_id assert burst.i_burst == i + assert burst.abs_orbit_number == 32518 assert burst.radar_center_frequency == 5405000454.33435 assert burst.wavelength == 0.05546576 From 83454cb548912907922938e854fe099c9d64d49d Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Thu, 25 Aug 2022 11:44:38 -0700 Subject: [PATCH 07/49] keeping the basename of the CADS and NADS for populating RTC metadata --- src/s1reader/s1_annotation.py | 17 +++++++++++++---- src/s1reader/s1_reader.py | 9 ++++++--- 2 files changed, 19 insertions(+), 7 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 793ee561..fc49c4dd 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -5,6 +5,7 @@ from dataclasses import dataclass import datetime +import os import xml.etree.ElementTree as ET import numpy as np @@ -134,6 +135,7 @@ def _parse_vectorlist(cls, name_vector_list: str, name_vector: str, str_type: st @dataclass class CalibrationAnnotation(AnnotationBase): '''Reader for Calibration Annotation Data Set (CADS)''' + basename_annotation: str list_azimuth_time: np.ndarray list_line: list list_pixel: None @@ -143,7 +145,7 @@ class CalibrationAnnotation(AnnotationBase): list_dn: list @classmethod - def from_et(cls, et_in: ET): + def from_et(cls, et_in: ET, path_annotation: str): '''Extracts the list of calibration informaton from etree from CADS Parameters: ----------- @@ -157,6 +159,7 @@ def from_et(cls, et_in: ET): ''' cls.xml_et = et_in + cls.basename_annotation = os.path.basename(path_annotation) cls.list_azimuth_time = cls._parse_vectorlist('calibrationVectorList', 'azimuthTime', 'datetime') cls.list_line = cls._parse_vectorlist('calibrationVectorList', 'line', 'scalar_int') cls.list_pixel = cls._parse_vectorlist('calibrationVectorList', 'pixel', 'vector_int') @@ -174,6 +177,7 @@ class NoiseAnnotation(AnnotationBase): # REF: .../isce2/components/isceobj/Sensor/GRD/Sentinel1.py # : 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 @@ -186,7 +190,7 @@ class NoiseAnnotation(AnnotationBase): az_noise_azimuth_lut: np.ndarray @classmethod - def from_et(cls,et_in: ET, ipf_version: version.Version): + def from_et(cls,et_in: ET, ipf_version: version.Version, path_annotation: str): '''Extracts list of noise information from etree Parameter @@ -201,6 +205,7 @@ def from_et(cls,et_in: ET, ipf_version: version.Version): ''' cls.xml_et = et_in + cls.basename_annotation = os.path.basename(path_annotation) if ipf_version < version_threshold_azimuth_noise_vector: #legacy SAFE data cls.rg_list_azimuth_time = cls._parse_vectorlist('noiseVectorList', 'azimuthTime', 'datetime') @@ -379,6 +384,7 @@ def closest_block_to_azimuth_time(vector_azimuth_time: np.ndarray, @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 @@ -419,6 +425,7 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: ''' + 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] @@ -447,7 +454,7 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: azimuth_line = None azimuth_lut = None - return cls(range_azimith_time, range_line, range_pixel, range_lut, + 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, @@ -458,6 +465,7 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: class BurstCalibration: '''Calibration information for Sentinel-1 IW SLC burst ''' + basename_cads: str azimuth_time: datetime.datetime = None line: float = None pixel: np.ndarray = None @@ -485,6 +493,7 @@ def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotati EAP correction info for the burst ''' + basename_cads = calibration_annotation.basename_annotation id_closest = closest_block_to_azimuth_time(calibration_annotation.list_azimuth_time, azimuth_time) @@ -504,7 +513,7 @@ def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotati # TODO Switch to LUT-based method when there is significant changes in the array beta_naught = np.mean(matrix_beta_naught) - return cls(azimuth_time, line, pixel, + return cls(basename_cads, azimuth_time, line, pixel, sigma_naught, beta_naught, gamma, dn) diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index 50ad389f..e229085c 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -390,13 +390,16 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, annotation_path.replace('annotation/', 'annotation/calibration/calibration-') with open_method(calibration_annotation_path, 'r') as f_cads: tree_cads = ET.parse(f_cads) - calibration_annotation = s1_annotation.CalibrationAnnotation.from_et(tree_cads) + calibration_annotation =\ + s1_annotation.CalibrationAnnotation.from_et(tree_cads, + calibration_annotation_path) # 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 = s1_annotation.NoiseAnnotation.from_et(tree_nads, ipf_version) + noise_annotation = s1_annotation.NoiseAnnotation.from_et(tree_nads, ipf_version, + noise_annotation_path) # load AUX_CAL annotation flag_apply_eap = is_eap_correction_necessary(ipf_version) @@ -404,7 +407,7 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, path_aux_cal = get_path_aux_cal(f'{os.path.dirname(s1_annotation.__file__)}/data/aux_cal', platform_id, product_annotation.inst_config_id) - + with zipfile.ZipFile(path_aux_cal,'r') as zipfile_aux_cal: str_safe_aux_cal = os.path.basename(path_aux_cal).replace('.zip','') with zipfile_aux_cal.open(f'{str_safe_aux_cal}/data/s1a-aux-cal.xml','r') as f_aux_cal: From b29ffe2150a9e25af3e50011e700a500056ed559 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 12:28:16 -0700 Subject: [PATCH 08/49] Update src/s1reader/s1_annotation.py Readability improvement on equation Co-authored-by: Liang Yu --- src/s1reader/s1_annotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index fc49c4dd..28533e0b 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -337,7 +337,7 @@ def from_et(cls, et_in: ET, pol: str, str_swath: str): if n_val == len(arr_eap_val): #Provided in real numbers: In case of AUX_CAL for old IPFs. cls.elevation_antenna_pattern = arr_eap_val elif n_val*2 == len(arr_eap_val): #Provided in complex numbers: In case of recent IPFs e.g. 3.10 - cls.elevation_antenna_pattern = arr_eap_val[0::2]+arr_eap_val[1::2]*1.0j + cls.elevation_antenna_pattern = arr_eap_val[0::2] + arr_eap_val[1::2] * 1.0j else: raise ValueError(f'The number of values does not match. n_val={n_val}, #values in elevationAntennaPattern/values={len(arr_eap_val)}') From 2ae4fa281a479b41424d43129d2e26137b6a1e0d Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 12:29:27 -0700 Subject: [PATCH 09/49] Update src/s1reader/s1_annotation.py Removing commented out code Co-authored-by: Liang Yu --- src/s1reader/s1_annotation.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 28533e0b..6d7cdbc9 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -564,9 +564,7 @@ def from_product_annotation_and_aux_cal(cls, product_annotation: ProductAnnotati eta_start = azimuth_time tau_0 = product_annotation.slant_range_time tau_sub = product_annotation.antenna_pattern_slant_range_time[id_closest] - #theta_sub = product_annotation.antenna_pattern_elevation_pattern[id_closest] theta_sub = product_annotation.antenna_pattern_elevation_angle[id_closest] - #self.theta_am = product_annotation.antenna_pattern_elevation_angle G_eap = aux_cal.elevation_antenna_pattern delta_theta = aux_cal.elevation_angle_increment From 79a4f2078e1e66bc14f6be68196de5ccba1ece31 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 12:38:15 -0700 Subject: [PATCH 10/49] Update src/s1reader/s1_annotation.py Reverting the docstring to be split into two lines for PEP8 compliance Co-authored-by: Liang Yu --- src/s1reader/s1_annotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 6d7cdbc9..ce48b706 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -519,7 +519,7 @@ def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotati @dataclass class BurstEAP: - '''EAP correction information for Sentinel-1 IW SLC burst + '''Elevation antenna pattern (EAP) correction information for Sentinel-1 IW SLC burst ''' #from LADS Ns: int # number of samples From 79032a49e5691962b01930808a6015e69c222821 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 12:58:08 -0700 Subject: [PATCH 11/49] Update src/s1reader/s1_annotation.py Improving docstring of the code copied from isce2 Co-authored-by: Liang Yu --- src/s1reader/s1_annotation.py | 36 +++++++++++++++++++++-------------- 1 file changed, 22 insertions(+), 14 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index ce48b706..d56740ed 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -577,27 +577,35 @@ def from_product_annotation_and_aux_cal(cls, product_annotation: ProductAnnotati def _anx2roll(self, delta_anx): ''' - Borrowed from ISCE2. The original docstring as below: - Returns the Platform nominal roll as function of elapsed time from - ascending node crossing time (ANX). - Straight from S1A documentation. + ascending node crossing time (ANX). (Implemented from S1A documentation.) + + Code copied from ISCE2. + + Parameters + ---------- + delta_anx: float + elapsed time from ascending node crossing time + + Returns + ------- + _: float + Estimated nominal roll (degrees) ''' - - ####Estimate altitude based on time elapsed since ANX + # Estimate altitude based on time elapsed since ANX altitude = self._anx2height(delta_anx) - ####Reference altitude - href = 711.700 #;km + # Reference altitude (km) + href = 711.700 - ####Reference boresight at reference altitude - boresight_ref = 29.450 # ; deg + # Reference boresight at reference altitude (degrees) + boresight_ref = 29.450 - ####Partial derivative of roll vs altitude - alpha_roll = 0.0566 # ;deg/km + # Partial derivative of roll vs altitude (degrees/km) + alpha_roll = 0.0566 - ####Estimate nominal roll - nominal_roll = boresight_ref - alpha_roll * (altitude/1000.0 - href) #Theta off nadir + # Estimate nominal roll i.e. theta off nadir (degrees) + nominal_roll = boresight_ref - alpha_roll * (altitude/1000.0 - href) return nominal_roll From de49003405db451c0497389c849ab24680aeefe0 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 12:59:42 -0700 Subject: [PATCH 12/49] Update src/s1reader/s1_burst_slc.py Removing the commented out codes Co-authored-by: Liang Yu --- src/s1reader/s1_burst_slc.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 5cf50d29..be14aa5d 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -626,8 +626,6 @@ def swath_name(self): @property def thermal_noise_lut(self): '''Returns the burst-sized LUT for thermal noise correction''' - #ncols = self.azimuth_last_range_sample - self.azimuth_first_range_sample + 1 - #nrows = self.line_to - self.line_from + 1 ncols = self.shape[1] nrows = self.shape[0] From 9013e4f395b248119078fadc2898b75c172962fc Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 13:00:30 -0700 Subject: [PATCH 13/49] Update src/s1reader/s1_burst_slc.py improvement on code brevity Co-authored-by: Liang Yu --- src/s1reader/s1_burst_slc.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index be14aa5d..8e3d0803 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -626,8 +626,7 @@ def swath_name(self): @property def thermal_noise_lut(self): '''Returns the burst-sized LUT for thermal noise correction''' - ncols = self.shape[1] - nrows = self.shape[0] + nrows, ncols = self.shape # Interpolate the range noise vector intp_rg_lut = InterpolatedUnivariateSpline(self.burst_noise.range_pixel, From 64f34e0059b271eda474dc208bd40092b697517d Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 13:04:09 -0700 Subject: [PATCH 14/49] Update src/s1reader/s1_burst_slc.py renaming variable for better clarity Co-authored-by: Liang Yu --- src/s1reader/s1_burst_slc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 8e3d0803..9de15621 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -629,7 +629,7 @@ def thermal_noise_lut(self): nrows, ncols = self.shape # Interpolate the range noise vector - intp_rg_lut = InterpolatedUnivariateSpline(self.burst_noise.range_pixel, + rg_lut_interp_obj = InterpolatedUnivariateSpline(self.burst_noise.range_pixel, self.burst_noise.range_lut, k=1) if self.burst_noise.azimuth_last_range_sample is not None: From e7df7b07b5b049f172d4cacbbb6a5d322bea7df4 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 13:10:30 -0700 Subject: [PATCH 15/49] Update src/s1reader/s1_burst_slc.py renaming variable name for clarity Co-authored-by: Liang Yu --- src/s1reader/s1_burst_slc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 9de15621..c8aa2fab 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -642,7 +642,7 @@ def thermal_noise_lut(self): if (self.burst_noise.azimuth_line is None) or (self.burst_noise.azimuth_lut is None): az_lut_interp = np.ones(nrows).reshape((nrows, 1)) else: # IPF >= 2.90 - intp_az_lut = InterpolatedUnivariateSpline(self.burst_noise.azimuth_line, + az_lut_intp_obj = InterpolatedUnivariateSpline(self.burst_noise.azimuth_line, self.burst_noise.azimuth_lut, k=1) grid_az = np.arange(self.burst_noise.line_from, self.burst_noise.line_to + 1) From bf029cc993b1adc28b2ecd587c6631c9d05a1b84 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 13:37:30 -0700 Subject: [PATCH 16/49] Update src/s1reader/s1_burst_slc.py variable name revised for clarity Co-authored-by: Liang Yu --- src/s1reader/s1_burst_slc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index c8aa2fab..07e6ea15 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -636,7 +636,7 @@ def thermal_noise_lut(self): grid_rg = np.arange(self.burst_noise.azimuth_last_range_sample + 1) else: grid_rg = np.arange(ncols) - rg_lut_interp = intp_rg_lut(grid_rg).reshape((1, ncols)) + rg_lut_interpolated = intp_rg_lut(grid_rg).reshape((1, ncols)) # Interpolate the azimuth noise vector if (self.burst_noise.azimuth_line is None) or (self.burst_noise.azimuth_lut is None): From b3f36127cd99d3f88c5046928873f10e8dcce546 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 13:39:34 -0700 Subject: [PATCH 17/49] Update src/s1reader/s1_burst_slc.py variable renamed for clarity Co-authored-by: Liang Yu --- src/s1reader/s1_burst_slc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 07e6ea15..cc089635 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -640,7 +640,7 @@ def thermal_noise_lut(self): # Interpolate the azimuth noise vector if (self.burst_noise.azimuth_line is None) or (self.burst_noise.azimuth_lut is None): - az_lut_interp = np.ones(nrows).reshape((nrows, 1)) + az_lut_interpolated = np.ones(nrows).reshape((nrows, 1)) else: # IPF >= 2.90 az_lut_intp_obj = InterpolatedUnivariateSpline(self.burst_noise.azimuth_line, self.burst_noise.azimuth_lut, From f0b0b82fdaf809bb63bc2dcd2849da5e0cb95502 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 13:40:59 -0700 Subject: [PATCH 18/49] Update src/s1reader/s1_burst_slc.py improvement on docstring Co-authored-by: Liang Yu --- src/s1reader/s1_burst_slc.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index cc089635..bf4dcaf5 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -625,7 +625,12 @@ def swath_name(self): @property def thermal_noise_lut(self): - '''Returns the burst-sized LUT for thermal noise correction''' + '''Returns the burst-sized LUT for thermal noise correction + + Returns + ------- + arr_lut_total: np.array + 2d array containing thermal noise correction look up table values nrows, ncols = self.shape # Interpolate the range noise vector From 4109740f079c65346fb97b4c696e9086322ffd28 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 13:44:27 -0700 Subject: [PATCH 19/49] Update src/s1reader/s1_burst_slc.py Readability improvement of equation Co-authored-by: Liang Yu --- src/s1reader/s1_burst_slc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index bf4dcaf5..81a96e42 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -666,7 +666,7 @@ def eap_compensation_lut(self): n_elt = len(self.burst_eap.G_eap) - theta_am = (np.arange(n_elt)-(n_elt-1)/2)* self.burst_eap.delta_theta + theta_am = (np.arange(n_elt) - (n_elt - 1) / 2) * self.burst_eap.delta_theta delta_anx = self.burst_eap.eta_start-self.burst_eap.ascending_node_time theta_offnadir = self.burst_eap._anx2roll(delta_anx.seconds + delta_anx.microseconds * 1.0e-6) From 6d62feb63423041aad082bef12e7ce2990c01568 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 13:51:59 -0700 Subject: [PATCH 20/49] addressing comments bt @LiangJYu --- src/s1reader/s1_annotation.py | 45 ++++++++++++++++++++++------------- src/s1reader/s1_burst_slc.py | 29 ++++++++++++---------- 2 files changed, 44 insertions(+), 30 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index fc49c4dd..9c17c76c 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -519,11 +519,12 @@ def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotati @dataclass class BurstEAP: - '''EAP correction information for Sentinel-1 IW SLC burst + '''Elevation antenna pattern (EAP) correction information for + Sentinel-1 IW SLC burst ''' #from LADS - Ns: int # number of samples - fs: float # range sampling rate + num_sample: int # number of samples + freq_sampling: float # range sampling rate eta_start: datetime.datetime tau_0: float # imageInformation/slantRangeTime tau_sub: np.ndarray # antennaPattern/slantRangeTime @@ -532,7 +533,7 @@ class BurstEAP: ascending_node_time: datetime.datetime #from AUX_CAL - G_eap: np.ndarray # elevationAntennaPattern + gain_eap: np.ndarray # elevationAntennaPattern delta_theta:float # elavationAngleIncrement @classmethod @@ -558,23 +559,22 @@ def from_product_annotation_and_aux_cal(cls, product_annotation: ProductAnnotati A burst-wide information for EAP correction ''' - id_closest = closest_block_to_azimuth_time(product_annotation.antenna_pattern_azimuth_time, azimuth_time) - Ns = product_annotation.number_of_samples - fs = product_annotation.range_sampling_rate + id_closest = closest_block_to_azimuth_time(product_annotation.antenna_pattern_azimuth_time, + azimuth_time) + num_sample = product_annotation.number_of_samples + freq_sampling = product_annotation.range_sampling_rate eta_start = azimuth_time tau_0 = product_annotation.slant_range_time tau_sub = product_annotation.antenna_pattern_slant_range_time[id_closest] - #theta_sub = product_annotation.antenna_pattern_elevation_pattern[id_closest] theta_sub = product_annotation.antenna_pattern_elevation_angle[id_closest] - #self.theta_am = product_annotation.antenna_pattern_elevation_angle - G_eap = aux_cal.elevation_antenna_pattern + gain_eap = aux_cal.elevation_antenna_pattern delta_theta = aux_cal.elevation_angle_increment ascending_node_time = product_annotation.ascending_node_time - return cls(Ns, fs, eta_start, tau_0, tau_sub, theta_sub, + return cls(num_sample, freq_sampling, eta_start, tau_0, tau_sub, theta_sub, azimuth_time, ascending_node_time, - G_eap, delta_theta) + gain_eap, delta_theta) def _anx2roll(self, delta_anx): @@ -606,12 +606,24 @@ def _anx2roll(self, delta_anx): def _anx2height(self, delta_anx): ''' - Borrowed from ISCE2. The original docstring as below: - Returns the platform nominal height as function of elapse time from ascending node crossing time (ANX). Straight from S1A documention. + Code copied from ISCE2. + + + Parameters: + ----------- + delta_anx: float + elapsed time from ANX time + + + Returns: + -------- + _: float + nominal height of the platform + ''' ###Average height @@ -624,14 +636,13 @@ def _anx2height(self, delta_anx): phi = np.array([3.1495, -1.5655 , -3.1297, 4.7222]) #;radians ###Orbital time period in seconds - Torb = (12*24*60*60) / 175. + t_orb = (12*24*60*60) / 175. ###Angular velocity - worb = 2*np.pi / Torb + worb = 2*np.pi / t_orb ####Evaluation of series h_t = h_0 - #for i in range(len(h)): for i, h_i in enumerate(h): h_t += h_i * np.sin((i+1) * worb * delta_anx + phi[i]) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 5cf50d29..b94d2b67 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -636,10 +636,10 @@ def thermal_noise_lut(self): self.burst_noise.range_lut, k=1) if self.burst_noise.azimuth_last_range_sample is not None: - grid_rg = np.arange(self.burst_noise.azimuth_last_range_sample + 1) + vec_rg = np.arange(self.burst_noise.azimuth_last_range_sample + 1) else: - grid_rg = np.arange(ncols) - rg_lut_interp = intp_rg_lut(grid_rg).reshape((1, ncols)) + vec_rg = np.arange(ncols) + rg_lut_interp = intp_rg_lut(vec_rg).reshape((1, ncols)) # Interpolate the azimuth noise vector if (self.burst_noise.azimuth_line is None) or (self.burst_noise.azimuth_lut is None): @@ -651,7 +651,7 @@ def thermal_noise_lut(self): grid_az = np.arange(self.burst_noise.line_from, self.burst_noise.line_to + 1) az_lut_interp = intp_az_lut(grid_az).reshape((nrows, 1)) - arr_lut_total = np.matmul(az_lut_interp, rg_lut_interp) + arr_lut_total = np.matmul(az_lut_interp, rg_lut_interpolated) return arr_lut_total @@ -659,10 +659,13 @@ def thermal_noise_lut(self): def eap_compensation_lut(self): '''Returns LUT for EAP compensation. Based on ESA docuemnt : - "Impact of the Elevation Antenna Pattern Phase Compensation on the Interferometric Phase Preservation" + "Impact of the Elevation Antenna Pattern Phase Compensation + on the Interferometric Phase Preservation" + Document URL: + https://sentinel.esa.int/documents/247904/1653440/Sentinel-1-IPF_EAP_Phase_correction ''' - - n_elt = len(self.burst_eap.G_eap) + + n_elt = len(self.burst_eap.gain_eap) theta_am = (np.arange(n_elt)-(n_elt-1)/2)* self.burst_eap.delta_theta @@ -671,14 +674,14 @@ def eap_compensation_lut(self): theta_eap = theta_am + theta_offnadir - tau = self.burst_eap.tau_0 + np.arange(self.burst_eap.Ns) / self.burst_eap.fs + tau = self.burst_eap.tau_0 + np.arange(self.burst_eap.num_sample) / self.burst_eap.freq_sampling theta = np.interp(tau, self.burst_eap.tau_sub, self.burst_eap.theta_sub) - interpolator_G = interp1d(theta_eap, self.burst_eap.G_eap) - G_eap_interpolated = interpolator_G(theta) - phi_EAP = np.angle(G_eap_interpolated) + interpolator_gain = interp1d(theta_eap, self.burst_eap.gain_eap) + gain_eap_interpolated = interpolator_gain(theta) + phi_eap = np.angle(gain_eap_interpolated) cJ = np.complex64(1.0j) - G_EAP = np.exp(cJ * phi_EAP) + gain_eap = np.exp(cJ * phi_eap) - return G_EAP + return gain_eap From 1af92bf962542598aa38a10c15e7b0b4d11a6eb2 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 14:08:42 -0700 Subject: [PATCH 21/49] docstring fix; variables renamed for clarity --- src/s1reader/s1_burst_slc.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 24d3d164..6acc4b5e 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -626,11 +626,13 @@ def swath_name(self): @property def thermal_noise_lut(self): '''Returns the burst-sized LUT for thermal noise correction - + Returns ------- arr_lut_total: np.array 2d array containing thermal noise correction look up table values + ''' + nrows, ncols = self.shape # Interpolate the range noise vector @@ -641,19 +643,20 @@ def thermal_noise_lut(self): vec_rg = np.arange(self.burst_noise.azimuth_last_range_sample + 1) else: vec_rg = np.arange(ncols) - rg_lut_interp = intp_rg_lut(vec_rg).reshape((1, ncols)) + rg_lut_interpolated = rg_lut_interp_obj(vec_rg).reshape((1, ncols)) + # Interpolate the azimuth noise vector if (self.burst_noise.azimuth_line is None) or (self.burst_noise.azimuth_lut is None): az_lut_interpolated = np.ones(nrows).reshape((nrows, 1)) else: # IPF >= 2.90 - az_lut_intp_obj = InterpolatedUnivariateSpline(self.burst_noise.azimuth_line, + az_lut_interp_obj = InterpolatedUnivariateSpline(self.burst_noise.azimuth_line, self.burst_noise.azimuth_lut, k=1) - grid_az = np.arange(self.burst_noise.line_from, self.burst_noise.line_to + 1) - az_lut_interp = intp_az_lut(grid_az).reshape((nrows, 1)) + vec_az = np.arange(self.burst_noise.line_from, self.burst_noise.line_to + 1) + az_lut_interpolated = az_lut_interp_obj(vec_az).reshape((nrows, 1)) - arr_lut_total = np.matmul(az_lut_interp, rg_lut_interpolated) + arr_lut_total = np.matmul(az_lut_interpolated, rg_lut_interpolated) return arr_lut_total From d9a76d90c7589ef8e6ee9675f0f4a8a1a1368c65 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 15:25:42 -0700 Subject: [PATCH 22/49] implemented s1_annotation.AucCal.load_from_zip_file() --- src/s1reader/s1_annotation.py | 12 ++++++++---- src/s1reader/s1_reader.py | 17 ++++++++--------- 2 files changed, 16 insertions(+), 13 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 85614218..e29ac4a3 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -7,6 +7,7 @@ import datetime import os import xml.etree.ElementTree as ET +import zipfile import numpy as np @@ -304,13 +305,13 @@ class AuxCal(AnnotationBase): @classmethod - def from_et(cls, et_in: ET, pol: str, str_swath: str): + def load_from_zip_file(cls, path_aux_cal_zip: str, pol: str, str_swath: str): '''A class method that Extracts list of information AUX_CAL from the input ET. Parameters --------- - et_in : ET - ET from AUX_CAL + et_in : path_aux_cal_zip + Path to the AUX_CAL .zip file pol: str {'vv','vh','hh','hv'} Polarization of interest str_swath: {'iw1','iw2','iw3'} @@ -321,6 +322,10 @@ def from_et(cls, et_in: ET, pol: str, str_swath: str): cls: AuxCal class populated by et_in in the parameter ''' + str_safe_aux_cal = os.path.basename(path_aux_cal_zip).replace('.zip','') + with zipfile.ZipFile(path_aux_cal_zip, 'r') as zipfile_aux_cal: + with zipfile_aux_cal.open(f'{str_safe_aux_cal}/data/s1a-aux-cal.xml','r') as f_aux_cal: + et_in = ET.parse(f_aux_cal) calibration_params_list = et_in.find('calibrationParamsList') for calibration_params in calibration_params_list: @@ -358,7 +363,6 @@ def from_et(cls, et_in: ET, pol: str, str_swath: str): return cls - def closest_block_to_azimuth_time(vector_azimuth_time: np.ndarray, azimuth_time_burst: datetime.datetime) -> int: '''Find the id of the closest data block in annotation. diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index e229085c..fc8765d2 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -404,15 +404,14 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, # load AUX_CAL annotation flag_apply_eap = is_eap_correction_necessary(ipf_version) if flag_apply_eap.phase_correction: - path_aux_cal = get_path_aux_cal(f'{os.path.dirname(s1_annotation.__file__)}/data/aux_cal', - platform_id, - product_annotation.inst_config_id) - - with zipfile.ZipFile(path_aux_cal,'r') as zipfile_aux_cal: - str_safe_aux_cal = os.path.basename(path_aux_cal).replace('.zip','') - with zipfile_aux_cal.open(f'{str_safe_aux_cal}/data/s1a-aux-cal.xml','r') as f_aux_cal: - tree_aux_cal = ET.parse(f_aux_cal) - aux_cal_subswath = s1_annotation.AuxCal.from_et(tree_aux_cal, pol, subswath_id) + path_aux_cal_zip = get_path_aux_cal(f'{os.path.dirname(s1_annotation.__file__)}' + '/data/aux_cal', + platform_id, + product_annotation.inst_config_id) + + aux_cal_subswath = s1_annotation.AuxCal.load_from_zip_file(path_aux_cal_zip, + pol, + subswath_id) else: #No need to load aux_cal (not applying EAP correction) aux_cal_subswath = None From 3504967289de8326e746f1bce9c9f192b07780bf Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 15:26:04 -0700 Subject: [PATCH 23/49] readability improvement --- src/s1reader/s1_burst_slc.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 6acc4b5e..474b60ce 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -637,8 +637,8 @@ def thermal_noise_lut(self): # Interpolate the range noise vector rg_lut_interp_obj = InterpolatedUnivariateSpline(self.burst_noise.range_pixel, - self.burst_noise.range_lut, - k=1) + self.burst_noise.range_lut, + k=1) if self.burst_noise.azimuth_last_range_sample is not None: vec_rg = np.arange(self.burst_noise.azimuth_last_range_sample + 1) else: @@ -651,8 +651,8 @@ def thermal_noise_lut(self): az_lut_interpolated = np.ones(nrows).reshape((nrows, 1)) else: # IPF >= 2.90 az_lut_interp_obj = InterpolatedUnivariateSpline(self.burst_noise.azimuth_line, - self.burst_noise.azimuth_lut, - k=1) + self.burst_noise.azimuth_lut, + k=1) vec_az = np.arange(self.burst_noise.line_from, self.burst_noise.line_to + 1) az_lut_interpolated = az_lut_interp_obj(vec_az).reshape((nrows, 1)) @@ -666,6 +666,7 @@ def eap_compensation_lut(self): Based on ESA docuemnt : "Impact of the Elevation Antenna Pattern Phase Compensation on the Interferometric Phase Preservation" + Document URL: https://sentinel.esa.int/documents/247904/1653440/Sentinel-1-IPF_EAP_Phase_correction ''' @@ -686,7 +687,7 @@ def eap_compensation_lut(self): interpolator_gain = interp1d(theta_eap, self.burst_eap.gain_eap) gain_eap_interpolated = interpolator_gain(theta) phi_eap = np.angle(gain_eap_interpolated) - cJ = np.complex64(1.0j) - gain_eap = np.exp(cJ * phi_eap) + c_j = np.complex64(1.0j) + gain_eap = np.exp(c_j * phi_eap) return gain_eap From 6b2a07aacb43bcb2ab7a4a82be931e2a9fc7fab7 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 16:15:17 -0700 Subject: [PATCH 24/49] s1_annotation.py - code cleanup; excention handling for AUX_CAL; PEP8 compliance --- src/s1reader/s1_annotation.py | 38 +++++++++++++++++++---------------- 1 file changed, 21 insertions(+), 17 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index e29ac4a3..014d913c 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -12,7 +12,6 @@ import numpy as np from packaging import version -#from scipy.interpolate import InterpolatedUnivariateSpline, interp1d #A IPF version from which has azimuth noise vector version_threshold_azimuth_noise_vector = version.parse('2.90') @@ -322,7 +321,11 @@ def load_from_zip_file(cls, path_aux_cal_zip: str, pol: str, str_swath: str): cls: AuxCal class populated by et_in in the parameter ''' - str_safe_aux_cal = os.path.basename(path_aux_cal_zip).replace('.zip','') + if os.path.exists(path_aux_cal_zip): + str_safe_aux_cal = os.path.basename(path_aux_cal_zip).replace('.zip','') + else: + raise ValueError(f'Cannot find AUX_CAL .zip file: {path_aux_cal_zip}') + with zipfile.ZipFile(path_aux_cal_zip, 'r') as zipfile_aux_cal: with zipfile_aux_cal.open(f'{str_safe_aux_cal}/data/s1a-aux-cal.xml','r') as f_aux_cal: et_in = ET.parse(f_aux_cal) @@ -430,7 +433,8 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: ''' basename_nads = noise_annotation.basename_annotation - id_closest = closest_block_to_azimuth_time(noise_annotation.rg_list_azimuth_time, azimuth_time) + 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] @@ -443,10 +447,11 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: azimuth_last_range_sample = noise_annotation.az_last_range_sample if ipf_version >= version_threshold_azimuth_noise_vector: - #Azimuth noise LUT exists - crop to the extent of the burst + # 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 + + # put some margin when possible if id_top > 0: id_top -= 1 if id_bottom < len(noise_annotation.az_line)-1: @@ -505,7 +510,6 @@ def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotati line = calibration_annotation.list_line[id_closest] pixel = calibration_annotation.list_pixel[id_closest] sigma_naught = calibration_annotation.list_sigma_nought[id_closest] - #beta_naught = calibration_annotation.list_beta_nought[id_closest] gamma = calibration_annotation.list_gamma[id_closest] dn = calibration_annotation.list_dn[id_closest] @@ -526,7 +530,7 @@ class BurstEAP: '''Elevation antenna pattern (EAP) correction information for Sentinel-1 IW SLC burst ''' - #from LADS + # from LADS num_sample: int # number of samples freq_sampling: float # range sampling rate eta_start: datetime.datetime @@ -536,7 +540,7 @@ class BurstEAP: azimuth_time: datetime.datetime ascending_node_time: datetime.datetime - #from AUX_CAL + # from AUX_CAL gain_eap: np.ndarray # elevationAntennaPattern delta_theta:float # elavationAngleIncrement @@ -585,14 +589,14 @@ def _anx2roll(self, delta_anx): ''' Returns the Platform nominal roll as function of elapsed time from ascending node crossing time (ANX). (Implemented from S1A documentation.) - + Code copied from ISCE2. - + Parameters ---------- delta_anx: float elapsed time from ascending node crossing time - + Returns ------- _: float @@ -638,22 +642,22 @@ def _anx2height(self, delta_anx): ''' - ###Average height + ### Average height h_0 = 707714.8 #;m - ####Perturbation amplitudes + #### Perturbation amplitudes h = np.array([8351.5, 8947.0, 23.32, 11.74]) #;m - ####Perturbation phases + #### Perturbation phases phi = np.array([3.1495, -1.5655 , -3.1297, 4.7222]) #;radians - ###Orbital time period in seconds + ###O rbital time period in seconds t_orb = (12*24*60*60) / 175. - ###Angular velocity + ### Angular velocity worb = 2*np.pi / t_orb - ####Evaluation of series + #### Evaluation of series h_t = h_0 for i, h_i in enumerate(h): h_t += h_i * np.sin((i+1) * worb * delta_anx + phi[i]) From d5cd5de4f9fa8e7dcb4a14d8e492688707580f1d Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 16:15:57 -0700 Subject: [PATCH 25/49] docstring for `s1_burst_slc.eap_compensation_lut()` --- src/s1reader/s1_burst_slc.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 474b60ce..160dc016 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -669,6 +669,11 @@ def eap_compensation_lut(self): Document URL: https://sentinel.esa.int/documents/247904/1653440/Sentinel-1-IPF_EAP_Phase_correction + + Returns: + ------- + gain_eap: EAP phase for the burst to be compensated + ''' n_elt = len(self.burst_eap.gain_eap) From db59419e57e5d07d024aed0262c80b6ee2db9a1c Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 16:18:03 -0700 Subject: [PATCH 26/49] class import scheme changed --- src/s1reader/s1_reader.py | 42 ++++++++++++++++++++------------------- 1 file changed, 22 insertions(+), 20 deletions(-) diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index fc8765d2..e9e85ff6 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -13,8 +13,13 @@ import shapely 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_burst_slc import Doppler, Sentinel1BurstSlc -from s1reader import s1_annotation esa_track_burst_id_file = f"{os.path.dirname(os.path.realpath(__file__))}/data/sentinel1_track_burst_id.txt" @@ -383,7 +388,7 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, # Load the Product annotation - for EAP calibration with open_method(annotation_path, 'r') as f_lads: tree_lads = ET.parse(f_lads) - product_annotation = s1_annotation.ProductAnnotation.from_et(tree_lads) + product_annotation = ProductAnnotation.from_et(tree_lads) # load the Calibraton annotation calibration_annotation_path =\ @@ -391,15 +396,15 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, with open_method(calibration_annotation_path, 'r') as f_cads: tree_cads = ET.parse(f_cads) calibration_annotation =\ - s1_annotation.CalibrationAnnotation.from_et(tree_cads, - calibration_annotation_path) + CalibrationAnnotation.from_et(tree_cads, + calibration_annotation_path) # 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 = s1_annotation.NoiseAnnotation.from_et(tree_nads, ipf_version, - noise_annotation_path) + noise_annotation = NoiseAnnotation.from_et(tree_nads, ipf_version, + noise_annotation_path) # load AUX_CAL annotation flag_apply_eap = is_eap_correction_necessary(ipf_version) @@ -409,11 +414,11 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, platform_id, product_annotation.inst_config_id) - aux_cal_subswath = s1_annotation.AuxCal.load_from_zip_file(path_aux_cal_zip, - pol, - subswath_id) + aux_cal_subswath = AuxCal.load_from_zip_file(path_aux_cal_zip, + pol, + subswath_id) else: - #No need to load aux_cal (not applying EAP correction) + # No need to load aux_cal (not applying EAP correction) aux_cal_subswath = None # Nearly all metadata loaded here is common to all bursts in annotation XML @@ -535,17 +540,14 @@ 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 = \ - s1_annotation.BurstCalibration.from_calibration_annotation(calibration_annotation, - sensing_start) - burst_noise = \ - s1_annotation.BurstNoise.from_noise_annotation(noise_annotation, sensing_start, - i*n_lines, (i+1)*n_lines-1, ipf_version) + 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) if flag_apply_eap.phase_correction: - burst_aux_cal = \ - s1_annotation.BurstEAP.from_product_annotation_and_aux_cal(product_annotation, - aux_cal_subswath, - sensing_start) + burst_aux_cal = BurstEAP.from_product_annotation_and_aux_cal(product_annotation, + aux_cal_subswath, + sensing_start) else: # Not applying EAP correction; No need to fill in `burst_aux_cal` burst_aux_cal = None From e908892e5c03c46fb3a1c0314ea4d796040eed17 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sun, 28 Aug 2022 16:29:08 -0700 Subject: [PATCH 27/49] PEP8 compliance --- src/s1reader/s1_annotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 014d913c..f343b5fc 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -271,7 +271,7 @@ def from_et(cls, et_in: ET): cls.xml_et = et_in cls.image_information_slant_range_time = cls._parse_scalar('imageAnnotation/imageInformation/slantRangeTime', 'scalar_float') - cls.inst_config_id = cls._parse_scalar('generalAnnotation/downlinkInformationList/downlinkInformation/downlinkValues/instrumentConfigId', 'scalar_int') + cls.inst_config_id = cls._parse_scalar('generalAnnotation/downlinkInformationList/downlinkInformation/downlinkValues/instrumentConfigId', 'scalar_int') cls.antenna_pattern_azimuth_time = cls._parse_vectorlist('antennaPattern/antennaPatternList', 'azimuthTime', 'datetime') cls.antenna_pattern_slant_range_time = cls._parse_vectorlist('antennaPattern/antennaPatternList', 'slantRangeTime', 'vector_float') cls.antenna_pattern_elevation_angle = cls._parse_vectorlist('antennaPattern/antennaPatternList', 'elevationAngle', 'vector_float') From a9149e1d69467077576fea63cef57ff64cc3acc6 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 31 Aug 2022 08:42:17 -0700 Subject: [PATCH 28/49] Update src/s1reader/s1_annotation.py Renaming variable for clarity Co-authored-by: Liang Yu --- src/s1reader/s1_annotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index f343b5fc..7522de44 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -240,7 +240,7 @@ class ProductAnnotation(AnnotationBase): image_information_slant_range_time: float # Attributes to be used when determining what AUX_CAL to load - inst_config_id: int + instrument_cfg_id: int #elevation_angle: antenna_pattern_azimuth_time: list From 7fb70a4e77803e89658273f2734da3bc3a57f101 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 31 Aug 2022 08:45:32 -0700 Subject: [PATCH 29/49] Update src/s1reader/s1_annotation.py Fix on docstring for `s1_annotation.AuxCal.load_from_zip_file()` Co-authored-by: Liang Yu --- src/s1reader/s1_annotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 7522de44..d15de025 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -309,7 +309,7 @@ def load_from_zip_file(cls, path_aux_cal_zip: str, pol: str, str_swath: str): Parameters --------- - et_in : path_aux_cal_zip + path_aux_cal_zip : str Path to the AUX_CAL .zip file pol: str {'vv','vh','hh','hv'} Polarization of interest From 53ffce4a2af1f178e42a4a22a10010e0d2644f96 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 31 Aug 2022 08:49:05 -0700 Subject: [PATCH 30/49] Update src/s1reader/s1_annotation.py Change on docstring for s1_annotation.BurstEAP._anx2height Co-authored-by: Liang Yu --- src/s1reader/s1_annotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index d15de025..a5f2666e 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -624,7 +624,7 @@ def _anx2height(self, delta_anx): ''' Returns the platform nominal height as function of elapse time from ascending node crossing time (ANX). - Straight from S1A documention. + Implementation from S1A documention. Code copied from ISCE2. From 0700e74f29775eed815e588ab2076b7d9b7ec715 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 31 Aug 2022 09:00:56 -0700 Subject: [PATCH 31/49] Update src/s1reader/s1_reader.py updates on s1_reader.get_path_aux_cal() Co-authored-by: Liang Yu --- src/s1reader/s1_reader.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index e9e85ff6..87dc6d2c 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -285,9 +285,9 @@ def get_ipf_version(tree: ET): return ipf_version -def get_path_aux_cal(directory_aux_cal: str, str_platform: str, inst_config_id: int): +def get_path_aux_cal(directory_aux_cal: str, str_platform: str, instrument_cfg_id: int): ''' - Decide what aux_cal to load + Decide which aux_cal to load Parameters: ----------- @@ -295,13 +295,12 @@ def get_path_aux_cal(directory_aux_cal: str, str_platform: str, inst_config_id: Directory for the AUX_CAL .zip files str_platform: str 'S1A' or 'S1B' - insc_config_id: int + instrument_cfg_id: int Instrument configuration ID Returns: -------- path_aux_cal: AUX_CAL file that corresponds to the criteria provided - ''' list_aux_cal = glob.glob(f'{directory_aux_cal}/{str_platform}_AUX_CAL_V*.SAFE.zip') list_aux_cal.sort() From fe8252e7fa8b9ba97163510f263c9fcbd62ece35 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 31 Aug 2022 09:03:35 -0700 Subject: [PATCH 32/49] Update src/s1reader/s1_reader.py Readability improvement Co-authored-by: Liang Yu --- src/s1reader/s1_reader.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index 87dc6d2c..895dc5b3 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -408,8 +408,8 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, # load AUX_CAL annotation flag_apply_eap = is_eap_correction_necessary(ipf_version) if flag_apply_eap.phase_correction: - path_aux_cal_zip = get_path_aux_cal(f'{os.path.dirname(s1_annotation.__file__)}' - '/data/aux_cal', + path_aux_cals = f'{os.path.dirname(s1_annotation.__file__)}/data/aux_cal' + path_aux_cal_zip = get_path_aux_cal(path_aux_cals, platform_id, product_annotation.inst_config_id) From 90c09bd367f499739f5daaa7b1f32b7e458e33ae Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 31 Aug 2022 12:37:23 -0700 Subject: [PATCH 33/49] command and variable name revised for clarity --- src/s1reader/s1_annotation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index a5f2666e..2fb2091c 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -13,8 +13,8 @@ from packaging import version -#A IPF version from which has azimuth noise vector -version_threshold_azimuth_noise_vector = version.parse('2.90') +#IPF version from which the NADS has azimuth noise vector annotation +ipf_version_az_noise_vector_available_from = version.parse('2.90') @dataclass class AnnotationBase: From d87c234e039889b66ee861f001b2a52d245310cd Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 31 Aug 2022 12:48:05 -0700 Subject: [PATCH 34/49] Raise error when `burst_noise` is not defined but user attempts to retrieve thermal correction LUT --- src/s1reader/s1_burst_slc.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 160dc016..84dee9a1 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -633,6 +633,9 @@ def thermal_noise_lut(self): 2d array containing thermal noise correction look up table values ''' + if self.burst_noise is None: + raise ValueError('burst_noise is not defined for this burst.') + nrows, ncols = self.shape # Interpolate the range noise vector From 08acee0432e80244d8523ba5043c9d538e69b428 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 31 Aug 2022 13:01:20 -0700 Subject: [PATCH 35/49] fixing variable names that causing CircieCI to fail --- src/s1reader/s1_annotation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 2fb2091c..081971b4 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -207,7 +207,7 @@ def from_et(cls,et_in: ET, ipf_version: version.Version, path_annotation: str): cls.xml_et = et_in cls.basename_annotation = os.path.basename(path_annotation) - if ipf_version < version_threshold_azimuth_noise_vector: #legacy SAFE data + if ipf_version < ipf_version_az_noise_vector_available_from: #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') @@ -446,7 +446,7 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: azimuth_last_azimuth_line = noise_annotation.az_last_azimuth_line azimuth_last_range_sample = noise_annotation.az_last_range_sample - if ipf_version >= version_threshold_azimuth_noise_vector: + if ipf_version >= ipf_version_az_noise_vector_available_from: # 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)) From 27f6998009f123d9410577f6e402d02f1bb0c75c Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Tue, 6 Sep 2022 15:18:12 -0700 Subject: [PATCH 36/49] LUT retrieval scheme changed; updates in accordance with the comments / suggestions from @Liang Yu --- src/s1reader/s1_annotation.py | 90 +++++++++++++++++++++++++++++++++-- src/s1reader/s1_burst_slc.py | 71 ++++----------------------- src/s1reader/s1_reader.py | 3 +- 3 files changed, 96 insertions(+), 68 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 081971b4..5f1ca8af 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -11,6 +11,7 @@ import numpy as np +from scipy.interpolate import InterpolatedUnivariateSpline, interp1d from packaging import version #IPF version from which the NADS has azimuth noise vector annotation @@ -470,6 +471,49 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: 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 @@ -531,7 +575,6 @@ class BurstEAP: Sentinel-1 IW SLC burst ''' # from LADS - num_sample: int # number of samples freq_sampling: float # range sampling rate eta_start: datetime.datetime tau_0: float # imageInformation/slantRangeTime @@ -569,7 +612,6 @@ def from_product_annotation_and_aux_cal(cls, product_annotation: ProductAnnotati ''' id_closest = closest_block_to_azimuth_time(product_annotation.antenna_pattern_azimuth_time, azimuth_time) - num_sample = product_annotation.number_of_samples freq_sampling = product_annotation.range_sampling_rate eta_start = azimuth_time tau_0 = product_annotation.slant_range_time @@ -580,11 +622,49 @@ def from_product_annotation_and_aux_cal(cls, product_annotation: ProductAnnotati ascending_node_time = product_annotation.ascending_node_time - return cls(num_sample, freq_sampling, eta_start, tau_0, tau_sub, theta_sub, + return cls(freq_sampling, eta_start, tau_0, tau_sub, theta_sub, azimuth_time, ascending_node_time, gain_eap, delta_theta) + def compute_eap_compensation_lut(self, num_sample): + '''Returns LUT for EAP compensation whose size is `num_sample`. + Based on ESA docuemnt : + "Impact of the Elevation Antenna Pattern Phase Compensation + on the Interferometric Phase Preservation" + + Document URL: + https://sentinel.esa.int/documents/247904/1653440/Sentinel-1-IPF_EAP_Phase_correction + + Parameter: + num_sample: int + Size of the output LUT + + Return: + ------- + gain_eap_interpolatd: EAP phase for the burst to be compensated + + ''' + + n_elt = len(self.gain_eap) + + theta_am = (np.arange(n_elt) - (n_elt - 1) / 2) * self.delta_theta + + delta_anx = self.eta_start - self.ascending_node_time + theta_offnadir = self._anx2roll(delta_anx.seconds + delta_anx.microseconds * 1.0e-6) + + theta_eap = theta_am + theta_offnadir + + tau = self.tau_0 + np.arange(num_sample) / self.freq_sampling + + theta = np.interp(tau, self.tau_sub, self.theta_sub) + + interpolator_gain = interp1d(theta_eap, self.gain_eap) + gain_eap_interpolated = interpolator_gain(theta) + + return gain_eap_interpolated + + def _anx2roll(self, delta_anx): ''' Returns the Platform nominal roll as function of elapsed time from @@ -599,7 +679,7 @@ def _anx2roll(self, delta_anx): Returns ------- - _: float + nominal_roll: float Estimated nominal roll (degrees) ''' # Estimate altitude based on time elapsed since ANX @@ -637,7 +717,7 @@ def _anx2height(self, delta_anx): Returns: -------- - _: float + h_t: float nominal height of the platform ''' diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 84dee9a1..18335186 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -8,7 +8,6 @@ import isce3 import numpy as np from osgeo import gdal -from scipy.interpolate import InterpolatedUnivariateSpline, interp1d from s1reader import s1_annotation @@ -625,77 +624,25 @@ def swath_name(self): @property def thermal_noise_lut(self): - '''Returns the burst-sized LUT for thermal noise correction - - Returns - ------- - arr_lut_total: np.array - 2d array containing thermal noise correction look up table values ''' - + 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.') - nrows, ncols = self.shape - - # Interpolate the range noise vector - rg_lut_interp_obj = InterpolatedUnivariateSpline(self.burst_noise.range_pixel, - self.burst_noise.range_lut, - k=1) - if self.burst_noise.azimuth_last_range_sample is not None: - vec_rg = np.arange(self.burst_noise.azimuth_last_range_sample + 1) - else: - vec_rg = np.arange(ncols) - rg_lut_interpolated = rg_lut_interp_obj(vec_rg).reshape((1, ncols)) - - - # Interpolate the azimuth noise vector - if (self.burst_noise.azimuth_line is None) or (self.burst_noise.azimuth_lut is None): - az_lut_interpolated = np.ones(nrows).reshape((nrows, 1)) - else: # IPF >= 2.90 - az_lut_interp_obj = InterpolatedUnivariateSpline(self.burst_noise.azimuth_line, - self.burst_noise.azimuth_lut, - k=1) - vec_az = np.arange(self.burst_noise.line_from, self.burst_noise.line_to + 1) - az_lut_interpolated = az_lut_interp_obj(vec_az).reshape((nrows, 1)) - - arr_lut_total = np.matmul(az_lut_interpolated, rg_lut_interpolated) - - return arr_lut_total + return self.burst_noise.compute_thermal_noise_lut(self.shape) @property def eap_compensation_lut(self): '''Returns LUT for EAP compensation. - Based on ESA docuemnt : - "Impact of the Elevation Antenna Pattern Phase Compensation - on the Interferometric Phase Preservation" - - Document URL: - https://sentinel.esa.int/documents/247904/1653440/Sentinel-1-IPF_EAP_Phase_correction - + Returns: ------- - gain_eap: EAP phase for the burst to be compensated + _: Interpolated EAP gain for the burst's lines ''' + if self.burst_eap is None: + raise ValueError('burst_eap is not defined for this burst.' + f' IPF version = {self.ipf_version}') - n_elt = len(self.burst_eap.gain_eap) - - theta_am = (np.arange(n_elt) - (n_elt - 1) / 2) * self.burst_eap.delta_theta - - delta_anx = self.burst_eap.eta_start-self.burst_eap.ascending_node_time - theta_offnadir = self.burst_eap._anx2roll(delta_anx.seconds + delta_anx.microseconds * 1.0e-6) - - theta_eap = theta_am + theta_offnadir - - tau = self.burst_eap.tau_0 + np.arange(self.burst_eap.num_sample) / self.burst_eap.freq_sampling - - theta = np.interp(tau, self.burst_eap.tau_sub, self.burst_eap.theta_sub) - - interpolator_gain = interp1d(theta_eap, self.burst_eap.gain_eap) - gain_eap_interpolated = interpolator_gain(theta) - phi_eap = np.angle(gain_eap_interpolated) - c_j = np.complex64(1.0j) - gain_eap = np.exp(c_j * phi_eap) - - return gain_eap + return self.burst_eap.compute_eap_compensation_lut(self.width) diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index 895dc5b3..f2247226 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -301,10 +301,11 @@ def get_path_aux_cal(directory_aux_cal: str, str_platform: str, instrument_cfg_i Returns: -------- path_aux_cal: AUX_CAL file that corresponds to the criteria provided + ''' list_aux_cal = glob.glob(f'{directory_aux_cal}/{str_platform}_AUX_CAL_V*.SAFE.zip') list_aux_cal.sort() - path_aux_cal = list_aux_cal[inst_config_id-1] + path_aux_cal = list_aux_cal[instrument_cfg_id-1] return path_aux_cal From 637a7d5472dbee124844696a1f1fdf931b42fac7 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Tue, 6 Sep 2022 15:41:22 -0700 Subject: [PATCH 37/49] Changing the length unit in BurstEAP._anx2roll() to [m] --- src/s1reader/s1_annotation.py | 10 +++++----- src/s1reader/s1_burst_slc.py | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 5f1ca8af..577166d5 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -685,17 +685,17 @@ def _anx2roll(self, delta_anx): # Estimate altitude based on time elapsed since ANX altitude = self._anx2height(delta_anx) - # Reference altitude (km) - href = 711.700 + # Reference altitude (m) + href = 711700.0 # Reference boresight at reference altitude (degrees) boresight_ref = 29.450 - # Partial derivative of roll vs altitude (degrees/km) - alpha_roll = 0.0566 + # Partial derivative of roll vs altitude (degrees/m) + alpha_roll = 0.0566 / 1000.0 # Estimate nominal roll i.e. theta off nadir (degrees) - nominal_roll = boresight_ref - alpha_roll * (altitude/1000.0 - href) + nominal_roll = boresight_ref - alpha_roll * (altitude - href) return nominal_roll diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 18335186..18214edf 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -635,7 +635,7 @@ def thermal_noise_lut(self): @property def eap_compensation_lut(self): '''Returns LUT for EAP compensation. - + Returns: ------- _: Interpolated EAP gain for the burst's lines From 9c13a08d45752cde3c95ae3d32e186832326d0cb Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Tue, 6 Sep 2022 15:46:28 -0700 Subject: [PATCH 38/49] change on variable name for clarity --- src/s1reader/s1_annotation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 577166d5..8d1236f3 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -15,7 +15,7 @@ from packaging import version #IPF version from which the NADS has azimuth noise vector annotation -ipf_version_az_noise_vector_available_from = version.parse('2.90') +min_ipf_version_az_noise_vector = version.parse('2.90') @dataclass class AnnotationBase: @@ -208,7 +208,7 @@ def from_et(cls,et_in: ET, ipf_version: version.Version, path_annotation: str): cls.xml_et = et_in cls.basename_annotation = os.path.basename(path_annotation) - if ipf_version < ipf_version_az_noise_vector_available_from: #legacy SAFE data + 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') @@ -447,7 +447,7 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: azimuth_last_azimuth_line = noise_annotation.az_last_azimuth_line azimuth_last_range_sample = noise_annotation.az_last_range_sample - if ipf_version >= ipf_version_az_noise_vector_available_from: + 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)) From 265a6c275a9011c5c56629de49bcee00b864b1ee Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Tue, 13 Sep 2022 10:26:56 -0700 Subject: [PATCH 39/49] Addressing PEP8 issue (R0201) --- src/s1reader/s1_annotation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 8d1236f3..490b3c55 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -699,8 +699,8 @@ def _anx2roll(self, delta_anx): return nominal_roll - - def _anx2height(self, delta_anx): + @classmethod + def _anx2height(cls, delta_anx): ''' Returns the platform nominal height as function of elapse time from ascending node crossing time (ANX). From 4aac1536101b93ef5abff24d9b522867c1bfac38 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Tue, 13 Sep 2022 10:38:09 -0700 Subject: [PATCH 40/49] Reverting the unit in `_anx2roll()` into the original documentataion; Comment cleanup --- src/s1reader/s1_annotation.py | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 490b3c55..02adf2c2 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -672,6 +672,9 @@ def _anx2roll(self, delta_anx): Code copied from ISCE2. + The units in this function is based on the reference documentation in the URL below: + https://sentinel.esa.int/documents/247904/1653440/Sentinel-1-IPF_EAP_Phase_correction + Parameters ---------- delta_anx: float @@ -685,17 +688,17 @@ def _anx2roll(self, delta_anx): # Estimate altitude based on time elapsed since ANX altitude = self._anx2height(delta_anx) - # Reference altitude (m) - href = 711700.0 + # Reference altitude (km) + href = 711.700 # Reference boresight at reference altitude (degrees) boresight_ref = 29.450 # Partial derivative of roll vs altitude (degrees/m) - alpha_roll = 0.0566 / 1000.0 + alpha_roll = 0.0566 # Estimate nominal roll i.e. theta off nadir (degrees) - nominal_roll = boresight_ref - alpha_roll * (altitude - href) + nominal_roll = boresight_ref - alpha_roll * (altitude/1000.0 - href) return nominal_roll @@ -722,22 +725,22 @@ def _anx2height(cls, delta_anx): ''' - ### Average height + # Average height (m) h_0 = 707714.8 #;m - #### Perturbation amplitudes - h = np.array([8351.5, 8947.0, 23.32, 11.74]) #;m + # Perturbation amplitudes (m) + h = np.array([8351.5, 8947.0, 23.32, 11.74]) - #### Perturbation phases - phi = np.array([3.1495, -1.5655 , -3.1297, 4.7222]) #;radians + # Perturbation phases (radians) + phi = np.array([3.1495, -1.5655 , -3.1297, 4.7222]) - ###O rbital time period in seconds + # Orbital time period (seconds) t_orb = (12*24*60*60) / 175. - ### Angular velocity + # Angular velocity (rad/sec) worb = 2*np.pi / t_orb - #### Evaluation of series + # Evaluation of series h_t = h_0 for i, h_i in enumerate(h): h_t += h_i * np.sin((i+1) * worb * delta_anx + phi[i]) From 1b7f805e09e574b347fe349fd2de5995203e4f5a Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 14 Sep 2022 07:14:49 -0700 Subject: [PATCH 41/49] Update src/s1reader/s1_annotation.py Clarification of the description Co-authored-by: Gustavo H. X. Shiroma <52007211+gshiroma@users.noreply.github.com> --- src/s1reader/s1_annotation.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 02adf2c2..7e6fa68a 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -14,7 +14,8 @@ from scipy.interpolate import InterpolatedUnivariateSpline, interp1d from packaging import version -#IPF version from which the NADS has azimuth noise vector annotation +# 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') @dataclass From 14e046c6739273e90802cd850b0883108a781edf Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 14 Sep 2022 07:15:22 -0700 Subject: [PATCH 42/49] Update src/s1reader/s1_annotation.py Clarification on docsctring Co-authored-by: Gustavo H. X. Shiroma <52007211+gshiroma@users.noreply.github.com> --- src/s1reader/s1_annotation.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 7e6fa68a..120710a0 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -148,7 +148,8 @@ class CalibrationAnnotation(AnnotationBase): @classmethod def from_et(cls, et_in: ET, path_annotation: str): - '''Extracts the list of calibration informaton from etree from CADS + '''Extracts the list of calibration informaton from etree from + the Calibration Annotation Data Set (CADS). Parameters: ----------- et_in: ET From 69072d94fe3f9e46cc3d8f232482b25e5cf5e475 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 14 Sep 2022 07:16:03 -0700 Subject: [PATCH 43/49] Update src/s1reader/s1_annotation.py fix on docstring Co-authored-by: Gustavo H. X. Shiroma <52007211+gshiroma@users.noreply.github.com> --- src/s1reader/s1_annotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 120710a0..0f3cdfaa 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -158,7 +158,7 @@ def from_et(cls, et_in: ET, path_annotation: str): Returns: -------- cls: CalibrationAnnotation - instance from CalibrationAnnotation initialized by the input parameter + Instance of CalibrationAnnotation initialized by the input parameter ''' cls.xml_et = et_in From de4096e0967e6cca6045218606b6b77b79ee9fc3 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 14 Sep 2022 07:17:04 -0700 Subject: [PATCH 44/49] Update src/s1reader/s1_annotation.py unnecessary comment removed Co-authored-by: Gustavo H. X. Shiroma <52007211+gshiroma@users.noreply.github.com> --- src/s1reader/s1_annotation.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 0f3cdfaa..0d980a27 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -177,7 +177,6 @@ def from_et(cls, et_in: ET, path_annotation: str): @dataclass class NoiseAnnotation(AnnotationBase): '''Reader for Noise Annotation Data Set (NADS) for IW SLC - # REF: .../isce2/components/isceobj/Sensor/GRD/Sentinel1.py # : ESA Documentation "Thermal Denoising of Products Generated by the S-1 IPF" ''' basename_annotation: str From 663c73514393c52a2206386e6c3331865e56e591 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 14 Sep 2022 09:13:22 -0700 Subject: [PATCH 45/49] Docstring clarification --- src/s1reader/s1_reader.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index f2247226..149485c1 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -310,7 +310,9 @@ def get_path_aux_cal(directory_aux_cal: str, str_platform: str, instrument_cfg_i return path_aux_cal def is_eap_correction_necessary(ipf_version: version.Version) -> SimpleNamespace : - '''Examines if what level of EAP correction is necessary, based on the IPF version. + ''' + Examines if what level of elevation antenna pattern (EAP) correction is necessary. + based on the IPF version. Based on the comment on PR: https://github.com/opera-adt/s1-reader/pull/48#discussion_r926138372 Parameter From 90aa6d758f0d1b46de2d1842099d536684c20c8b Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 14 Sep 2022 09:13:45 -0700 Subject: [PATCH 46/49] Addressing PEP8 C0301 --- src/s1reader/s1_annotation.py | 274 +++++++++++++++++++++++++--------- 1 file changed, 205 insertions(+), 69 deletions(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 0d980a27..545ab50a 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -47,6 +47,7 @@ def _parse_scalar(cls, path_field: str, str_type: str): val_out becomes np.array when str_type is vector* ''' + elem_field = cls.xml_et.find(path_field) if str_type == 'datetime': val_out = datetime.datetime.strptime(elem_field.text, '%Y-%m-%dT%H:%M:%S.%f') @@ -137,6 +138,7 @@ def _parse_vectorlist(cls, name_vector_list: str, name_vector: str, str_type: st @dataclass class CalibrationAnnotation(AnnotationBase): '''Reader for Calibration Annotation Data Set (CADS)''' + basename_annotation: str list_azimuth_time: np.ndarray list_line: list @@ -148,7 +150,8 @@ class CalibrationAnnotation(AnnotationBase): @classmethod def from_et(cls, et_in: ET, path_annotation: str): - '''Extracts the list of calibration informaton from etree from + ''' + Extracts the list of calibration informaton from etree from the Calibration Annotation Data Set (CADS). Parameters: ----------- @@ -162,23 +165,48 @@ def from_et(cls, et_in: ET, path_annotation: str): ''' cls.xml_et = et_in - cls.basename_annotation = os.path.basename(path_annotation) - cls.list_azimuth_time = cls._parse_vectorlist('calibrationVectorList', 'azimuthTime', 'datetime') - cls.list_line = cls._parse_vectorlist('calibrationVectorList', 'line', 'scalar_int') - cls.list_pixel = cls._parse_vectorlist('calibrationVectorList', 'pixel', 'vector_int') - cls.list_sigma_nought = cls._parse_vectorlist('calibrationVectorList', 'sigmaNought', 'vector_float') - cls.list_beta_nought = cls._parse_vectorlist('calibrationVectorList', 'betaNought', 'vector_float') - cls.list_gamma = cls._parse_vectorlist('calibrationVectorList', 'gamma', 'vector_float') - cls.list_dn = cls._parse_vectorlist('calibrationVectorList', 'dn', 'vector_float') + cls.basename_annotation = \ + os.path.basename(path_annotation) + + cls.list_azimuth_time = \ + cls._parse_vectorlist('calibrationVectorList', + 'azimuthTime', + 'datetime') + cls.list_line = \ + cls._parse_vectorlist('calibrationVectorList', + 'line', + 'scalar_int') + cls.list_pixel = \ + cls._parse_vectorlist('calibrationVectorList', + 'pixel', + 'vector_int') + cls.list_sigma_nought = \ + cls._parse_vectorlist('calibrationVectorList', + 'sigmaNought', + 'vector_float') + cls.list_beta_nought = \ + cls._parse_vectorlist('calibrationVectorList', + 'betaNought', + 'vector_float') + cls.list_gamma = \ + cls._parse_vectorlist('calibrationVectorList', + 'gamma', + 'vector_float') + cls.list_dn = \ + cls._parse_vectorlist('calibrationVectorList', + 'dn', + 'vector_float') return cls @dataclass class NoiseAnnotation(AnnotationBase): - '''Reader for Noise Annotation Data Set (NADS) for IW SLC - # : ESA Documentation "Thermal Denoising of Products Generated by the S-1 IPF" ''' + 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 @@ -193,7 +221,8 @@ class NoiseAnnotation(AnnotationBase): @classmethod def from_et(cls,et_in: ET, ipf_version: version.Version, path_annotation: str): - '''Extracts list of noise information from etree + ''' + Extracts list of noise information from etree Parameter ---------- @@ -209,11 +238,24 @@ def from_et(cls,et_in: ET, ipf_version: version.Version, path_annotation: str): 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') + 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 @@ -222,29 +264,63 @@ def from_et(cls,et_in: ET, ipf_version: version.Version, path_annotation: str): 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] + 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): - '''For L1 SLC product annotation file. For EAP correction.''' + ''' + Reader for L1 Product annotation for IW SLC + For Elevation Antenna Pattern (EAP) correction + ''' + image_information_slant_range_time: float # Attributes to be used when determining what AUX_CAL to load instrument_cfg_id: int - #elevation_angle: + # elevation_angle: antenna_pattern_azimuth_time: list antenna_pattern_slant_range_time: list antenna_pattern_elevation_angle: list @@ -258,8 +334,10 @@ class ProductAnnotation(AnnotationBase): @classmethod def from_et(cls, et_in: ET): - '''Extracts list of product information from etree from L1 annotation data set (LADS) - Parameter + ''' + Extracts list of product information from etree from + L1 annotation data set (LADS) Parameter + ---------- et_in : xml.etree.ElementTree Parsed LADS annotation .xml @@ -272,25 +350,54 @@ def from_et(cls, et_in: ET): cls.xml_et = et_in - cls.image_information_slant_range_time = cls._parse_scalar('imageAnnotation/imageInformation/slantRangeTime', 'scalar_float') - cls.inst_config_id = cls._parse_scalar('generalAnnotation/downlinkInformationList/downlinkInformation/downlinkValues/instrumentConfigId', 'scalar_int') - cls.antenna_pattern_azimuth_time = cls._parse_vectorlist('antennaPattern/antennaPatternList', 'azimuthTime', 'datetime') - cls.antenna_pattern_slant_range_time = cls._parse_vectorlist('antennaPattern/antennaPatternList', 'slantRangeTime', 'vector_float') - cls.antenna_pattern_elevation_angle = cls._parse_vectorlist('antennaPattern/antennaPatternList', 'elevationAngle', 'vector_float') - cls.antenna_pattern_elevation_pattern = cls._parse_vectorlist('antennaPattern/antennaPatternList', 'elevationPattern', 'vector_float') - cls.ascending_node_time = cls._parse_scalar('imageAnnotation/imageInformation/ascendingNodeTime', 'datetime') - cls.number_of_samples = cls._parse_scalar('imageAnnotation/imageInformation/numberOfSamples', 'scalar_int') - cls.number_of_samples = cls._parse_scalar('imageAnnotation/imageInformation/numberOfSamples', 'scalar_int') - cls.range_sampling_rate = cls._parse_scalar('generalAnnotation/productInformation/rangeSamplingRate', 'scalar_float') - - cls.slant_range_time = cls._parse_scalar('imageAnnotation/imageInformation/slantRangeTime','scalar_float') + cls.antenna_pattern_azimuth_time = \ + cls._parse_vectorlist('antennaPattern/antennaPatternList', + 'azimuthTime', + 'datetime') + cls.antenna_pattern_slant_range_time = \ + cls._parse_vectorlist('antennaPattern/antennaPatternList', + 'slantRangeTime', + 'vector_float') + cls.antenna_pattern_elevation_angle = \ + cls._parse_vectorlist('antennaPattern/antennaPatternList', + 'elevationAngle', + 'vector_float') + cls.antenna_pattern_elevation_pattern = \ + cls._parse_vectorlist('antennaPattern/antennaPatternList', + 'elevationPattern', + 'vector_float') + + cls.image_information_slant_range_time = \ + cls._parse_scalar('imageAnnotation/imageInformation/slantRangeTime', + 'scalar_float') + cls.ascending_node_time = \ + cls._parse_scalar('imageAnnotation/imageInformation/ascendingNodeTime', + 'datetime') + cls.number_of_samples = \ + cls._parse_scalar('imageAnnotation/imageInformation/numberOfSamples', + 'scalar_int') + cls.number_of_samples = \ + cls._parse_scalar('imageAnnotation/imageInformation/numberOfSamples', + 'scalar_int') + cls.range_sampling_rate = \ + cls._parse_scalar('generalAnnotation/productInformation/rangeSamplingRate', + 'scalar_float') + cls.slant_range_time = \ + cls._parse_scalar('imageAnnotation/imageInformation/slantRangeTime', + 'scalar_float') + + cls.inst_config_id = \ + cls._parse_scalar('generalAnnotation/downlinkInformationList/downlinkInformation/' + 'downlinkValues/instrumentConfigId', + 'scalar_int') return cls @dataclass class AuxCal(AnnotationBase): - '''AUX_CAL information for elevation antenna pattern (EAP) correction''' + '''AUX_CAL information for EAP correction''' + beam_nominal_near_range: float beam_nominal_far_range: float @@ -307,7 +414,8 @@ class AuxCal(AnnotationBase): @classmethod def load_from_zip_file(cls, path_aux_cal_zip: str, pol: str, str_swath: str): - '''A class method that Extracts list of information AUX_CAL from the input ET. + ''' + A class method that extracts list of information AUX_CAL from the input ET. Parameters --------- @@ -337,29 +445,50 @@ def load_from_zip_file(cls, path_aux_cal_zip: str, pol: str, str_swath: str): swath_xml = calibration_params.find('swath').text polarisation_xml = calibration_params.find('polarisation').text if polarisation_xml == pol.upper() and swath_xml==str_swath.upper(): - print(f'Found a calibration parameters for swath {str_swath} and polarization {pol}.') - cls.beam_nominal_near_range = float(calibration_params.find('elevationAntennaPattern/beamNominalNearRange').text) - cls.beam_nominal_far_range = float(calibration_params.find('elevationAntennaPattern/beamNominalFarRange').text) - cls.elevation_angle_increment = float(calibration_params.find('elevationAntennaPattern/elevationAngleIncrement').text) - - n_val = int(calibration_params.find('elevationAntennaPattern/values').attrib['count']) - arr_eap_val = np.array([float(token_val) for token_val in calibration_params.find('elevationAntennaPattern/values').text.split()]) - if n_val == len(arr_eap_val): #Provided in real numbers: In case of AUX_CAL for old IPFs. + cls.beam_nominal_near_range = \ + float(calibration_params. + find('elevationAntennaPattern/beamNominalNearRange').text) + cls.beam_nominal_far_range = \ + float(calibration_params. + find('elevationAntennaPattern/beamNominalFarRange').text) + cls.elevation_angle_increment = \ + float(calibration_params. + find('elevationAntennaPattern/elevationAngleIncrement').text) + + n_val = \ + int(calibration_params. + find('elevationAntennaPattern/values').attrib['count']) + arr_eap_val = \ + np.array([float(token_val) for \ + token_val in calibration_params. + find('elevationAntennaPattern/values').text.split()]) + + if n_val == len(arr_eap_val): + # Provided in real numbers: In case of AUX_CAL for old IPFs. cls.elevation_antenna_pattern = arr_eap_val - elif n_val*2 == len(arr_eap_val): #Provided in complex numbers: In case of recent IPFs e.g. 3.10 + elif n_val*2 == len(arr_eap_val): + # Provided in complex numbers: In case of recent IPFs e.g. 3.10 cls.elevation_antenna_pattern = arr_eap_val[0::2] + arr_eap_val[1::2] * 1.0j else: - raise ValueError(f'The number of values does not match. n_val={n_val}, #values in elevationAntennaPattern/values={len(arr_eap_val)}') + raise ValueError('The number of values does not match. ' + f'n_val={n_val}, ' + f'#len(elevationAntennaPattern/values)={len(arr_eap_val)}') cls.azimuth_angle_increment = \ - float(calibration_params.find('azimuthAntennaPattern/azimuthAngleIncrement').text) + float(calibration_params. + find('azimuthAntennaPattern/azimuthAngleIncrement').text) cls.azimuth_antenna_pattern = \ - np.array([float(token_val) for token_val in calibration_params.find('azimuthAntennaPattern/values').text.split()]) + np.array([float(token_val) for \ + token_val in calibration_params. + find('azimuthAntennaPattern/values').text.split()]) cls.azimuth_antenna_element_pattern_increment = \ - float(calibration_params.find('azimuthAntennaElementPattern/azimuthAngleIncrement').text) + float(calibration_params. + find('azimuthAntennaElementPattern/azimuthAngleIncrement').text) cls.azimuth_antenna_element_pattern = \ - np.array([float(token_val) for token_val in calibration_params.find('azimuthAntennaElementPattern/values').text.split()]) + np.array([float(token_val) for \ + token_val in calibration_params. + find('azimuthAntennaElementPattern/values').text.split()]) cls.absolute_calibration_constant = \ float(calibration_params.find('absoluteCalibrationConstant').text) @@ -370,7 +499,8 @@ def load_from_zip_file(cls, path_aux_cal_zip: str, pol: str, str_swath: str): def closest_block_to_azimuth_time(vector_azimuth_time: np.ndarray, azimuth_time_burst: datetime.datetime) -> int: - '''Find the id of the closest data block in annotation. + ''' + Find the id of the closest data block in annotation. To be used when populating BurstNoise, BurstCalibration, and BurstEAP. Parameters @@ -410,9 +540,13 @@ class BurstNoise: 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 + 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 ---------- @@ -430,7 +564,7 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: Returns ------- cls: BurstNoise - An instance from BurstNoise initialized by the input parameters + Instance of BurstNoise initialized by the input parameters ''' @@ -473,7 +607,8 @@ def from_noise_annotation(cls, noise_annotation: NoiseAnnotation, azimuth_time: def compute_thermal_noise_lut(self, shape_lut): - '''Calculate thermal noise LUT whose shape is `shape_lut` + ''' + Calculate thermal noise LUT whose shape is `shape_lut` Parameter: ---------- @@ -544,7 +679,7 @@ def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotati Returns ------- cls: BurstCalibration - EAP correction info for the burst + Radiometric correction information for the burst ''' basename_cads = calibration_annotation.basename_annotation @@ -572,8 +707,7 @@ def from_calibration_annotation(cls, calibration_annotation: CalibrationAnnotati @dataclass class BurstEAP: - '''Elevation antenna pattern (EAP) correction information for - Sentinel-1 IW SLC burst + '''EAP correction information for Sentinel-1 IW SLC burst ''' # from LADS freq_sampling: float # range sampling rate @@ -629,7 +763,8 @@ def from_product_annotation_and_aux_cal(cls, product_annotation: ProductAnnotati def compute_eap_compensation_lut(self, num_sample): - '''Returns LUT for EAP compensation whose size is `num_sample`. + ''' + Returns LUT for EAP compensation whose size is `num_sample`. Based on ESA docuemnt : "Impact of the Elevation Antenna Pattern Phase Compensation on the Interferometric Phase Preservation" @@ -686,6 +821,7 @@ def _anx2roll(self, delta_anx): nominal_roll: float Estimated nominal roll (degrees) ''' + # Estimate altitude based on time elapsed since ANX altitude = self._anx2height(delta_anx) From 2ce289a9ac26290d50f156ac9e033f6e6e73b20b Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 14 Sep 2022 15:49:52 -0700 Subject: [PATCH 47/49] AuxCal.get_aux_cal_instrument_config_id() - check the instrument configuration id of the AUX_CAL .zip --- src/s1reader/s1_annotation.py | 38 ++++++++++++++++++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 545ab50a..884dac6e 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -4,6 +4,7 @@ ''' from dataclasses import dataclass + import datetime import os import xml.etree.ElementTree as ET @@ -411,6 +412,38 @@ class AuxCal(AnnotationBase): absolute_calibration_constant: float noise_calibration_factor: float + @classmethod + def get_aux_cal_instrument_config_id(cls, path_aux_cal_zip: str): + ''' + Returns the instrument configuration ID of the input AUX_CAL .zip file + + Parameter: + ---------- + path_aux_cal_zip: str + Path to the input aux_cal .zip file + + Return: + ------- + instrument_config_id_aux_cal: int + instrument configuration ID of the input AUX_CAL file + + ''' + + # Check the instrument ID in the manifest. + s1auxsar = '{http://www.esa.int/safe/sentinel-1.0/sentinel-1/auxiliary/sar}' + str_safe_aux_cal = os.path.basename(path_aux_cal_zip).replace('.zip','') + + path_config_id = ('metadataSection/metadataObject/metadataWrap/xmlData' + f'/{s1auxsar}standAloneProductInformation' + f'/{s1auxsar}instrumentConfigurationId') + + with zipfile.ZipFile(path_aux_cal_zip, 'r') as zipfile_aux_cal: + with zipfile_aux_cal.open(f'{str_safe_aux_cal}/manifest.safe','r') as f_manifest: + et_manifest = ET.parse(f_manifest) + instrument_config_id_aux_cal = int(et_manifest.find(path_config_id).text) + + return instrument_config_id_aux_cal + @classmethod def load_from_zip_file(cls, path_aux_cal_zip: str, pol: str, str_swath: str): @@ -431,6 +464,10 @@ def load_from_zip_file(cls, path_aux_cal_zip: str, pol: str, str_swath: str): cls: AuxCal class populated by et_in in the parameter ''' + + if not path_aux_cal_zip.endswith('.zip'): + raise ValueError('Only AUX_CAL files in .zip format are accepted.') + if os.path.exists(path_aux_cal_zip): str_safe_aux_cal = os.path.basename(path_aux_cal_zip).replace('.zip','') else: @@ -649,7 +686,6 @@ def compute_thermal_noise_lut(self, shape_lut): return arr_lut_total - @dataclass class BurstCalibration: '''Calibration information for Sentinel-1 IW SLC burst From 2913b291b9665a2dd24576b67652ef52b1858b1f Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Wed, 14 Sep 2022 15:50:41 -0700 Subject: [PATCH 48/49] get_path_aux_cal() - check if the instrument configuration ID matches --- src/s1reader/s1_reader.py | 27 ++++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index 149485c1..97ad76b1 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -296,17 +296,30 @@ def get_path_aux_cal(directory_aux_cal: str, str_platform: str, instrument_cfg_i str_platform: str 'S1A' or 'S1B' instrument_cfg_id: int - Instrument configuration ID + Instrument configuration ID. - Returns: + Return: -------- - path_aux_cal: AUX_CAL file that corresponds to the criteria provided + path_aux_cal: str + Path to the AUX_CAL file that corresponds to the criteria provided + None if the matching AUX_CAL is not found in `directory_aux_cal` ''' list_aux_cal = glob.glob(f'{directory_aux_cal}/{str_platform}_AUX_CAL_V*.SAFE.zip') list_aux_cal.sort() + + # Initial guess of the matching aux_cal path_aux_cal = list_aux_cal[instrument_cfg_id-1] + if instrument_cfg_id != AuxCal.get_aux_cal_instrument_config_id(path_aux_cal): + # Failed to guess the instrument ID based on the files' orders in the list. + # Go through the all AUX_CAL files in the list to find the matching file + path_aux_cal = None + for path_aux_cal in list_aux_cal: + if instrument_cfg_id == AuxCal.get_aux_cal_instrument_config_id(path_aux_cal): + return path_aux_cal + + return path_aux_cal def is_eap_correction_necessary(ipf_version: version.Version) -> SimpleNamespace : @@ -412,11 +425,11 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, flag_apply_eap = is_eap_correction_necessary(ipf_version) if flag_apply_eap.phase_correction: path_aux_cals = f'{os.path.dirname(s1_annotation.__file__)}/data/aux_cal' - path_aux_cal_zip = get_path_aux_cal(path_aux_cals, - platform_id, - product_annotation.inst_config_id) + path_aux_cal = get_path_aux_cal(path_aux_cals, + platform_id, + product_annotation.inst_config_id) - aux_cal_subswath = AuxCal.load_from_zip_file(path_aux_cal_zip, + aux_cal_subswath = AuxCal.load_from_zip_file(path_aux_cal, pol, subswath_id) else: From f019bae3139be725b7d46f46e7eb2bd1fdd6b8a5 Mon Sep 17 00:00:00 2001 From: Liang Yu Date: Wed, 21 Sep 2022 11:36:05 -0700 Subject: [PATCH 49/49] untested refactor to isolate noise loading from annotation --- src/s1reader/s1_annotation.py | 237 --------------------------------- src/s1reader/s1_noise.py | 156 ++++++++++++++++++++++ src/s1reader/s1_reader.py | 41 ++---- src/s1reader/utils/__init__.py | 0 src/s1reader/utils/utils.py | 25 ++++ 5 files changed, 192 insertions(+), 267 deletions(-) create mode 100644 src/s1reader/s1_noise.py create mode 100644 src/s1reader/utils/__init__.py create mode 100644 src/s1reader/utils/utils.py diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 884dac6e..6fe39d6b 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): ''' @@ -557,135 +449,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_noise.py b/src/s1reader/s1_noise.py new file mode 100644 index 00000000..75bf3896 --- /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 97ad76b1..bfc419c9 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. @@ -416,10 +396,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 flag_apply_eap = is_eap_correction_necessary(ipf_version) @@ -557,8 +535,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 flag_apply_eap.phase_correction: burst_aux_cal = BurstEAP.from_product_annotation_and_aux_cal(product_annotation, aux_cal_subswath, diff --git a/src/s1reader/utils/__init__.py b/src/s1reader/utils/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/src/s1reader/utils/utils.py b/src/s1reader/utils/utils.py new file mode 100644 index 00000000..63160c60 --- /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)