Merge branch 'dswx-ni-calval' into dswx_ni_masking_fast

Author: oberonia78

src/dswx_sar/common/_dswx_geogrid.py

@@ -1,6 +1,7 @@
 
 from dataclasses import dataclass
-
+from typing import Optional, Tuple
+import math
 import numpy as np
 from osgeo import osr, gdal
 
@@ -102,45 +103,194 @@ def from_geotiff(cls, geotiff_path):
         return cls(start_x, start_y, end_x, end_y, spacing_x, spacing_y,
                    length, width, epsg)
 
-    def update_geogrid(self, geotiff_path):
+    def update_geogrid(self, geotiff_path: str) -> None:
+        """
+        Update this geogrid to be the union of itself and the geogrid from geotiff_path.
+        Works for both north-up rasters (spacing_y < 0) and south-up (spacing_y > 0).
+
+        Conventions:
+        - We union using true bounds (xmin/xmax/ymin/ymax) independent of sign.
+        - We keep this object's spacing sign if already set; otherwise adopt new's.
+        - start_x/start_y represent the "origin corner" implied by spacing signs:
+            * spacing_x > 0  => start_x is xmin
+            * spacing_x < 0  => start_x is xmax
+            * spacing_y < 0  => start_y is ymax (north-up)
+            * spacing_y > 0  => start_y is ymin (south-up)
+        """
+        new = DSWXGeogrid.from_geotiff(geotiff_path)
+
+        # EPSG check
+        if (not np.isnan(self.epsg)) and int(self.epsg) != int(new.epsg):
+            raise ValueError(
+                f"EPSG codes do not match: existing={self.epsg}, new={new.epsg}"
+            )
+
+        # Decide spacings (magnitude + sign)
+        # If self spacing is unset, inherit from new; otherwise keep self.
+        if np.isnan(self.spacing_x):
+            self.spacing_x = new.spacing_x
+        if np.isnan(self.spacing_y):
+            self.spacing_y = new.spacing_y
+
+        # Helper: bounds independent of sign
+        def _bounds(g):
+            xmin = min(g.start_x, g.end_x)
+            xmax = max(g.start_x, g.end_x)
+            ymin = min(g.start_y, g.end_y)
+            ymax = max(g.start_y, g.end_y)
+            return xmin, ymin, xmax, ymax
+
+        # Union bounds (handle NaNs gracefully)
+        xmin0, ymin0, xmax0, ymax0 = _bounds(self) if not np.isnan(self.start_x) else (np.nan, np.nan, np.nan, np.nan)
+        xmin1, ymin1, xmax1, ymax1 = _bounds(new)
+
+        xs = [v for v in (xmin0, xmin1) if not np.isnan(v)]
+        ys = [v for v in (ymin0, ymin1) if not np.isnan(v)]
+        xe = [v for v in (xmax0, xmax1) if not np.isnan(v)]
+        ye = [v for v in (ymax0, ymax1) if not np.isnan(v)]
+
+        xmin_u = min(xs)
+        ymin_u = min(ys)
+        xmax_u = max(xe)
+        ymax_u = max(ye)
+
+        # Reconstruct start/end consistent with spacing sign convention
+        sx = self.spacing_x
+        sy = self.spacing_y
+
+        # X
+        if sx >= 0:
+            self.start_x = xmin_u
+            self.end_x   = xmax_u
+        else:
+            self.start_x = xmax_u
+            self.end_x   = xmin_u
+
+        # Y
+        if sy < 0:  # north-up
+            self.start_y = ymax_u
+            self.end_y   = ymin_u
+        else:       # south-up or unusual
+            self.start_y = ymin_u
+            self.end_y   = ymax_u
+
+        # width/length (ensure positive)
+        if not np.isnan(sx) and sx != 0:
+            self.width = int(round((self.end_x - self.start_x) / sx))
+            self.width = abs(self.width)
+
+        if not np.isnan(sy) and sy != 0:
+            self.length = int(round((self.end_y - self.start_y) / sy))
+            self.length = abs(self.length)
+
+        # EPSG
+        if np.isnan(self.epsg):
+            self.epsg = new.epsg
+
+
+    @staticmethod
+    def _snap_bounds_to_spacing(
+        xmin: float, ymin: float, xmax: float, ymax: float,
+        spacing_x: float, spacing_y: float,
+    ) -> Tuple[float, float, float, float]:
+        """
+        Snap bounds to pixel grid so width/length become integers.
+        Uses spacing magnitude; preserves original spacing sign when writing geogrid.
+        """
+        sx = abs(spacing_x)
+        sy = abs(spacing_y)
+        if sx == 0 or sy == 0 or math.isnan(sx) or math.isnan(sy):
+            raise ValueError("Invalid spacing for snapping bounds.")
+
+        # snap outward (cover the requested bbox)
+        xmin_s = math.floor(xmin / sx) * sx
+        ymin_s = math.floor(ymin / sy) * sy
+        xmax_s = math.ceil(xmax / sx) * sx
+        ymax_s = math.ceil(ymax / sy) * sy
+        return xmin_s, ymin_s, xmax_s, ymax_s
+
+    def clip_to_bbox(
+        self,
+        bbox: Tuple[float, float, float, float],
+        bbox_epsg: Optional[int] = None,
+        snap: bool = True,
+        snap_outward: bool = True,  # keep outward snap as default
+    ) -> None:
         """
-        Update the existing geographical grid parameters based on a new
-        GeoTIFF file, extending the grid to encompass both the existing
-        and new grid areas.
+        Intersect current geogrid with bbox (must be in same epsg).
+        Updates start/end/width/length to the intersection.
         """
