From 2e1a46708f8e0693cb70482a69bcb95387e36253 Mon Sep 17 00:00:00 2001
From: Jungkyo Jung <jungkyoj@nisar-adt-dev-5.jpl.nasa.gov>
Date: Sat, 11 Apr 2026 19:33:57 +0000
Subject: [PATCH 1/6] fix a bug in projecing water to cryosphere region

---
 python/packages/isce3/unwrap/preprocess.py    | 264 ++++++++++++++----
 python/packages/nisar/workflows/insar.py      |  10 +-
 python/packages/nisar/workflows/ionosphere.py |  79 ++++--
 python/packages/nisar/workflows/unwrap.py     |   1 +
 4 files changed, 272 insertions(+), 82 deletions(-)

diff --git a/python/packages/isce3/unwrap/preprocess.py b/python/packages/isce3/unwrap/preprocess.py
index 618dbe3e5..a256ae900 100644
--- a/python/packages/isce3/unwrap/preprocess.py
+++ b/python/packages/isce3/unwrap/preprocess.py
@@ -1,7 +1,7 @@
 import pathlib
 import journal
 import numpy as np
-from osgeo import gdal
+from osgeo import gdal, osr
 
 from scipy.ndimage import median_filter, map_coordinates
 
@@ -222,66 +222,227 @@ def _get_gdal_raster_shape_type(raster_path):
     return data_shape, np_data_type
 
 
-def _read_gdal_with_bbox(gdal_raster, bbox):
-    '''
-    Extract image from the gdal-supported file with bbox
+def _read_gdal_with_bbox(input_raster, bbox, bbox_epsg=4326):
+    """
+    Read only the raster subset intersecting the input bbox, and return it
+    reprojected to bbox_epsg.
 
     Parameters
     ----------
-    gdal_raster: osgeo.gdal.Dataset
-        gdal dataset to extract the subset image
-    bbox: list
-        list of [xmin, ymin, xmax, ymax]
+    input_raster : gdal.Dataset
+        Input GDAL raster
+    bbox : list[float]
+        [xmin, ymin, xmax, ymax]
+    bbox_epsg : int
+        EPSG of bbox coordinates
 
     Returns
     -------
-    subset_data: numpy.ndarray
-        Median absolute deviation of `arr`
-    [sub_x0, sub_y0, dx, dy]: list
-        sub_x0: x coordinate of upper left
-        sub_y0: y coordinate of upper left
-        dx: x spacing
-        dy: y spacing
-    '''
-    xmin, ymin, xmax, ymax = bbox
+    arr : numpy.ndarray
+        Raster subset reprojected to bbox_epsg
+    raster_info : list[float]
+        [block_x0, block_y0, block_dx, block_dy] in bbox_epsg coordinates
+    """
+    gt = input_raster.GetGeoTransform()
+    proj = input_raster.GetProjection()
+    band = input_raster.GetRasterBand(1)
+
+    if band is None:
+        raise RuntimeError("Failed to access raster band.")
+
+    if gt is None:
+        raise RuntimeError("Raster geotransform is missing.")
+
+    # north-up only, same practical assumption as most GeoTIFF use cases here
+    if gt[2] != 0 or gt[4] != 0:
+        raise NotImplementedError(
+            "_read_gdal_with_bbox currently supports only north-up rasters."
+        )
+
+    raster_srs = osr.SpatialReference()
+    raster_srs.ImportFromWkt(proj)
+
+    try:
+        raster_srs.AutoIdentifyEPSG()
+    except Exception:
+        pass
+
+    raster_epsg = raster_srs.GetAuthorityCode(None)
+    if raster_epsg is None:
+        raise RuntimeError("Could not determine raster EPSG from projection.")
+    raster_epsg = int(raster_epsg)
+
+    xmin, ymin, xmax, ymax = map(float, bbox)
+
+    # target SRS = bbox_epsg
+    dst_srs = osr.SpatialReference()
+    dst_srs.ImportFromEPSG(int(bbox_epsg))
+    try:
+        dst_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+        raster_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    except Exception:
+        pass
+
+    # First, transform bbox corners to raster CRS only to check overlap
+    if bbox_epsg != raster_epsg:
+        tx_to_raster = osr.CoordinateTransformation(dst_srs, raster_srs)
+
+        corners = [
+            tx_to_raster.TransformPoint(xmin, ymin)[:2],
+            tx_to_raster.TransformPoint(xmin, ymax)[:2],
+            tx_to_raster.TransformPoint(xmax, ymin)[:2],
+            tx_to_raster.TransformPoint(xmax, ymax)[:2],
+        ]
+        xs = [p[0] for p in corners]
+        ys = [p[1] for p in corners]
+        xmin_src, xmax_src = min(xs), max(xs)
+        ymin_src, ymax_src = min(ys), max(ys)
+    else:
+        xmin_src, xmax_src = xmin, xmax
+        ymin_src, ymax_src = ymin, ymax
+
+    x0 = gt[0]
+    dx = gt[1]
+    y0 = gt[3]
+    dy = gt[5]
+
+    raster_width = input_raster.RasterXSize
+    raster_height = input_raster.RasterYSize
+
+    raster_xmin = min(x0, x0 + raster_width * dx)
+    raster_xmax = max(x0, x0 + raster_width * dx)
+    raster_ymin = min(y0, y0 + raster_height * dy)
+    raster_ymax = max(y0, y0 + raster_height * dy)
+
+    xmin_i = max(xmin_src, raster_xmin)
+    xmax_i = min(xmax_src, raster_xmax)
+    ymin_i = max(ymin_src, raster_ymin)
+    ymax_i = min(ymax_src, raster_ymax)
+
+    if xmin_i >= xmax_i or ymin_i >= ymax_i:
+        raise ValueError("Input bbox does not overlap raster.")
+
+    src_nodata = band.GetNoDataValue()
+    if src_nodata is None:
+        src_nodata = 0
+
+    # Use source pixel size magnitude to define output resolution in target CRS.
+    # This keeps behavior simple and avoids very strange defaults from Warp.
+    # For your case (bbox_epsg=4326), output grid becomes lon/lat.
+    if bbox_epsg == raster_epsg:
+        out_xres = abs(dx)
+        out_yres = abs(dy)
+    else:
+        tx_to_bbox = osr.CoordinateTransformation(raster_srs, dst_srs)
+
+        # transform two neighboring source points to estimate target resolution
+        p00 = tx_to_bbox.TransformPoint(x0, y0)[:2]
+        p10 = tx_to_bbox.TransformPoint(x0 + dx, y0)[:2]
+        p01 = tx_to_bbox.TransformPoint(x0, y0 + dy)[:2]
+
+        est_dx = abs(p10[0] - p00[0])
+        est_dy = abs(p01[1] - p00[1])
+
+        # fallback in case estimate becomes zero or numerically unstable
+        if not np.isfinite(est_dx) or est_dx <= 0:
+            est_dx = max((xmax - xmin) / 1000.0, 1e-6)
+        if not np.isfinite(est_dy) or est_dy <= 0:
+            est_dy = max((ymax - ymin) / 1000.0, 1e-6)
+
+        out_xres = est_dx
+        out_yres = est_dy
+
+    # Warp directly to the requested bbox / requested CRS
+    warped_ds = gdal.Warp(
+        "",
+        input_raster,
+        format="MEM",
+        dstSRS=dst_srs.ExportToWkt(),
+        outputBounds=(xmin, ymin, xmax, ymax),
+        outputBoundsSRS=dst_srs.ExportToWkt(),
+        xRes=out_xres,
+        yRes=out_yres,
+        resampleAlg=gdal.GRA_NearestNeighbour,
+        srcNodata=src_nodata,
+        dstNodata=src_nodata,
+        targetAlignedPixels=False,
+        multithread=False,
+    )
+
+    if warped_ds is None:
+        raise RuntimeError("gdal.Warp failed to create warped subset.")
+
+    warped_band = warped_ds.GetRasterBand(1)
+    if warped_band is None:
+        raise RuntimeError("Failed to access warped raster band.")
+
+    arr = warped_band.ReadAsArray()
+    if arr is None:
+        raise RuntimeError("Failed to read warped raster window.")
+
+    warped_gt = warped_ds.GetGeoTransform()
+    if warped_gt is None:
+        raise RuntimeError("Warped raster geotransform is missing.")
+
+    block_x0 = warped_gt[0]
+    block_y0 = warped_gt[3]
+    block_dx = warped_gt[1]
+    block_dy = warped_gt[5]
+
+    return arr, [block_x0, block_y0, block_dx, block_dy]
+
+
+def _find_rdr2geo_paths(scratch_path, freq):
+    """
+    Find x.rdr and y.rdr files for the given frequency inside scratch_path.
 
