Index

estimate_ionospheric_delay(ifg_file_list, slc_files, tec_files, geom_files, reference_point, output_dir, epsg, bounds, file_date_fmt='%Y%m%d')

Estimate the range delay (in meters) caused by ionosphere for each interferogram.

Note: Currently this workflow only supports ionospheric corrections on OPERA CSLC input datasets.

Parameters:

Name	Type	Description	Default
ifg_file_list	list[Path]	List of interferogram files.	required
slc_files	Dict[date, list[Filename]]	Dictionary of SLC files indexed by date.	required
tec_files	Dict[date, list[Filename]]	Dictionary of TEC files indexed by date.	required
geom_files	list[Path]	Dictionary of geometry files with height and incidence angle, or los_east and los_north.	required
reference_point	tuple[int, int]	Reference point (row, col) used during time series inversion. If not provided, no spatial referencing is performed.	required
output_dir	Path	Output directory.	required
epsg	int	the EPSG code of the input data	required
bounds	Bbox	Output bounds.	required
file_date_fmt	str	The format string to use when parsing the dates from the file names. Default is "%Y%m%d".	'%Y%m%d'

Returns:

Type	Description
list[Path]	List of newly created ionospheric phase delay corrections.

Source code in src/dolphin/atmosphere/ionosphere.py