-        new_geogrid = DSWXGeogrid.from_geotiff(geotiff_path)
-
-        if self.epsg != new_geogrid.epsg and not np.isnan(self.epsg):
-            raise ValueError("EPSG codes of the existing and "
-                             "new geogrids do not match.")
-        self.start_x = min(filter(lambda x: not np.isnan(x),
-                                  [self.start_x, new_geogrid.start_x]))
-        self.end_x = max(filter(lambda x: not np.isnan(x),
-                                [self.end_x, new_geogrid.end_x]))
-
-        if self.spacing_y > 0 or np.isnan(self.spacing_y):
-            self.end_y = max(filter(lambda x: not np.isnan(x),
-                                    [self.end_y, new_geogrid.end_y]))
-            self.start_y = min(filter(lambda x: not np.isnan(x),
-                                      [self.start_y, new_geogrid.start_y]))
+        xmin_b, ymin_b, xmax_b, ymax_b = bbox
+
+        epsg_ok = not (self.epsg is None or (isinstance(self.epsg, float) and np.isnan(self.epsg)))
+        if bbox_epsg is not None and epsg_ok and int(self.epsg) != int(bbox_epsg):
+            raise ValueError(f"EPSG mismatch: geogrid epsg={self.epsg}, bbox epsg={bbox_epsg}. "
+                             "Reproject bbox first or force EPSG earlier.")
+
+        xmin_g, ymin_g, xmax_g, ymax_g = self.bounds()
+
+        ixmin = max(xmin_g, xmin_b)
+        iymin = max(ymin_g, ymin_b)
+        ixmax = min(xmax_g, xmax_b)
+        iymax = min(ymax_g, ymax_b)
+
+        if ixmin >= ixmax or iymin >= iymax:
+            if ixmin >= ixmax:
+                print('error x', ixmin, ixmax)
+            if iymin >= iymax:
+                print('error y', iymin, iymax)
+            raise ValueError("DB bbox does not intersect the data extent (geogrid).")
+
+        if snap:
+            # snap outward so you don't accidentally clip requested area by <1 pixel
+            ixmin, iymin, ixmax, iymax = self._snap_bounds_to_spacing(
+                ixmin, iymin, ixmax, iymax, self.spacing_x, self.spacing_y
+            )
+
+        # Re-apply sign convention used in your geogrid:
+        # spacing_x typically +, spacing_y may be negative.
+        sx = self.spacing_x
+        sy = self.spacing_y
+
+        # start_x should be left edge
+        self.start_x = ixmin
+        self.end_x = ixmax
+
+        # start_y should match your spacing_y sign convention:
+        # If sy < 0, start_y is the TOP (max y). If sy > 0, start_y is BOTTOM (min y).
+        if sy < 0:
+            self.start_y = iymax
+            self.end_y = iymin
         else:
