diff --git a/src/s1reader/s1_annotation.py b/src/s1reader/s1_annotation.py index 38085c0..21d3981 100644 --- a/src/s1reader/s1_annotation.py +++ b/src/s1reader/s1_annotation.py @@ -15,6 +15,8 @@ from scipy.interpolate import InterpolatedUnivariateSpline, interp1d from packaging import version +from isce3.core import speed_of_light + # 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') @@ -333,6 +335,7 @@ class ProductAnnotation(AnnotationBase): slant_range_time: float + @classmethod def from_et(cls, et_in: ET): ''' @@ -897,3 +900,120 @@ def _anx2height(cls, delta_anx): h_t += h_i * np.sin((i+1) * worb * delta_anx + phi[i]) return h_t + + +@dataclass +class BurstExtendedCoeffs: + ''' + Segments of FM rate / Doppler centroid polynomial coefficients. + For (linear) interpolation of FM rate / Doppler Centroid along azimuth. + To be used for calculating azimuth FM rate mismatch mitigation + ''' + + # FM rate + fm_rate_aztime_vec: np.ndarray + fm_rate_coeff_arr: np.ndarray + fm_rate_tau0_vec:np.ndarray + + # Doppler centroid + dc_aztime_vec: np.ndarray + dc_coeff_arr: np.ndarray + dc_tau0_vec: np.ndarray + + @classmethod + def from_polynomial_lists(cls, + az_fm_rate_list: list, + doppler_centroid_list: list, + sensing_start: datetime.datetime, + sensing_end: datetime.datetime): + ''' + Extract coefficients from the list of the polynomial lists that fall within + the provided sensing start / end times of a burst. + + Parameters: + ----------- + az_fm_rate_list: list[isce3.core.Poly1d] + List of azimuth FM rate polynomials + doppler_centroid_list: list[isce3.core.Poly1d] + List of doppler centroid polynomials + sensing_start: datetime.datetime + Azimuth start time of the burst + sensing_end: datetime.datetime + Azimuth end time of the burst + ''' + + # Extract polynomial info for azimuth FM rate + (fm_rate_aztime_burst_vec, + fm_rate_coeff_burst_arr, + fm_rate_tau0_burst_vec) = cls.extract_polynomial_sequence(az_fm_rate_list, + sensing_start, + sensing_end) + + (dc_aztime_burst_vec, + dc_coeff_burst_arr, + dc_tau0_burst_vec) = cls.extract_polynomial_sequence(doppler_centroid_list, + sensing_start, + sensing_end) + + return cls(fm_rate_aztime_burst_vec, fm_rate_coeff_burst_arr, fm_rate_tau0_burst_vec, + dc_aztime_burst_vec, dc_coeff_burst_arr, dc_tau0_burst_vec) + + + @classmethod + def extract_polynomial_sequence(cls, polynomial_list: list, + datetime_start: datetime.datetime, + datetime_end: datetime.datetime): + ''' + Scan `vec_azimuth_time` end find indices of the vector + that covers the period defined with + `datetime_start` and `datetime_end` + + Parameters: + ----------- + polynomial_list: list + list of (azimuth_time, isce3.core.Poly1d) + datetime_start: datetime.datetime + Start time of the period + datetime_end: datetime.datetime + end time of the period + + Returns: + -------- + tuple + Tuple of (vec_aztime_sequence, + arr_coeff_sequence, + vec_tau0_sequence) + as a sequence of polynomial info that covers the period + defined in the parameters. + vec_aztime_sequence: azimuth time of each sample in the sequence + arr_coeff_sequence: N by 3 npy array whose row is coefficients of + each sample in the sequence + vec_tau0_sequence: Range start time of each sample in the sequence + ''' + + # NOTE: dt is defined as: [azimuth time] - [start/end time] + # find index of poly time closest to start time that is less than start time + dt_wrt_start = np.array([(poly[0] - datetime_start).total_seconds() for poly in polynomial_list]) + dt_wrt_start = np.ma.masked_array(dt_wrt_start, mask=dt_wrt_start > 0) + index_start = np.argmax(dt_wrt_start) + + # find index of poly time closest to end time that is greater than end time + dt_wrt_end = np.array([(poly[0] - datetime_end).total_seconds() for poly in polynomial_list]) + dt_wrt_end = np.ma.masked_array(dt_wrt_end, mask=dt_wrt_end < 0) + index_end = np.argmin(dt_wrt_end) + + # Done extracting the IDs. Extract the polynomial sequence + vec_aztime_sequence = [] + arr_coeff_sequence = [] + vec_tau0_sequence = [] + + # Scale factor to convert range (in meters) to seconds (tau) + range_to_tau = 2.0 / speed_of_light + for poly in polynomial_list[index_start:index_end+1]: + vec_aztime_sequence.append(poly[0]) + arr_coeff_sequence.append(poly[1].coeffs) + vec_tau0_sequence.append(poly[1].mean * range_to_tau) + + return (np.array(vec_aztime_sequence), + np.array(arr_coeff_sequence), + np.array(vec_tau0_sequence)) diff --git a/src/s1reader/s1_burst_slc.py b/src/s1reader/s1_burst_slc.py index 022df0a..a506806 100644 --- a/src/s1reader/s1_burst_slc.py +++ b/src/s1reader/s1_burst_slc.py @@ -8,7 +8,7 @@ import isce3 import numpy as np from osgeo import gdal - +from scipy.interpolate import InterpolatedUnivariateSpline from s1reader import s1_annotation from .s1_burst_id import S1BurstId @@ -106,6 +106,86 @@ class represents a polynomial function of range poly = isce3.core.Poly2d(coeffs, xmin, ymin, xnorm, ynorm) return poly +def _evaluate_polynomial_array(arr_polynomial, grid_tau, vec_tau0): + ''' + Evaluate the polynomials on the correction grid. + To be used for calculating azimuth FM mismatch rate. + + Parameters: + ----------- + arr_polynomial: np.ndarray + coefficients interpolated to azimuth times in each line + in correction grid + grid_tau: np.ndarray + 2d numpy array filled with range time in the correction grid + vec_tau0: np.ndarray + range start time for each line in the correction grid + + Return: + ------- + eval_out: np.ndarray + Evaluated values on the correction grid + ''' + + ncol = grid_tau.shape[1] + term_tau = grid_tau - vec_tau0 * np.ones(ncol)[np.newaxis, ...] + + eval_out = ( + arr_polynomial[:, 0][..., np.newaxis] * np.ones(ncol)[np.newaxis, ...] + + (arr_polynomial[:, 1][..., np.newaxis] * np.ones(ncol)[np.newaxis, ...]) + * term_tau + + (arr_polynomial[:, 2][..., np.newaxis] * np.ones(ncol)[np.newaxis, ...]) + * term_tau**2 + ) + + return eval_out + +def _llh_to_ecef(lat, lon, hgt, ellipsoid, in_degree=True): + ''' + Calculate cartesian coordinates in ECEF from + latitude, longitude, and altitude + + Parameters: + ----------- + lat: np.ndarray + latitude as numpy array + lon: np.ndarray + longitude as numpy array + hgt: np.ndarray + height as numpy array + ellipsoid: isce3.core.Ellipsoid + Definition of Ellipsoid + in_degree: bool + True if the units of lat and lon are degrees. + False if the units are radian. + + Return: + x_ecef : ndarray + ECEF X coordinate (in meters). + y_ecef : ndarray + ECEF X coordinate (in meters). + z_ecef : ndarray + ECEF X coordinate (in meters). + x, y, and z as a tuple of np.ndarray + + ''' + + if in_degree: + rad_lat = np.radians(lat) + rad_lon = np.radians(lon) + else: + rad_lat = lat + rad_lon = lon + + v_ellipsoid = ellipsoid.a / np.sqrt(1 - ellipsoid.e2 * np.sin(rad_lat) * np.sin(rad_lat)) + + x_ecef = (v_ellipsoid + hgt) * np.cos(rad_lat) * np.cos(rad_lon) + y_ecef = (v_ellipsoid + hgt) * np.cos(rad_lat) * np.sin(rad_lon) + z_ecef = (v_ellipsoid * (1 - ellipsoid.e2) + hgt) * np.sin(rad_lat) + + return (x_ecef, y_ecef, z_ecef) + + @dataclass class AzimuthCarrierComponents: kt: np.ndarray @@ -173,6 +253,10 @@ class Sentinel1BurstSlc: burst_noise: s1_annotation.BurstNoise # Thermal noise correction burst_eap: s1_annotation.BurstEAP # EAP correction + # Time series of FM rate / Doppler centroid polynomial coefficients + # for azimuth FM rate mismatch mitigation + extended_coeffs: s1_annotation.BurstExtendedCoeffs + def __str__(self): return f"Sentinel1BurstSlc: {self.burst_id} at {self.sensing_start}" @@ -621,6 +705,281 @@ def az_carrier_components(self, offset, position): return AzimuthCarrierComponents(kt, eta, eta_ref) + def az_fm_rate_mismatch_mitigation(self, + path_dem: str, + path_scratch: str = None, + range_step=None, + az_step=None, + threshold_rdr2geo=1e-8, + numiter_rdr2geo=25, + custom_radargrid=None): + ''' + Calculate azimuth FM rate mismatch mitigation + Based on ETAD-DLR-DD-0008, Algorithm Technical Baseline Document. + Available: https://sentinels.copernicus.eu/documents/247904/4629150/ + Sentinel-1-ETAD-Algorithm-Technical-Baseline-Document.pdf + + Parameters: + ----------- + path_dem: str + Path to the DEM to calculate the actual azimuth FM rate + path_scratch: str + Path to the scratch directory to store intermediate data + If None, the scratch files will be saved on temporary directory generated by system + range_step: float + Range step of the correction grid in meters + az_step: float + Azimuth step of the correction grid in seconds + threshold_rdr2geo: int + Threshold of the iteration for rdr2geo + numiter_rdr2geo: int + Maximum number of iteration for rdr2geo + custom_radargrid: isce3.product.RadarGridParameters + ISCE3 radar grid to define the correction grid. + If None, the full resolution radargrid of the burst will be used. + + Return: + ------- + _: isce3.core.LUT2d + azimuth FM rate mismatch rate in radar grid in seconds. + Examples + ---------- + >>> correction_grid = burst.az_fm_rate_mismatch_mitigation("my_dem.tif") + ''' + + # Create temporary directory for scratch when + # `path_scratch` is not provided. + if path_scratch is None: + temp_dir_obj = tempfile.TemporaryDirectory() + path_scratch = temp_dir_obj.name + else: + temp_dir_obj = None + + + os.makedirs(path_scratch, exist_ok=True) + + correction_radargrid = self.as_isce3_radargrid() + if custom_radargrid is not None: + correction_radargrid = custom_radargrid + + # Override the radargrid definition if `rg_step` and `az_step` is defined + if range_step and az_step: + if custom_radargrid is not None: + warnings.warn('range_step and az_step assigned. ' + 'Overriding the custom radargrid definition.') + + width_radargrid, length_radargrid = \ + [vec.size for vec in self._steps_to_vecs(range_step, az_step)] + sensing_start_radargrid = self.as_isce3_radargrid().sensing_start + correction_radargrid = isce3.product.RadarGridParameters( + sensing_start_radargrid, + self.wavelength, + 1/az_step, + self.starting_range, + range_step, + isce3.core.LookSide.Right, + length_radargrid, + width_radargrid, + self.as_isce3_radargrid().ref_epoch) + + # Define the correction grid from the radargrid + # Also define the staggered grid in azimuth to calculate acceeration + width_grid = correction_radargrid.width + length_grid = correction_radargrid.length + intv_t = 1.0 / correction_radargrid.prf + intv_tau = (correction_radargrid.range_pixel_spacing * 2.0 + / isce3.core.speed_of_light) + delta_sec_from_ref_epoch = ((correction_radargrid.ref_epoch + - self.orbit.reference_epoch).total_seconds() + + correction_radargrid.sensing_start) + vec_t = np.arange(length_grid) * intv_t + delta_sec_from_ref_epoch + vec_t_staggered = (np.arange(length_grid + 1) + * intv_t + + delta_sec_from_ref_epoch + - intv_t / 2) + + vec_tau = np.arange(width_grid)*intv_tau + + grid_tau, grid_t = np.meshgrid(vec_tau, vec_t) + + vec_position_intp = np.zeros((length_grid, 3)) + vec_vel_intp = np.zeros((length_grid, 3)) + vec_vel_intp_staggered = [None] * (length_grid + 1) + for i_azimuth, t_azimuth in enumerate(vec_t): + vec_position_intp[i_azimuth, :], vec_vel_intp[i_azimuth, :] = \ + self.orbit.interpolate(t_azimuth) + + # Calculate velocity on staggered grid + # to kinematically calculate acceleration + for i_azimuth, t_azimuth in enumerate(vec_t_staggered): + _, vec_vel_intp_staggered[i_azimuth] = self.orbit.interpolate(t_azimuth) + + vec_acceleration_intp = np.diff(vec_vel_intp_staggered, axis=0) / intv_t + + # convert azimuth time to seconds from the reference epoch of burst_in.orbit + fm_rate_aztime_sec_vec = [(isce3.core.DateTime(datetime_vec) \ + - self.orbit.reference_epoch).total_seconds() + for datetime_vec in self.extended_coeffs.fm_rate_aztime_vec] + + dc_aztime_sec_vec = [(isce3.core.DateTime(datetime_vec) \ + - self.orbit.reference_epoch).total_seconds() + for datetime_vec in self.extended_coeffs.dc_aztime_vec] + + # calculate splined interpolation of the coeffs. and tau_0s + interpolator_tau0_ka = InterpolatedUnivariateSpline( + fm_rate_aztime_sec_vec, + self.extended_coeffs.fm_rate_tau0_vec, + k=1) + tau0_ka_interp = interpolator_tau0_ka(vec_t)[..., np.newaxis] + + interpolator_tau0_fdc_interp = InterpolatedUnivariateSpline( + dc_aztime_sec_vec, + self.extended_coeffs.dc_tau0_vec, + k=1) + tau0_fdc_interp = interpolator_tau0_fdc_interp(vec_t)[..., np.newaxis] + + # add range time origin to vec_tau + grid_tau += tau0_ka_interp * np.ones(vec_tau.shape)[np.newaxis, ...] + + # Interpolate the DC and fm rate coeffs along azimuth time + def interp_coeffs(az_time, coeffs, az_time_interp): + '''Convenience function to interpolate DC and FM rate coeffiicients + ''' + interpolated_coeffs = [] + for i in range(3): + coeff_interpolator = InterpolatedUnivariateSpline( + az_time, coeffs[:, i], k=1) + interpolated_coeffs.append(coeff_interpolator(az_time_interp)) + return np.array(interpolated_coeffs).transpose() + + dc_coeffs = interp_coeffs(fm_rate_aztime_sec_vec, + self.extended_coeffs.dc_coeff_arr, + vec_t) + + fm_rate_coeffs = interp_coeffs(fm_rate_aztime_sec_vec, + self.extended_coeffs.fm_rate_coeff_arr, + vec_t) + + # Run topo on scratch directory + dem_raster = isce3.io.Raster(path_dem) + epsg = dem_raster.get_epsg() + proj = isce3.core.make_projection(epsg) + ellipsoid = proj.ellipsoid + grid_doppler = isce3.core.LUT2d() + + rdr2geo_obj = isce3.geometry.Rdr2Geo( + correction_radargrid, + self.orbit, + ellipsoid, + grid_doppler, + threshold=threshold_rdr2geo, + numiter=numiter_rdr2geo) + + str_datetime = datetime.datetime.now().strftime('%Y%m%d_%H%M%S%f') + list_filename_llh = [f'{path_scratch}/{llh}_{str_datetime}.rdr' + for llh + in ['lat', 'lon', 'hgt']] + lat_raster, lon_raster, hgt_raster = [isce3.io.Raster( + filename_llh, + correction_radargrid.width, + correction_radargrid.length, + 1, + gdal.GDT_Float64, + 'ENVI') + for filename_llh in list_filename_llh] + + rdr2geo_obj.topo(dem_raster, lon_raster, lat_raster, hgt_raster) + + # make sure that the ISCE3 rasters are written out to file system + lat_raster.close_dataset() + lon_raster.close_dataset() + hgt_raster.close_dataset() + + # Do the bunch of computation + kappa_steer_vec = (2 * (np.linalg.norm(vec_vel_intp, axis=1)) + / isce3.core.speed_of_light * self.radar_center_frequency + * self.azimuth_steer_rate) + + kappa_steer_grid = (kappa_steer_vec[..., np.newaxis] + * np.ones(grid_tau.shape[1])[np.newaxis, ...]) + + t_burst = (grid_t[0, 0] + grid_t[-1, 0]) / 2.0 + index_mid_burst_t = int(grid_t.shape[0] / 2 + 0.5) + tau_burst = (grid_tau[index_mid_burst_t, 0] + + grid_tau[index_mid_burst_t, -1]) / 2.0 + tau0_fdc_burst = tau0_fdc_interp[index_mid_burst_t] + tau0_fm_rate_burst = tau0_ka_interp[index_mid_burst_t] + + a0_burst, a1_burst, a2_burst = dc_coeffs[index_mid_burst_t, :] + b0_burst, b1_burst, b2_burst = fm_rate_coeffs[index_mid_burst_t, :] + + kappa_annotation_burst = (b0_burst + + b1_burst * (tau_burst - tau0_fm_rate_burst) + + b2_burst * (tau_burst - tau0_fm_rate_burst)**2) + + kappa_annotation_grid = _evaluate_polynomial_array(fm_rate_coeffs, + grid_tau, + tau0_ka_interp) + + grid_kappa_t = ( (kappa_annotation_grid * kappa_steer_grid) + / (kappa_annotation_grid - kappa_steer_grid)) + freq_dcg_burst = ( a0_burst + + a1_burst * (tau_burst - tau0_fdc_burst) + + a2_burst * (tau_burst - tau0_fdc_burst)**2) + + grid_freq_dcg = _evaluate_polynomial_array(dc_coeffs, + grid_tau, + tau0_fdc_interp) + + grid_freq_dc = (grid_freq_dcg + + grid_kappa_t * ((grid_t - t_burst) + + (freq_dcg_burst / kappa_annotation_burst + - grid_freq_dcg / kappa_annotation_grid))) + + lat_map = gdal.Open(list_filename_llh[0], gdal.GA_ReadOnly).ReadAsArray() + lon_map = gdal.Open(list_filename_llh[1], gdal.GA_ReadOnly).ReadAsArray() + hgt_map = gdal.Open(list_filename_llh[2], gdal.GA_ReadOnly).ReadAsArray() + + # Clean up the temporary directory in case it exists. + if temp_dir_obj is not None: + temp_dir_obj.cleanup() + + x_ecef, y_ecef, z_ecef = _llh_to_ecef(lat_map, + lon_map, + hgt_map, + ellipsoid) + + # Populate the position, velocity, and + # acceleration of the satellite to the correction grid + tau_ones = np.ones(grid_tau.shape[1])[np.newaxis, ...] + x_s, y_s, z_s = [vec_position_intp[:, i][..., np.newaxis] * tau_ones + for i in range(3)] + vx_s, vy_s, vz_s = [vec_vel_intp[:, i][..., np.newaxis] * tau_ones + for i in range(3)] + ax_s, ay_s, az_s = [vec_acceleration_intp[:, i][..., np.newaxis] * tau_ones + for i in range(3)] + + mag_xs_xg = np.sqrt( (x_s - x_ecef)**2 + + (y_s - y_ecef)**2 + + (z_s - z_ecef)**2) + + dotp_dxsg_acc = ( (x_s - x_ecef) * ax_s + + (y_s - y_ecef) * ay_s + + (z_s - z_ecef) * az_s) + + kappa_annotation_true = ( -(2 * self.radar_center_frequency) + / (isce3.core.speed_of_light * mag_xs_xg) + * (dotp_dxsg_acc + (vx_s**2 + vy_s**2 + vz_s**2))) + + delta_t_freq_mm = (grid_freq_dc + * (-1 / kappa_annotation_grid + + 1 / kappa_annotation_true)) + + # Prepare to export to LUT2d + vec_range = grid_tau[index_mid_burst_t, :] * isce3.core.speed_of_light / 2.0 + + return isce3.core.LUT2d(vec_range, vec_t, delta_t_freq_mm) + @property def sensing_mid(self): diff --git a/src/s1reader/s1_orbit.py b/src/s1reader/s1_orbit.py index 780a3e7..d1c804f 100644 --- a/src/s1reader/s1_orbit.py +++ b/src/s1reader/s1_orbit.py @@ -219,7 +219,6 @@ def get_file_name_tokens(zip_path: str) -> [str, list[datetime.datetime]]: item_valid = lambda item, sat_id: os.path.isfile(item) and sat_id in os.path.basename(item) - def get_orbit_file_from_dir(zip_path: str, orbit_dir: str, auto_download: bool = False) -> str: '''Get orbit state vector list for a given swath. @@ -274,7 +273,6 @@ def get_orbit_file_from_dir(zip_path: str, orbit_dir: str, auto_download: bool = return orbit_file - def get_orbit_file_from_list(zip_path: str, orbit_file_list: list) -> str: '''Get orbit file for a given S-1 swath from a list of files diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index b128dcf..7369a51 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -18,7 +18,8 @@ from s1reader import s1_annotation # to access __file__ from s1reader.s1_annotation import ProductAnnotation, NoiseAnnotation,\ CalibrationAnnotation, AuxCal, \ - BurstCalibration, BurstEAP, BurstNoise + BurstCalibration, BurstEAP, BurstNoise,\ + BurstExtendedCoeffs from s1reader.s1_burst_slc import Doppler, Sentinel1BurstSlc from s1reader.s1_burst_id import S1BurstId @@ -62,7 +63,8 @@ def parse_polynomial_element(elem, poly_name): half_c = 0.5 * isce3.core.speed_of_light r0 = half_c * float(elem.find('t0').text) - # NOTE Format of the azimuth FM rate polynomials has changed when IPF version was somewhere between 2.36 and 2.82 + # NOTE Format of the azimuth FM rate polynomials has changed + # when IPF version was somewhere between 2.36 and 2.82 if elem.find(poly_name) is None: # before the format change i.e. older IPF coeffs = [float(x.text) for x in elem[2:]] else: # after the format change i.e. newer IPF @@ -470,7 +472,8 @@ def get_track_burst_num(track_burst_num_file: str = esa_track_burst_id_file): return track_burst_num def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, - iw2_annotation_path: str, open_method=open, flag_apply_eap: bool = True): + iw2_annotation_path: str, open_method=open, + flag_apply_eap: bool = True): '''Parse bursts in Sentinel-1 annotation XML. Parameters: @@ -515,7 +518,8 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, # load the Calibraton annotation try: calibration_annotation_path =\ - annotation_path.replace('annotation/', 'annotation/calibration/calibration-') + 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 =\ @@ -526,7 +530,8 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, # load the Noise annotation try: - noise_annotation_path = annotation_path.replace('annotation/', 'annotation/calibration/noise-') + 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, @@ -643,10 +648,11 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, azimuth_time_interval) doppler = Doppler(poly1d, lut2d) + sensing_duration = datetime.timedelta( + seconds=n_lines * azimuth_time_interval) if len(osv_list) > 0: # get orbit from state vector list/element tree - sensing_duration = datetime.timedelta( - seconds=n_lines * azimuth_time_interval) + orbit = get_burst_orbit(sensing_start, sensing_start + sensing_duration, osv_list) else: @@ -671,12 +677,15 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, burst_calibration = None else: burst_calibration = BurstCalibration.from_calibration_annotation(calibration_annotation, - sensing_start) + sensing_start) if noise_annotation is None: burst_noise = None else: - burst_noise = BurstNoise.from_noise_annotation(noise_annotation, sensing_start, - i*n_lines, (i+1)*n_lines-1, ipf_version) + burst_noise = BurstNoise.from_noise_annotation(noise_annotation, + sensing_start, + i*n_lines, + (i+1)*n_lines-1, + ipf_version) if aux_cal_subswath is None: # Not applying EAP correction; (IPF high enough or user turned that off) # No need to fill in `burst_aux_cal` @@ -686,6 +695,13 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, aux_cal_subswath, sensing_start) + # Extended FM and DC coefficient information + extended_coeffs = BurstExtendedCoeffs.from_polynomial_lists( + az_fm_rate_list, + doppler_list, + sensing_start, + sensing_start + sensing_duration) + bursts[i] = Sentinel1BurstSlc(ipf_version, sensing_start, radar_freq, wavelength, azimuth_steer_rate, azimuth_time_interval, slant_range_time, starting_range, iw2_mid_range, @@ -698,7 +714,8 @@ 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, burst_noise, burst_aux_cal) + burst_calibration, burst_noise, burst_aux_cal, + extended_coeffs) return bursts @@ -745,6 +762,8 @@ def load_bursts(path: str, orbit_path: str, swath_num: int, pol: str = 'vv', returned. Default of None returns all bursts. If not all burst IDs are found, a list containing found bursts will be returned (empty list if none are found). + flag_apply_eap: bool + Turn on/off EAP related features (AUX_CAL loader) Returns: -------- @@ -772,9 +791,11 @@ def load_bursts(path: str, orbit_path: str, swath_num: int, pol: str = 'vv', if not os.path.exists(path): raise FileNotFoundError(f'{path} not found') elif os.path.isdir(path): - bursts = _burst_from_safe_dir(path, id_str, orbit_path, flag_apply_eap) + bursts = _burst_from_safe_dir(path, id_str, orbit_path, + flag_apply_eap) elif os.path.isfile(path): - bursts = _burst_from_zip(path, id_str, orbit_path, flag_apply_eap) + bursts = _burst_from_zip(path, id_str, orbit_path, + flag_apply_eap) else: raise ValueError(f'{path} is unsupported') @@ -799,7 +820,8 @@ def load_bursts(path: str, orbit_path: str, swath_num: int, pol: str = 'vv', return bursts -def _burst_from_zip(zip_path: str, id_str: str, orbit_path: str, flag_apply_eap: bool): +def _burst_from_zip(zip_path: str, id_str: str, orbit_path: str, + flag_apply_eap: bool): '''Find bursts in a Sentinel-1 zip file. Parameters: @@ -810,6 +832,8 @@ def _burst_from_zip(zip_path: str, id_str: str, orbit_path: str, flag_apply_eap: Identification of desired burst. Format: iw[swath_num]-slc-[pol] orbit_path : str Path the orbit file. + flag_apply_eap: bool + Turn on/off EAP related features (AUX_CAL loader) Returns: -------- @@ -841,8 +865,8 @@ def _burst_from_zip(zip_path: str, id_str: str, orbit_path: str, flag_apply_eap: flag_apply_eap=flag_apply_eap) return bursts - -def _burst_from_safe_dir(safe_dir_path: str, id_str: str, orbit_path: str, flag_apply_eap: bool): +def _burst_from_safe_dir(safe_dir_path: str, id_str: str, orbit_path: str, + flag_apply_eap: bool): '''Find bursts in a Sentinel-1 SAFE structured directory. Parameters: @@ -853,6 +877,8 @@ def _burst_from_safe_dir(safe_dir_path: str, id_str: str, orbit_path: str, flag_ Identification of desired burst. Format: iw[swath_num]-slc-[pol] orbit_path : str Path the orbit file. + flag_apply_eap: bool + Turn on/off EAP related features (AUX_CAL loader) Returns: --------