def estimate_ionospheric_delay(
    ifg_file_list: Sequence[Path],
    slc_files: Mapping[tuple[datetime.datetime], Sequence[Filename]],
    tec_files: Mapping[tuple[datetime.datetime], Sequence[Filename]],
    geom_files: dict[str, Path],
    reference_point: ReferencePoint | None,
    output_dir: Path,
    epsg: int,
    bounds: Bbox,
    file_date_fmt: str = "%Y%m%d",
) -> list[Path]:
    """Estimate the range delay (in meters) caused by ionosphere for each interferogram.

    Note: Currently this workflow only supports ionospheric corrections
    on OPERA CSLC input datasets.

    Parameters
    ----------
    ifg_file_list : list[Path]
        List of interferogram files.
    slc_files : Dict[datetime.date, list[Filename]]
        Dictionary of SLC files indexed by date.
    tec_files : Dict[datetime.date, list[Filename]]
        Dictionary of TEC files indexed by date.
    geom_files : list[Path]
        Dictionary of geometry files with height and incidence angle, or
        los_east and los_north.
    reference_point : tuple[int, int], optional
        Reference point (row, col) used during time series inversion.
        If not provided, no spatial referencing is performed.
    output_dir : Path
        Output directory.
    epsg : int
        the EPSG code of the input data
    bounds : Bbox
        Output bounds.
    file_date_fmt : str, optional
        The format string to use when parsing the dates from the file names.
        Default is "%Y%m%d".

    Returns
    -------
    list[Path]
        List of newly created ionospheric phase delay corrections.

    """
    if epsg != 4326:
        left, bottom, right, top = transform_bounds(
            CRS.from_epsg(epsg), CRS.from_epsg(4326), *bounds
        )
    else:
        left, bottom, right, top = bounds

    # Read the incidence angle
    if "los_east" in geom_files:
        # ISCE3 geocoded products
        los_east = io.load_gdal(geom_files["los_east"]).astype(np.float32)
        los_north = io.load_gdal(geom_files["los_north"]).astype(np.float32)
        inc_angle = np.arccos(np.sqrt(1 - los_east**2 - los_north**2)) * 180 / np.pi
    else:
        # ISCE2 radar coordinate
        inc_angle = io.load_gdal(geom_files["incidence_angle"]).astype(np.float32)

    iono_inc_angle = incidence_angle_ground_to_iono(inc_angle)

    # frequency
    for key in slc_files:
        if "compressed" not in str(slc_files[key][0]).lower():
            one_of_slcs = slc_files[key][0]
            break

    wavelength = oput.get_radar_wavelength(one_of_slcs)
    freq = SPEED_OF_LIGHT / wavelength

    output_iono = output_dir / "ionosphere"
    output_iono.mkdir(exist_ok=True)

    output_paths: list[Path] = []

    num_lats, num_lons = iono_inc_angle.shape

    downsample_factor: int = 10
    # Create downsampled grid for efficiency
    num_lats, num_lons = iono_inc_angle.shape
    ds_num_lats = max(5, num_lats // downsample_factor)
    ds_num_lons = max(5, num_lons // downsample_factor)

    # Create the downsampled grid
    ds_out_lats = np.linspace(top, bottom, ds_num_lats)
    ds_out_lons = np.linspace(left, right, ds_num_lons)
    ds_lat_grid, ds_lon_grid = np.meshgrid(ds_out_lats, ds_out_lons, indexing="ij")

    # Downsample the incidence angle
    ds_inc_angle = zoom(
        iono_inc_angle, (ds_num_lats / num_lats, ds_num_lons / num_lons), order=1
    )
    # Ensure same size as the lat/lon grids
    ds_inc_angle = ds_inc_angle[:ds_num_lats, :ds_num_lons]

    for ifg in ifg_file_list:
        ref_date, sec_date = get_dates(ifg, fmt=file_date_fmt)

        name = f"{format_date_pair(ref_date, sec_date)}_ionoDelay.tif"
        iono_delay_product_name = output_iono / name

        output_paths.append(iono_delay_product_name)
        if iono_delay_product_name.exists():
            logger.info(
                "Tropospheric correction for interferogram "
                f"{format_date_pair(ref_date, sec_date)} already exists, skipping"
            )
            continue

        # The keys in slc_files do not necessarily have one date,
        # it means with the production file naming convention,
        # there will be multiple dates stored in the keys from get_dates
        # function. So in the following if we set reference_date = (ref_date, ),
        # there will be an error that it does not find the key.
        # Examples of the file naming convention for CSLC and compressed cslc is:
        # OPERA_L2_CSLC-S1_T042-088905-IW1\
        #                       _20221119T000000Z_20221120T000000Z_S1A_VV_v1.0.h5
        # OPERA_L2_COMPRESSED-CSLC-S1_T042-088905-IW1_20221107T000000Z_\
        #           20221107T000000Z_20230506T000000Z_20230507T000000Z_VV_v1.0.h5
        reference_date = next(key for key in slc_files if ref_date in key)
        # temporary fix for compressed SLCs while the required metadata i not included
        # in them. there will be modification in a future PR to add metadata to
        # compressed SLCs
        secondary_date = next(
            key
            for key in slc_files
            if sec_date in key and "compressed" not in str(slc_files[key][0]).lower()
        )

        secondary_time = oput.get_zero_doppler_time(slc_files[secondary_date][0])
        if "compressed" in str(slc_files[reference_date][0]).lower():
            # this is for when we have compressed slcs but the actual
            # reference date does not exist in the input data
            reference_time = secondary_time
        else:
            reference_time = oput.get_zero_doppler_time(slc_files[reference_date][0])

        # Calculate interpolated range delays
        vtec_reference = read_zenith_tec(
            reference_time, tec_files[(ref_date,)][0], ds_lat_grid, ds_lon_grid
        )
        # Make the downsampled (ds_) version
        ds_range_delay_reference = vtec_to_range_delay(
            vtec_reference, ds_inc_angle, freq
        )
        vtec_secondary = read_zenith_tec(
            secondary_time, tec_files[(sec_date,)][0], ds_lat_grid, ds_lon_grid
        )
        ds_range_delay_secondary = vtec_to_range_delay(
            vtec_secondary, ds_inc_angle, freq
        )

        # For displacement output, positive corresponds to motion toward the satellite
        # If range_delay_secondary is smaller than range_delay_reference, then the
        # resulting delay difference is positive. A smaller excess secondary delay
        # looks the same as motion toward the satellite.
        ds_ifg_iono_range_delay = ds_range_delay_reference - ds_range_delay_secondary

        # Finally upsample the end result:
        ifg_iono_range_delay = zoom(
            ds_ifg_iono_range_delay,
            (num_lats / ds_num_lats, num_lons / ds_num_lons),
            order=1,
        )

        if reference_point is not None:
            ref_row, ref_col = reference_point
            ifg_iono_range_delay -= ifg_iono_range_delay[ref_row, ref_col]

        # Write 2D tropospheric correction layer to disc
        io.write_arr(
            arr=ifg_iono_range_delay,
            output_name=iono_delay_product_name,
            like_filename=ifg,
            units="meters",
        )

    return output_paths