-            self.start_y = max(filter(lambda x: not np.isnan(x),
-                                      [self.start_y, new_geogrid.start_y]))
-            self.end_y = min(filter(lambda x: not np.isnan(x),
-                                    [self.end_y, new_geogrid.end_y]))
-
-        self.spacing_x = new_geogrid.spacing_x \
-            if not np.isnan(new_geogrid.spacing_x) else self.spacing_x
-        self.spacing_y = new_geogrid.spacing_y \
-            if not np.isnan(new_geogrid.spacing_y) else self.spacing_y
-
-        if not np.isnan(self.start_x) and not np.isnan(self.end_x) and \
-           not np.isnan(self.spacing_x):
-            self.width = int((self.end_x - self.start_x) / self.spacing_x)
-
-        if not np.isnan(self.start_y) and not np.isnan(self.end_y) and \
-           not np.isnan(self.spacing_y):
-            self.length = int((self.end_y - self.start_y) / self.spacing_y)
-
-        self.epsg = new_geogrid.epsg \
-            if not np.isnan(new_geogrid.epsg) else self.epsg
+            self.start_y = iymin
+            self.end_y = iymax
+
+        self.width  = int(round((self.end_x - self.start_x) / self.spacing_x)) if sx != 0 else self.width
+        self.length = int(round((self.end_y - self.start_y) / self.spacing_y)) if sy != 0 else self.length
+
+    def bounds(self) -> Tuple[float, float, float, float]:
+        """
+        Return geogrid bounds as (xmin, ymin, xmax, ymax) in the geogrid EPSG.
+        Works regardless of spacing_y sign (north-up geotiffs usually have spacing_y < 0).
+        """
+        # Ensure initialized
+        required = {
+            "start_x": self.start_x,
+            "start_y": self.start_y,
+            "end_x": self.end_x,
+            "end_y": self.end_y,
+            "spacing_x": self.spacing_x,
+            "spacing_y": self.spacing_y,
+        }
+        missing = [k for k, v in required.items() if v is None or (isinstance(v, float) and np.isnan(v))]
+        if missing:
+            raise RuntimeError(
+                f"Geogrid not initialized; missing {missing}. "
+                "Call update_geogrid()/from_geotiff() before clip_to_bbox()."
+            )
+
+        xmin = min(self.start_x, self.end_x)
+        xmax = max(self.start_x, self.end_x)
+        ymin = min(self.start_y, self.end_y)
+        ymax = max(self.start_y, self.end_y)
+        return xmin, ymin, xmax, ymax

src/dswx_sar/common/_dswx_sar_util.py

@@ -12,6 +12,7 @@
 from osgeo import gdal, osr, ogr
 from pyproj import Transformer
 from scipy.signal import convolve2d
+from scipy.ndimage import uniform_filter
 
 gdal.DontUseExceptions()
 