-    geotransform = gdal_raster.GetGeoTransform()
-    x0 = geotransform[0]
-    y0 = geotransform[3]
-    dx = geotransform[1]
-    dy = geotransform[5]
+    The function searches recursively because the files may exist in several
+    possible directories such as:
 
-    idx_start = int(np.floor((xmin - x0) / dx))
-    idx_end = int(np.ceil((xmax - x0) / dx))
+        scratch/rdr2geo/freqA/
+        scratch/ionosphere/main_diff_ms_band/rdr2geo/freqB/
+        scratch/ionosphere/main_side_band/rdr2geo/freqB/
 
-    if dy > 0:
-        idy_start = int(np.floor((ymin - y0) / dy))
-        idy_end = int(np.ceil((ymax - y0) / dy))
-        sub_y0 = idy_start*dy + y0
-    else:
-        idy_start = int(np.floor((ymax - y0) / dy))
-        idy_end = int(np.ceil((ymin - y0) / dy))
+    Parameters
+    ----------
+    scratch_path : pathlib.Path
+    freq : str
+
+    Returns
+    -------
+    dict
+        {"x": path_to_x_rdr, "y": path_to_y_rdr}
+
+    Raises
+    ------
+    FileNotFoundError
+        If the rdr2geo files cannot be located.
+    """
+
+    scratch_path = pathlib.Path(scratch_path)
+
+    candidates = list(
+        scratch_path.glob(f"**/rdr2geo/freq{freq}/x.rdr")
+    )
 
-    idx_start = max(0, idx_start)
-    idy_start = max(0, idy_start)
+    if not candidates:
+        raise FileNotFoundError(
+            f"Could not find any x.rdr under {scratch_path} "
+            f"for frequency {freq}."
+        )
 
-    x_width = idx_end - idx_start
-    y_length = idy_end - idy_start
+    # choose the shallowest path (closest to scratch root)
+    candidates.sort(key=lambda p: len(p.parts))
+    x_path = candidates[0]
 
-    if x_width > gdal_raster.RasterXSize:
-        x_width = gdal_raster.RasterXSize
-    if y_length > gdal_raster.RasterYSize:
-        y_length = gdal_raster.RasterYSize
+    y_path = x_path.parent / "y.rdr"
 
-    sub_x0 = idx_start*dx + x0
-    sub_y0 = idy_start*dy + y0
-    raster_band = gdal_raster.GetRasterBand(1)
-    subset_data = raster_band.ReadAsArray(idx_start,
-                                          idy_start,
-                                          x_width,
-                                          y_length)
+    if not y_path.exists():
+        raise FileNotFoundError(
+            f"Found {x_path} but corresponding y.rdr does not exist."
+        )
 
-    return subset_data, [sub_x0, sub_y0, dx, dy]
+    return {"x": str(x_path), "y": str(y_path)}
 
 
 def project_map_to_radar(cfg, input_data_path, freq):
@@ -303,7 +464,6 @@ def project_map_to_radar(cfg, input_data_path, freq):
         projected data into radar grid  absolute
     '''
     scratch_path = pathlib.Path(cfg['product_path_group']['scratch_path'])
-    rdr2geo_path = f'{scratch_path}/rdr2geo'
 
     az_looks = cfg["processing"]["crossmul"]["azimuth_looks"]
     rg_looks = cfg["processing"]["crossmul"]["range_looks"]
@@ -314,10 +474,8 @@ def project_map_to_radar(cfg, input_data_path, freq):
     if unw_rg_looks != 1:
         rg_looks = unw_rg_looks
 
-    # prepare input paths
-    topo_paths = {xy: f'{rdr2geo_path}/freq{freq}/{xy}.rdr' for xy in 'xy'}
+    topo_paths = _find_rdr2geo_paths(scratch_path, freq)
 
-    # get input shape and type - input type also output type
     _, output_dtype = _get_gdal_raster_shape_type(input_data_path)
     geo_data_raster = gdal.Open(input_data_path)
 
@@ -374,13 +532,12 @@ def project_map_to_radar(cfg, input_data_path, freq):
                     cval=np.nan,
                     prefilter=False)
 
-    # stack to make whole then return
     return output_arrays
 
 
 def interpret_subswath_mask(subswath_mask, nodata=255):
     """
