Fix a bug in projecting water to cryosphere radargrid

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
+

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,