@@ -1905,27 +1906,31 @@ def _calculate_output_bounds(geotransform,
     list
         Adjusted bounding box coordinates [x_min, y_min, x_max, y_max].
     """
-    x_min = geotransform[0]
-    x_max = x_min + width * geotransform[1] 
+    gt0, gt1, gt2, gt3, gt4, gt5 = geotransform
 
-    if geotransform[5] < 0:
-        y_max = geotransform[3]
-        y_min = y_max + length * geotransform[5]
-    else:
-        y_min = geotransform[3]
-        y_max = y_min + length * geotransform[5]
+    # Corner coordinates derived from geotransform + raster size
+    # (no rotation assumed; if rotation exists, this needs a different approach)
+    x0 = gt0
+    y0 = gt3
+    x1 = gt0 + width  * gt1
+    y1 = gt3 + length * gt5
 
-    x_diff = x_max - x_min
-    y_diff = y_max - y_min
+    xmin = min(x0, x1)
+    xmax = max(x0, x1)
+    ymin = min(y0, y1)
+    ymax = max(y0, y1)
 
-    if x_diff % output_spacing != 0:
-        x_max = x_min + (x_diff // output_spacing + 1) * output_spacing
+    s = float(abs(output_spacing))
+    if not np.isfinite(s) or s <= 0:
+        raise ValueError(f"Invalid output_spacing: {output_spacing}")
 
-    output_spacing = -1 * np.abs(output_spacing)
-    if y_diff % output_spacing != 0:
-        y_min = y_max + (y_diff // np.abs(output_spacing) + 1) * output_spacing
+    # Snap outward so the bounds fully cover the original area
+    xmin_s = math.floor(xmin / s) * s
+    ymin_s = math.floor(ymin / s) * s
+    xmax_s = math.ceil (xmax / s) * s
+    ymax_s = math.ceil (ymax / s) * s
 
-    return [x_min, y_min, x_max, y_max]
+    return [xmin_s, ymin_s, xmax_s, ymax_s]
 
 
 def _perform_warp_in_memory(input_file,
@@ -2114,6 +2119,37 @@ def _aggregate_10m_to_30m_conv(image, ratio, normalize_flag):
 
     return aggregated_data
 
+
+def _aggregate_10m_to_30m_fast(image: np.ndarray, ratio: int, normalize_flag: bool):
+    """
+    Fast box-sum aggregation + stride sampling.
+    Matches: convolve2d(image, ones, mode='same') then [ratio//2::ratio, ratio//2::ratio]
+    """
+    # Make sure we are working in float for normalization robustness
+    img = image.astype(np.float32, copy=False)
+
+    # uniform_filter gives the mean; multiply by area to get sum
+    area = float(ratio * ratio)
+    summed = uniform_filter(img, size=ratio, mode="constant", cval=0.0) * area
+
+    s = ratio // 2
+    aggregated = summed[s::ratio, s::ratio]
+
+    if normalize_flag:
+        valid = (img > 0).astype(np.float32, copy=False)
+        count = uniform_filter(valid, size=ratio, mode="constant", cval=0.0) * area
+        count = count[s::ratio, s::ratio]
+
+        # avoid divide-by-zero
+        aggregated = np.divide(
+            aggregated, count,
+            out=np.zeros_like(aggregated, dtype=np.float32),
+            where=(count > 0)
+        )
+
+    return aggregated
+
+
 def majority_element(num_list):
     """
     Determine the majority element in a list

src/dswx_sar/common/_filter_SAR.py

@@ -358,7 +358,6 @@ def tv_bregman(X: np.ndarray, **kwargs) -> np.ndarray:
     https://scikit-image.org/docs/stable/api/skimage.restoration.html#skimage.restoration.denoise_tv_bregman
     """
     lamb = kwargs.get('lambda_value', -1)
-    print('bregman', kwargs)
     X_db = np.log10(X, out=np.full(X.shape, np.nan), where=(~np.isnan(X)))
     X_db[np.isnan(X)] = -30
     X_db_dspkl = denoise_tv_bregman(X_db, weight=lamb)

src/dswx_sar/common/_metadata.py

@@ -11,6 +11,7 @@
 from dswx_sar.common._dswx_sar_util import (band_assign_value_dict,
                                            read_geotiff)
 from dswx_sar.common.gcov_reader import RTCReader
+from dswx_sar.common.read_h5_s3 import open_h5
 
 # Constants
 UNKNOWN = 'UNKNOWN'
@@ -511,7 +512,7 @@ def collect_frame_id(h5_list):
     frame_id_list = []
     frame_path = '/science/LSAR/identification/frameNumber'
     for rtc_file in h5_list:
-        with h5py.File(rtc_file) as src:
+        with open_h5(rtc_file) as src:
             frame = src[frame_path][()]
             frame_id_list.append(frame)
     return list(set(frame_id_list))
@@ -543,7 +544,7 @@ def count_rfi_frames(h5_list, pol_list, rfi_likelihood_thresh):
     for input_idx, rtc_file in enumerate(h5_list):
         rfi_found_in_frame = False
 
-        with h5py.File(rtc_file) as src:
+        with open_h5(rtc_file) as src:
             if freq_path_list not in src:
                 # if listOfFrequencies missing, we can't evaluate; skip this file
                 continue