-    Interprets a subswath mask integer by decoding its digits into boolean 
+    Interprets a subswath mask integer by decoding its digits into boolean
     flags indicating reference validity, secondary validity, and water
     presence.
 
@@ -421,3 +578,4 @@ def interpret_subswath_mask(subswath_mask, nodata=255):
     water = np.where(nd, False, water)
 
     return reference_valid, secondary_valid, water
+    
\ No newline at end of file
diff --git a/python/packages/nisar/workflows/insar.py b/python/packages/nisar/workflows/insar.py
index aba22f796..807383848 100644
--- a/python/packages/nisar/workflows/insar.py
+++ b/python/packages/nisar/workflows/insar.py
@@ -46,12 +46,6 @@ def run(cfg: dict, out_paths: dict, run_steps: dict):
     if run_steps['geo2rdr']:
         geo2rdr.run(cfg)
 
-    # Remove the rdr2geo scratch folder
-    rdr2geo_scratch_path = pathlib.Path(f"{scratch_path}/rdr2geo")
-    _remove_intermediate_dir(rdr2geo_scratch_path,
-                             intermediate_files_removal_flag,
-                             info_channel)
-
     if run_steps['prepare_insar_hdf5']:
         prepare_insar_hdf5.run(cfg)
 
@@ -140,8 +134,8 @@ def run(cfg: dict, out_paths: dict, run_steps: dict):
         # Geocode ROFF
         geocode_insar.run(cfg, out_paths['ROFF'], out_paths['GOFF'], InputProduct.ROFF)
 
