diff --git a/src/s1reader/s1_orbit.py b/src/s1reader/s1_orbit.py index 2643e30..7bb0d1b 100644 --- a/src/s1reader/s1_orbit.py +++ b/src/s1reader/s1_orbit.py @@ -1,5 +1,8 @@ import datetime import os +import xml.etree.ElementTree as ET + +import isce3 # date format used in file names FMT = "%Y%m%dT%H%M%S" @@ -121,3 +124,42 @@ def get_orbit_file_from_dir(path: str, orbit_dir: str) -> str: path, [f'{orbit_dir}/{item}' for item in os.listdir(orbit_dir)]) return orbit_path + +def burst_orbit_from_file(sensing_start, sensing_stop, osv_list: ET.Element): + '''Init and return ISCE3 orbit from element in orbit XML file. + + Parameters: + ----------- + sensing_start : datetime.datetime + Sensing start of burst; taken from azimuth time + sensing_stop : datetime.datetime + Sensing stop of burst + osv_list : xml.etree.ElementTree.Element + ElementTree containing orbit state vectors + + Returns: + -------- + _ : datetime + Sensing mid as datetime object. + ''' + fmt = "UTC=%Y-%m-%dT%H:%M:%S.%f" + orbit_sv = [] + # add start & end padding to ensure sufficient number of orbit points + pad = datetime.timedelta(seconds=60) + for osv in osv_list: + t_orbit = datetime.datetime.strptime(osv[1].text, fmt) + + if t_orbit > sensing_stop + pad: + break + + if t_orbit > sensing_start - pad: + pos = [float(osv[i].text) for i in range(4,7)] + vel = [float(osv[i].text) for i in range(7,10)] + orbit_sv.append(isce3.core.StateVector(isce3.core.DateTime(t_orbit), + pos, vel)) + + # use list of stateVectors to init and return isce3.core.Orbit + time_delta = datetime.timedelta(days=2) + ref_epoch = isce3.core.DateTime(sensing_start - time_delta) + + return isce3.core.Orbit(orbit_sv, ref_epoch) diff --git a/src/s1reader/s1_reader.py b/src/s1reader/s1_reader.py index f469ac4..792c007 100644 --- a/src/s1reader/s1_reader.py +++ b/src/s1reader/s1_reader.py @@ -10,9 +10,10 @@ import isce3 from nisar.workflows.stage_dem import check_dateline from s1reader.s1_burst_slc import Doppler, Sentinel1BurstSlc +from s1reader.s1_orbit import burst_orbit_from_file # 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"): +def as_datetime(t_str, fmt="%Y-%m-%dT%H:%M:%S.%f"): '''Parse given time string to datetime.datetime object. Parameters: @@ -126,7 +127,7 @@ def doppler_poly1d_to_lut2d(doppler_poly1d, starting_slant_range, return isce3.core.LUT2d(slant_ranges, az_times, np.vstack((freq_1d, freq_1d))) -def get_burst_orbit(sensing_start, sensing_stop, osv_list: ET.Element): +def parse_annotation_orbit(sensing_start, sensing_stop, osv_list: ET.Element): '''Init and return ISCE3 orbit. Parameters: @@ -143,21 +144,18 @@ def get_burst_orbit(sensing_start, sensing_stop, osv_list: ET.Element): _ : datetime Sensing mid as datetime object. ''' - fmt = "UTC=%Y-%m-%dT%H:%M:%S.%f" - orbit_sv = [] + fmt = "%Y-%m-%dT%H:%M:%S.%f" + n_svs = int(osv_list.attrib['count']) + orbit_sv = [[]] * n_svs # add start & end padding to ensure sufficient number of orbit points pad = datetime.timedelta(seconds=60) - for osv in osv_list: - t_orbit = datetime.datetime.strptime(osv[1].text, fmt) + for i, osv in enumerate(osv_list): + t_orbit = datetime.datetime.strptime(osv[0].text, fmt) - if t_orbit > sensing_stop + pad: - break - - if t_orbit > sensing_start - pad: - pos = [float(osv[i].text) for i in range(4,7)] - vel = [float(osv[i].text) for i in range(7,10)] - orbit_sv.append(isce3.core.StateVector(isce3.core.DateTime(t_orbit), - pos, vel)) + pos = [float(osv[2][i].text) for i in range(3)] + vel = [float(osv[3][i].text) for i in range(3)] + orbit_sv[i] = isce3.core.StateVector(isce3.core.DateTime(t_orbit), + pos, vel) # use list of stateVectors to init and return isce3.core.Orbit time_delta = datetime.timedelta(days=2) @@ -336,8 +334,13 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, iw2_mid_range = iw2_starting_range + 0.5 * iw2_n_samples * range_pxl_spacing # find orbit state vectors in 'Data_Block/List_of_OSVs' - orbit_tree = ET.parse(orbit_path) - osv_list = orbit_tree.find('Data_Block/List_of_OSVs') + if orbit_path is None: + osv_list = tree.find('generalAnnotation/orbitList') + get_orbit = parse_annotation_orbit + else: + orbit_tree = ET.parse(orbit_path) + osv_list = orbit_tree.find('Data_Block/List_of_OSVs') + get_orbit = burst_orbit_from_file # load individual burst burst_list_elements = tree.find('swathTiming/burstList') @@ -369,8 +372,8 @@ def burst_from_xml(annotation_path: str, orbit_path: str, tiff_path: str, # get orbit from state vector list/element tree sensing_duration = datetime.timedelta( seconds=n_samples * azimuth_time_interval) - orbit = get_burst_orbit(sensing_start, sensing_start + sensing_duration, - osv_list) + orbit = get_orbit(sensing_start, sensing_start + sensing_duration, + osv_list) # determine burst offset and dimensions # TODO move to own method @@ -425,7 +428,7 @@ def _is_zip_annotation_xml(path: str, id_str: str) -> bool: return True return False -def load_bursts(path: str, orbit_path: str, swath_num: int, pol: str='vv', +def load_bursts(path: str, swath_num: int, pol: str='vv', orbit_path: str=None, burst_ids: list[str]=None): '''Find bursts in a Sentinel-1 zip file or a SAFE structured directory. @@ -433,12 +436,12 @@ def load_bursts(path: str, orbit_path: str, swath_num: int, pol: str='vv', ----------- path : str Path to Sentinel-1 zip file or SAFE directory - orbit_path : str - Path the orbit file. swath_num : int Integer of subswath of desired burst. {1, 2, 3} pol : str Polarization of desired burst. {hh, vv, hv, vh} + orbit_path : str + Path the orbit file. If None, annotation XML orbit is parsed. burst_ids : list[str] List of burst IDs for which their Sentinel1BurstSlc objects will be returned. Default of None returns all bursts. Empty list returned if @@ -466,6 +469,11 @@ def load_bursts(path: str, orbit_path: str, swath_num: int, pol: str='vv', if pol not in pols: raise ValueError(f"polarization not in {pols}") + if orbit_path is None: + warn_str = 'No external orbit given. ' + warn_str += 'Orbit will be constructed from annotation XML.' + warnings.warn(warn_str) + id_str = f'iw{swath_num}-slc-{pol}' if not os.path.exists(path): @@ -506,7 +514,7 @@ def _burst_from_zip(zip_path: str, id_str: str, orbit_path: str): id_str: str Identifcation of desired burst. Format: iw[swath_num]-slc-[pol] orbit_path : str - Path the orbit file. + Path the orbit file. If None, annotation XML orbit is parsed. Returns: -------- @@ -547,7 +555,7 @@ def _burst_from_safe_dir(safe_dir_path: str, id_str: str, orbit_path: str): id_str: str Identifcation of desired burst. Format: iw[swath_num]-slc-[pol] orbit_path : str - Path the orbit file. + Path the orbit file. If None, annotation XML orbit is parsed. Returns: --------