-    # Remove the 'ionosphere' and 'geocode_corrections' scratch folders
-    for workflow_name in ['ionosphere','geocode_corrections']:
+    # Remove the 'ionosphere', 'rdr2geo' and 'geocode_corrections' scratch folders
+    for workflow_name in ['ionosphere', 'geocode_corrections', 'rdr2geo']:
         workflow_scratch_path = pathlib.Path(f"{scratch_path}/{workflow_name}")
         _remove_intermediate_dir(workflow_scratch_path,
                                  intermediate_files_removal_flag,
diff --git a/python/packages/nisar/workflows/ionosphere.py b/python/packages/nisar/workflows/ionosphere.py
index 539de3e6e..a74512ec6 100644
--- a/python/packages/nisar/workflows/ionosphere.py
+++ b/python/packages/nisar/workflows/ionosphere.py
@@ -537,8 +537,10 @@ def _to_complex_if_needed(arr):
 
                     # Ensure complex (convert real-valued phase in radians to
                     # complex phase)
-                    first_data_block = _to_complex_if_needed(first_data_block)
-                    second_data_block = _to_complex_if_needed(second_data_block)
+                    use_complex = np.iscomplexobj(first_data_block) or np.iscomplexobj(second_data_block)
+                    if use_complex:
+                        first_data_block = _to_complex_if_needed(first_data_block)
+                        second_data_block = _to_complex_if_needed(second_data_block)
 
                     # Resample first dataset (frequency A / low subband) to
                     # match the grid of the second dataset (frequency B / high subband).
@@ -580,7 +582,10 @@ def _to_complex_if_needed(arr):
                         invalid = None
 
                     # Compute the differential phase
-                    diff_phase = first_data_block * np.conj(second_data_block)
+                    if use_complex:
+                        diff_phase = first_data_block * np.conj(second_data_block)
+                    else:
+                        diff_phase = first_data_block - second_data_block
                     if invalid is not None:
                         diff_phase[invalid] = invalid_fill_value
 
@@ -883,13 +888,16 @@ def insar_ionosphere_pair(original_cfg, runw_hdf5):
                                             'ionosphere',
                                             'main_diff_ms_band',
                                             'RIFG.h5')
+                unwrap_needed = True
+                diff_phase_output = pathlib.Path(diff_dir, 'RIFG.h5')
             else:
                 out_paths = original_out_paths
                 new_scratch = orig_scratch_path
-                phase_second = out_paths['RIFG']
+                phase_second = out_paths['RUNW']
                 additional_runw = out_paths['RUNW']
 
-            diff_phase_output = pathlib.Path(diff_dir, 'RIFG.h5')
+                unwrap_needed = False
+                diff_phase_output = pathlib.Path(diff_dir, 'RUNW.h5')
             iono_insar_cfg['product_path_group'][
                 'scratch_path'] = diff_dir
             iono_insar_cfg['product_path_group'][
@@ -925,12 +933,18 @@ def insar_ionosphere_pair(original_cfg, runw_hdf5):
 
             second_data_path = []
             for pol_b in pol_list_b:
+                if rerun_insar_pairs > 0:
+                    dest_freq_path = f"{swath_path}/frequencyB"
+                    dest_pol_path = f"{dest_freq_path}/interferogram/{pol_b}"
+                    rifg_path_freq = f"{dest_pol_path}/wrappedInterferogram"
 
-                dest_freq_path = f"{swath_path}/frequencyB"
-                dest_pol_path = f"{dest_freq_path}/interferogram/{pol_b}"
-                rifg_path_freq = f"{dest_pol_path}/wrappedInterferogram"
+                    second_data_path.append(rifg_path_freq)
+                else:
+                    dest_freq_path = f"{runw_swath_path}/frequencyB"
+                    dest_pol_path = f"{dest_freq_path}/interferogram/{pol_b}"
+                    runw_path_b_freq = f"{dest_pol_path}/unwrappedPhase"
 
-                second_data_path.append(rifg_path_freq)
+                    second_data_path.append(runw_path_b_freq)
             second_slant_path = f"{dest_freq_path}/interferogram/slantRange"
             second_mask_path = f"{dest_freq_path}/interferogram/mask"
 
@@ -951,8 +965,8 @@ def insar_ionosphere_pair(original_cfg, runw_hdf5):
             # Since main_diff_low_high_subband method does not need to
             # unwrap low and high subband interferogram, but need to
             # unwrap the difference between low and high subband interferogram
-            unwrap.run(iono_insar_cfg, out_paths['RIFG'], out_paths['RUNW'])
-
+            if unwrap_needed:
+                unwrap.run(iono_insar_cfg, out_paths['RIFG'], out_paths['RUNW'])
     # restore original paths
     original_cfg['input_file_group']['reference_rslc_file'] = \
         partial_orig_cfg_dict['reference_rslc_file']
@@ -1022,7 +1036,8 @@ def run_insar_workflow(iono_insar_cfg, original_dict, out_paths,
     # decimate offsets for frequency B and create ionosphere layers
     if 'B' in iono_freq_pol:
         decimate_freq_a_offset(iono_insar_cfg, original_dict)
-
+    print(iono_insar_cfg['processing']['input_subset'][
+                        'list_of_frequencies'])
     if iono_insar_cfg['processing']['fine_resample']['enabled']:
         resample_slc_v2.run(iono_insar_cfg, 'fine')
     else:
@@ -1205,6 +1220,20 @@ def run(cfg: dict, runw_hdf5: str):
         filling_method=filling_method,
         outputdir=os.path.join(iono_path, iono_method))
 
+    # compute the full water masks
+    water_mask_b_blk = None
+    water_mask_a_blk = None
+    if "water" in mask_type and filter_bool:
+        water_mask_path = cfg[
+            "dynamic_ancillary_file_group"]["water_mask_file"]
+
+        water_distance_a = project_map_to_radar(cfg, water_mask_path, 'A')
+        water_mask_a = (water_distance_a == 0)
+
+        if iono_method in iono_method_sideband:
+            water_distance_b = project_map_to_radar(cfg, water_mask_path, 'B')
+            water_mask_b = (water_distance_b == 0)
+
     # pull parameters for polarizations
     pol_list_a = list(iono_freq_pols['A'])
     if iono_method in iono_method_sideband:
@@ -1378,6 +1407,10 @@ def run(cfg: dict, runw_hdf5: str):
                         [block_rows_data, cols_main],
                         dtype=int)
 
+                if "water" in mask_type:
+                    water_mask_a_blk = None if water_mask_a is None else \
+                        water_mask_a[row_start:row_start + block_rows_data, :]
+
                 if iono_method == 'main_diff_low_high_subband':
                     main_image = np.empty([block_rows_data, cols_main],
                                           dtype=float)
@@ -1544,6 +1577,15 @@ def run(cfg: dict, runw_hdf5: str):
                         [block_rows_data, cols_side],
                         dtype=int)
 
+                if "water" in mask_type:
+                    water_mask_a_blk = None if water_mask_a is None \
+                        else water_mask_a[row_start:row_start +
+                                          block_rows_data, :]
+
+                    water_mask_b_blk = None if water_mask_b is None \
+                        else water_mask_b[row_start:row_start +
+                                          block_rows_data, :]
+
                 with HDF5OptimizedReader(
                         name=runw_freq_a_str, mode='r',
                         libver='latest', swmr=True) as src_main_h5, \
@@ -1825,15 +1867,10 @@ def run(cfg: dict, runw_hdf5: str):
                     # boundary of the water bodies. The values 0-100 represent
                     # the distance from the coastline and values from 101-200
                     # represent the distance from inland water boundaries.
-                    water_mask_path = \
-                            cfg["dynamic_ancillary_file_group"][
-                                "water_mask_file"]
-                    water_distance = project_map_to_radar(
-                        cfg,
-                        water_mask_path,
-                        'A')
-                    mask_image = water_distance[
-                        row_start:row_start + block_rows_data, :] == 0
+                    if iono_method in iono_method_sideband:
+                        mask_image = water_mask_b_blk
+                    else:
+                        mask_image = water_mask_a_blk
                     mask_array = mask_array & mask_image
 
                 valid_area = iono_phase_obj.get_valid_area(
diff --git a/python/packages/nisar/workflows/unwrap.py b/python/packages/nisar/workflows/unwrap.py
index 2efc6ec66..3b526baa3 100644
--- a/python/packages/nisar/workflows/unwrap.py
+++ b/python/packages/nisar/workflows/unwrap.py
@@ -323,6 +323,7 @@ def run(cfg: dict, input_hdf5: str, output_hdf5: str):
 
                 if bridge_cfg['enabled']:
                     unwrapped_phase = dst_h5[unw_path][()]
+                    unwrapped_phase[mask] = 0
                     unwrapped_phase = bridge_unwrapped_phase(
                         unwrapped_phase,
                         radius=bridge_cfg['bridge_radius'],

From 729f6adf320d46e5c5549808fb65230303d4180c Mon Sep 17 00:00:00 2001
From: Jungkyo Jung <jungkyoj@nisar-adt-dev-5.jpl.nasa.gov>
Date: Sat, 11 Apr 2026 22:17:10 +0000
Subject: [PATCH 2/6] fix issue in unitttest

---
 python/packages/nisar/workflows/unwrap.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/python/packages/nisar/workflows/unwrap.py b/python/packages/nisar/workflows/unwrap.py
index 3b526baa3..d3a3a7fd7 100644
--- a/python/packages/nisar/workflows/unwrap.py
+++ b/python/packages/nisar/workflows/unwrap.py
@@ -323,7 +323,9 @@ def run(cfg: dict, input_hdf5: str, output_hdf5: str):
 
                 if bridge_cfg['enabled']:
                     unwrapped_phase = dst_h5[unw_path][()]
-                    unwrapped_phase[mask] = 0
+                    if unwrap_args["preprocess_wrapped_phase"]["enabled"]:
+                        if mask is not None:
+                            unwrapped_phase[mask] = 0
                     unwrapped_phase = bridge_unwrapped_phase(
                         unwrapped_phase,
                         radius=bridge_cfg['bridge_radius'],

From 37646188f124a0616b28905c8643b163300c7713 Mon Sep 17 00:00:00 2001
From: Seongsu Jeong <27862199+seongsujeong@users.noreply.github.com>
Date: Thu, 21 May 2026 03:25:03 -0700
Subject: [PATCH 3/6] extend the input raster's EPSG beyond 4326

---
 python/packages/isce3/unwrap/preprocess.py | 77 ++++++++++++++++++++--
 1 file changed, 72 insertions(+), 5 deletions(-)

diff --git a/python/packages/isce3/unwrap/preprocess.py b/python/packages/isce3/unwrap/preprocess.py
index a256ae900..a8cd2ebff 100644
--- a/python/packages/isce3/unwrap/preprocess.py
+++ b/python/packages/isce3/unwrap/preprocess.py
@@ -1,8 +1,9 @@
 import pathlib
 import journal
 import numpy as np
-from osgeo import gdal, osr
 
+from pyproj import Transformer
+from osgeo import gdal, osr
 from scipy.ndimage import median_filter, map_coordinates
 
 
@@ -392,6 +393,47 @@ def _read_gdal_with_bbox(input_raster, bbox, bbox_epsg=4326):
     return arr, [block_x0, block_y0, block_dx, block_dy]
 
 
+def _get_epsg_from_gdal_dataset(dataset):
+    """
+    Detect EPSG code from a GDAL dataset.
+
+    Parameters
+    ----------
+    dataset : gdal.Dataset
+        Input GDAL dataset
+
+    Returns
+    -------
+    int or None
+        EPSG code if successfully detected, None otherwise
+    """
+    proj = dataset.GetProjection()
+
+    if not proj:
+        return None
+
+    srs = osr.SpatialReference()
+    try:
+        srs.ImportFromWkt(proj)
+    except Exception:
+        return None
+
+    try:
+        srs.AutoIdentifyEPSG()
+    except Exception:
+        pass
+
+    epsg_code = srs.GetAuthorityCode(None)
+    if epsg_code is not None:
+        try:
+            return int(epsg_code)
+        except (ValueError, TypeError):
+            print(f"Warning: Failed to detect EPSG code. Dataset: {dataset.GetDescription()}, projection: {proj}")
+            return None
+
+    return None
+
+
 def _find_rdr2geo_paths(scratch_path, freq):
     """
     Find x.rdr and y.rdr files for the given frequency inside scratch_path.
@@ -479,9 +521,22 @@ def project_map_to_radar(cfg, input_data_path, freq):
     _, output_dtype = _get_gdal_raster_shape_type(input_data_path)
     geo_data_raster = gdal.Open(input_data_path)
 
+    # Determine the EPSG code from the input watermask projection
+    # The coordinate values in x.rdr and y.rdr should match this projection
+    bbox_epsg = _get_epsg_from_gdal_dataset(geo_data_raster)
+
+    if bbox_epsg is None:
+        error_channel = journal.error('unwrap.preprocess.project_map_to_radar')
+        err_str = (f"Could not determine EPSG code from input raster: "
+                   f"{input_data_path}. Please ensure the raster has valid "
+                   f"projection information.")
+        error_channel.log(err_str)
+        raise ValueError(err_str)
+
     # for both x and y rasters, decimate and get extents
     decimated_blocks = {}
     decimated_extents = {}
+
     for xy, input_path in topo_paths.items():
         # open input raster for reading
         input_data_raster = gdal.Open(input_path)
@@ -503,18 +558,31 @@ def project_map_to_radar(cfg, input_data_path, freq):
                        slice_rg_start:slice_rg_end:rg_looks]
 
         # save decimated extents and array for current axis
-        decimated_extents[xy] = [np.nanmin(decimated_arr),
-                                 np.nanmax(decimated_arr)]
+
         decimated_blocks[xy] = decimated_arr
         del input_data
 
+    # Reproject coordinates in `decimated_blocks` when `geo_data_raster` is not in 4326
+    # NOTE: transforming (5000, 5000) points took about 2.5 seconds on M1 pro, so
+    # this should not be a bottleneck.
+    if bbox_epsg != 4326:
+        transformer_4326_to_watermask = Transformer.from_crs(4326, bbox_epsg, always_xy=True)
+        decimated_blocks['x'], decimated_blocks['y'] = transformer_4326_to_watermask.transform(
+                decimated_blocks['x'], decimated_blocks['y'])
+
+    # update decimated extents after reprojection
+    for xy in ['x', 'y']:
+        decimated_extents[xy] = [np.nanmin(decimated_blocks[xy]),
+                                 np.nanmax(decimated_blocks[xy])]
+
     # get bounding for decimated extents
     bbox = [decimated_extents['x'][0], decimated_extents['y'][0],
             decimated_extents['x'][1], decimated_extents['y'][1]]
 
     # read map bounded by decimated extents of xy block
+    # Pass the detected EPSG code so bbox coordinates are interpreted correctly
     input_arr_block, [block_x0, block_y0, block_dx, block_dy] = \
-        _read_gdal_with_bbox(geo_data_raster, bbox)
+        _read_gdal_with_bbox(geo_data_raster, bbox, bbox_epsg=bbox_epsg)
 
     # prepare output array
     output_arrays = np.zeros(decimated_blocks['y'].shape,
@@ -578,4 +646,3 @@ def interpret_subswath_mask(subswath_mask, nodata=255):
     water = np.where(nd, False, water)
 
     return reference_valid, secondary_valid, water
-    
\ No newline at end of file

From f8d37c73f26ee32c41335ea01ab3cdce978b648a Mon Sep 17 00:00:00 2001
From: Seongsu Jeong <27862199+seongsujeong@users.noreply.github.com>
Date: Wed, 27 May 2026 14:12:10 -0700
Subject: [PATCH 4/6] Apply suggestion from code review

---
 python/packages/isce3/unwrap/preprocess.py    | 26 +++++++++----------
 python/packages/nisar/workflows/ionosphere.py | 15 ++++-------
 2 files changed, 17 insertions(+), 24 deletions(-)

diff --git a/python/packages/isce3/unwrap/preprocess.py b/python/packages/isce3/unwrap/preprocess.py
index a8cd2ebff..ab88651e9 100644
--- a/python/packages/isce3/unwrap/preprocess.py
+++ b/python/packages/isce3/unwrap/preprocess.py
@@ -257,13 +257,15 @@ def _read_gdal_with_bbox(input_raster, bbox, bbox_epsg=4326):
     # north-up only, same practical assumption as most GeoTIFF use cases here
     if gt[2] != 0 or gt[4] != 0:
         raise NotImplementedError(
-            "_read_gdal_with_bbox currently supports only north-up rasters."
+            "_read_gdal_with_bbox currently does not support affine transformation."
         )
 
     raster_srs = osr.SpatialReference()
     raster_srs.ImportFromWkt(proj)
 
     try:
+        # Attempt to identify and SET EPSG code for raster SRS, especially
+        # in case that the EPSG code is missing in authority info
         raster_srs.AutoIdentifyEPSG()
     except Exception:
         pass
@@ -275,7 +277,6 @@ def _read_gdal_with_bbox(input_raster, bbox, bbox_epsg=4326):
 
     xmin, ymin, xmax, ymax = map(float, bbox)
 
-    # target SRS = bbox_epsg
     dst_srs = osr.SpatialReference()
     dst_srs.ImportFromEPSG(int(bbox_epsg))
     try:
@@ -298,6 +299,11 @@ def _read_gdal_with_bbox(input_raster, bbox, bbox_epsg=4326):
         ys = [p[1] for p in corners]
         xmin_src, xmax_src = min(xs), max(xs)
         ymin_src, ymax_src = min(ys), max(ys)
+
+        # handle antimeridian crossing case for bbox in lat/lon
+        if bbox_epsg == 4326 and (xmax_src - xmin_src) > 180:
+            xmin_src, xmax_src = (xmax_src, xmin_src + 360.0)
+
     else:
         xmin_src, xmax_src = xmin, xmax
         ymin_src, ymax_src = ymin, ymax
@@ -341,17 +347,8 @@ def _read_gdal_with_bbox(input_raster, bbox, bbox_epsg=4326):
         p10 = tx_to_bbox.TransformPoint(x0 + dx, y0)[:2]
         p01 = tx_to_bbox.TransformPoint(x0, y0 + dy)[:2]
 
-        est_dx = abs(p10[0] - p00[0])
-        est_dy = abs(p01[1] - p00[1])
-
-        # fallback in case estimate becomes zero or numerically unstable
-        if not np.isfinite(est_dx) or est_dx <= 0:
-            est_dx = max((xmax - xmin) / 1000.0, 1e-6)
-        if not np.isfinite(est_dy) or est_dy <= 0:
-            est_dy = max((ymax - ymin) / 1000.0, 1e-6)
-
-        out_xres = est_dx
-        out_yres = est_dy
+        out_xres = abs(p10[0] - p00[0])
+        out_yres = abs(p01[1] - p00[1])
 
     # Warp directly to the requested bbox / requested CRS
     warped_ds = gdal.Warp(
@@ -407,6 +404,7 @@ def _get_epsg_from_gdal_dataset(dataset):
     int or None
         EPSG code if successfully detected, None otherwise
     """
+    warning_channel = journal.warning('unwrap._get_epsg_from_gdal_dataset')
     proj = dataset.GetProjection()
 
     if not proj:
@@ -428,7 +426,7 @@ def _get_epsg_from_gdal_dataset(dataset):
         try:
             return int(epsg_code)
         except (ValueError, TypeError):
-            print(f"Warning: Failed to detect EPSG code. Dataset: {dataset.GetDescription()}, projection: {proj}")
+            warning_channel.log(f"Failed to detect EPSG code. Dataset: {dataset.GetDescription()}, projection: {proj}")
             return None
 
     return None
diff --git a/python/packages/nisar/workflows/ionosphere.py b/python/packages/nisar/workflows/ionosphere.py
index a74512ec6..12f39ac4b 100644
--- a/python/packages/nisar/workflows/ionosphere.py
+++ b/python/packages/nisar/workflows/ionosphere.py
@@ -537,10 +537,8 @@ def _to_complex_if_needed(arr):
 
                     # Ensure complex (convert real-valued phase in radians to
                     # complex phase)
-                    use_complex = np.iscomplexobj(first_data_block) or np.iscomplexobj(second_data_block)
-                    if use_complex:
-                        first_data_block = _to_complex_if_needed(first_data_block)
-                        second_data_block = _to_complex_if_needed(second_data_block)
+                    first_data_block = _to_complex_if_needed(first_data_block)
+                    second_data_block = _to_complex_if_needed(second_data_block)
 
                     # Resample first dataset (frequency A / low subband) to
                     # match the grid of the second dataset (frequency B / high subband).
@@ -582,10 +580,8 @@ def _to_complex_if_needed(arr):
                         invalid = None
 
                     # Compute the differential phase
-                    if use_complex:
-                        diff_phase = first_data_block * np.conj(second_data_block)
-                    else:
-                        diff_phase = first_data_block - second_data_block
+                    diff_phase = first_data_block * np.conj(second_data_block)
+
                     if invalid is not None:
                         diff_phase[invalid] = invalid_fill_value
 
@@ -1036,8 +1032,7 @@ def run_insar_workflow(iono_insar_cfg, original_dict, out_paths,
     # decimate offsets for frequency B and create ionosphere layers
     if 'B' in iono_freq_pol:
         decimate_freq_a_offset(iono_insar_cfg, original_dict)
-    print(iono_insar_cfg['processing']['input_subset'][
-                        'list_of_frequencies'])
+
     if iono_insar_cfg['processing']['fine_resample']['enabled']:
         resample_slc_v2.run(iono_insar_cfg, 'fine')
     else:

From bf9e946dbcc3e19a1f06c50d0f6f36b3a313bd51 Mon Sep 17 00:00:00 2001
From: Seongsu Jeong <27862199+seongsujeong@users.noreply.github.com>
Date: Thu, 28 May 2026 14:16:06 -0700
Subject: [PATCH 5/6] Apply suggestion from PR review

---
 python/packages/nisar/workflows/ionosphere.py | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/python/packages/nisar/workflows/ionosphere.py b/python/packages/nisar/workflows/ionosphere.py
index 12f39ac4b..5c576ebb3 100644
--- a/python/packages/nisar/workflows/ionosphere.py
+++ b/python/packages/nisar/workflows/ionosphere.py
@@ -884,16 +884,15 @@ def insar_ionosphere_pair(original_cfg, runw_hdf5):
                                             'ionosphere',
                                             'main_diff_ms_band',
                                             'RIFG.h5')
-                unwrap_needed = True
-                diff_phase_output = pathlib.Path(diff_dir, 'RIFG.h5')
+
             else:
                 out_paths = original_out_paths
                 new_scratch = orig_scratch_path
-                phase_second = out_paths['RUNW']
+                phase_second = out_paths['RIFG']
                 additional_runw = out_paths['RUNW']
 
-                unwrap_needed = False
-                diff_phase_output = pathlib.Path(diff_dir, 'RUNW.h5')
+            diff_phase_output = pathlib.Path(diff_dir, 'RIFG.h5')
+
             iono_insar_cfg['product_path_group'][
                 'scratch_path'] = diff_dir
             iono_insar_cfg['product_path_group'][
@@ -940,7 +939,7 @@ def insar_ionosphere_pair(original_cfg, runw_hdf5):
                     dest_pol_path = f"{dest_freq_path}/interferogram/{pol_b}"
                     runw_path_b_freq = f"{dest_pol_path}/unwrappedPhase"
 
-                    second_data_path.append(runw_path_b_freq)
+                second_data_path.append(rifg_path_freq)
             second_slant_path = f"{dest_freq_path}/interferogram/slantRange"
             second_mask_path = f"{dest_freq_path}/interferogram/mask"
 
@@ -961,8 +960,8 @@ def insar_ionosphere_pair(original_cfg, runw_hdf5):
             # Since main_diff_low_high_subband method does not need to
             # unwrap low and high subband interferogram, but need to
             # unwrap the difference between low and high subband interferogram
-            if unwrap_needed:
-                unwrap.run(iono_insar_cfg, out_paths['RIFG'], out_paths['RUNW'])
+            unwrap.run(iono_insar_cfg, out_paths['RIFG'], out_paths['RUNW'])
+
     # restore original paths
     original_cfg['input_file_group']['reference_rslc_file'] = \
         partial_orig_cfg_dict['reference_rslc_file']

From 5915644ed47239d3e5eb54d40cdc069e919d401d Mon Sep 17 00:00:00 2001
From: Seongsu Jeong <27862199+seongsujeong@users.noreply.github.com>
Date: Thu, 28 May 2026 14:26:04 -0700
Subject: [PATCH 6/6] apply suggestion from PR review

---
 python/packages/nisar/workflows/ionosphere.py | 14 ++++----------
 1 file changed, 4 insertions(+), 10 deletions(-)

diff --git a/python/packages/nisar/workflows/ionosphere.py b/python/packages/nisar/workflows/ionosphere.py
index 5c576ebb3..f82949da2 100644
--- a/python/packages/nisar/workflows/ionosphere.py
+++ b/python/packages/nisar/workflows/ionosphere.py
@@ -928,18 +928,12 @@ def insar_ionosphere_pair(original_cfg, runw_hdf5):
 
             second_data_path = []
             for pol_b in pol_list_b:
-                if rerun_insar_pairs > 0:
-                    dest_freq_path = f"{swath_path}/frequencyB"
-                    dest_pol_path = f"{dest_freq_path}/interferogram/{pol_b}"
-                    rifg_path_freq = f"{dest_pol_path}/wrappedInterferogram"
-
-                    second_data_path.append(rifg_path_freq)
-                else:
-                    dest_freq_path = f"{runw_swath_path}/frequencyB"
-                    dest_pol_path = f"{dest_freq_path}/interferogram/{pol_b}"
-                    runw_path_b_freq = f"{dest_pol_path}/unwrappedPhase"
+                dest_freq_path = f"{swath_path}/frequencyB"
+                dest_pol_path = f"{dest_freq_path}/interferogram/{pol_b}"
+                rifg_path_freq = f"{dest_pol_path}/wrappedInterferogram"
 
                 second_data_path.append(rifg_path_freq)
+
             second_slant_path = f"{dest_freq_path}/interferogram/slantRange"
             second_mask_path = f"{dest_freq_path}/interferogram/mask"