Error message on stacked_dem_fill, re-write of raster_to_node for speed and accuracy -- previous refactor broke it.

Eli · Eli · commit e5c6e7fd2281 · 2025-09-07T19:22:31.000-07:00
diff --git a/schimpy/prepare_schism.py b/schimpy/prepare_schism.py
@@ -134,10 +134,19 @@ def create_hgrid(s, inputs, logger):
                 optimized_elevation = optimizer.optimize(opt_params)
                 s.mesh.nodes[:, 2] = np.negative(optimized_elevation)
             elif method=="volume_tvd":
-                dem_spec = section.get("dem_list")
-                if dem_spec is None:
+                dem_list = section.get("dem_list")
+                if dem_list is None:
                     raise ValueError("For method=volume_tvd, provide dem_list")
-
+                s.dem_list = dem_list
+                logger.info(f"Using stacked DEM fill for depth initialization. Length of DEM list: {len(dem_list)}")
+                s.mesh.nodes[:, 2] = stacked_dem_fill(
+                            dem_list,
+                            s.mesh.nodes[:, :2],
+                            inputs["prepro_output_dir"],
+                            require_all=False,
+                            na_fill=default_depth_for_missing_dem,
+                            negate=True
+                        )
                 sl = opt_params.get("shoreline")
                 shoreline = None
                 if sl is not None:
@@ -200,12 +209,12 @@ def create_hgrid(s, inputs, logger):
                 logger.info("Start optimizing the mesh ('volume_tvd').")
                 _ = refine_volume_tvd(
                     s.mesh,
-                    dem_spec=dem_spec,
+                    dem_spec=dem_list,
                     out_dir=inputs["prepro_output_dir"],
                     shoreline=shoreline,
                     floor=floor,
                     tvd=tvd,
-                    cache_dir=os.path.join(inputs["prepro_output_dir"], "dem_cache"),
+                    cache_dir=os.path.join(inputs["prepro_output_dir"], ".dem_cache"),
                     logger=logger
                 )
             else:
@@ -535,13 +544,7 @@ def process_output_dir(inputs):
         outdir = inputs["prepro_output_dir"]
         force = True
     else:
-        warnings.warn(
-            "No output_dir specification. This will not be allowed in the future. \n"
-            + "Using '.' but specification of a separate directory recommended. ",
-            FutureWarning,
-        )
-        outdir = "."
-        inputs["prepro_output_dir"] = outdir
+        raise ValueError("Inplace preprocessing not allowed. prepro_output_dir must be specified in launch yaml file. \nConvention is 'prepro_out'")
     if os.path.exists(outdir):
         created = False
     else:
diff --git a/schimpy/raster_to_nodes.py b/schimpy/raster_to_nodes.py
@@ -1,12 +1,42 @@
-# -*- coding: utf-8 -*-
-"""Optimized raster_to_nodes function without rasterstats, using rasterio directly."""
-
 import numpy as np
-import rasterio
-from rasterio.mask import mask
-from shapely.geometry import Polygon, mapping
 import itertools
+import rasterio
+from rasterio.features import rasterize
+from rasterio.windows import from_bounds
+from shapely.geometry import Polygon, box
+
+# ---------- helpers (minimal) ----------
+
+def _union_missing_mask(ma, extra_mask_values):
+    """Union raster mask with explicit sentinels (e.g., 0.0, -9999)."""
+    if extra_mask_values is None:
+        return ma
+    if not isinstance(extra_mask_values, (list, tuple, np.ndarray)):
+        extra_mask_values = [extra_mask_values]
+    data = ma.data
+    add_mask = np.zeros(data.shape, dtype=bool)
+    for mv in extra_mask_values:
+        add_mask |= (data == mv)
+    return np.ma.array(data, mask=(np.asarray(ma.mask, bool) | add_mask), copy=False)
+
+def _elements_bounds_in_raster(elements_polys, raster_bounds):
+    """BBox of union(elements ∩ raster) or None if empty."""
+    rb = box(*raster_bounds)
+    any_hit = False
+    xmin, ymin, xmax, ymax = +np.inf, +np.inf, -np.inf, -np.inf
+    for poly in elements_polys:
+        if poly.is_empty: 
+            continue
+        inter = poly.intersection(rb)
+        if inter.is_empty: 
+            continue
+        bxmin, bymin, bxmax, bymax = inter.bounds
+        xmin, ymin = min(xmin, bxmin), min(ymin, bymin)
+        xmax, ymax = max(xmax, bxmax), max(ymax, bymax)
+        any_hit = True
+    return (xmin, ymin, xmax, ymax) if any_hit else None
 
+# ---------- main (rasterize-only) ----------
 
 def raster_to_nodes(
     mesh,
@@ -17,116 +47,192 @@ def raster_to_nodes(
     mask_value=None,
     fill_value=None,
     band=1,
+    missing_policy="exclude",
+    rasterize_all_touched=True,
 ):
     """
-    Applies raster values to mesh by calculating means over surrounding elements.
+    Map raster values onto mesh nodes by averaging values over surrounding elements.
+
+    This function samples a raster onto an unstructured mesh. Each element polygon
+    is intersected with the raster grid, values are averaged per element, and then
+    element values are aggregated to nodes using an area-weighted mean.
 
     Parameters
     ----------
-    (same as original)
+    mesh : object
+        Mesh object with the following interface:
+          - ``mesh.nodes`` : array of node coordinates (N, >=2), x and y in the first two columns.
+          - ``mesh.elem(i)`` : return node indices for element ``i``.
+          - ``mesh.get_elems_i_from_node(node_i)`` : return element indices connected to node ``node_i``.
+          - ``mesh.areas()`` : array of element areas indexed by element id.
+    nodes_sel : sequence of int
+        Node indices at which values will be computed.
+    path_raster : str
+        Path to a raster file (e.g., GeoTIFF) readable by rasterio.
+    bins : sequence of float, optional
+        Bin edges used to classify raster values. If provided, the raster values
+        are digitized into bins before aggregation. Length must be ``M``.
+    mapped_values : sequence of float, optional
+        Values assigned to each bin index. Must have length ``M+1`` if ``bins`` is provided.
+        Typically used to map classes into the interval [0, 1].
+    mask_value : scalar or sequence of scalars, optional
+        Additional raster values to treat as missing (in addition to raster
+        NoData). Example: ``[0.0, -9999]``.
+    fill_value : float, optional
+        Default value for elements with no valid raster pixels in the continuous
+        (non-classified) case. If omitted, defaults to 0.0. Ignored when
+        classification is active.
+    band : int, default=1
+        1-based band index in the raster to sample.
+    missing_policy : {"exclude", "as_zero"}, default="exclude"
+        Policy for handling missing pixels in the classified case:
+        
+        - ``"exclude"`` : compute class averages using only valid pixels. Elements
+          with no valid pixels receive the default class value (0.0).
+        - ``"as_zero"`` : missing pixels count as class 0 by including them in the
+          denominator. This dilutes averages toward 0 when coverage is sparse.
+          
+        Has no effect in the continuous (non-classified) case.
+    rasterize_all_touched : bool, default=True
+        If True, count all raster cells touched by element polygons. If False,
+        count only cells whose centers fall within the polygons.
 
     Returns
     -------
-    numpy.array
-        means of raster values of the element balls around the nodes
+    ndarray of float
+        Array of values at each node in ``nodes_sel``. Length matches ``len(nodes_sel)``.
+
+    Notes
+    -----
+    - Continuous mode (no ``bins``): averages raster values directly.
+    - Classified mode (with ``bins`` and ``mapped_values``): raster values are
+      digitized into bins, mapped to user-provided class values, then averaged.
+    - All aggregation is area-weighted: element values are computed from raster
+      pixels, and node values are weighted by connected element areas.
     """
 
-    # Step 1: Prepare binning if needed
-    classify_raster = False
-    if bins is not None:
+    classify = bins is not None
+    if classify:
         if mapped_values is None:
             raise ValueError("mapped_values must be provided if bins are used.")
         if len(mapped_values) != len(bins) + 1:
             raise ValueError("mapped_values must be one longer than bins.")
-        classify_raster = True
+        if missing_policy not in ("exclude", "as_zero"):
+            raise ValueError("missing_policy must be 'exclude' or 'as_zero'.")
 
-    # Step 2: Precompute node balls (cache)
-    elements_in_balls = {
-        node_i: list(mesh.get_elems_i_from_node(node_i)) for node_i in nodes_sel
-    }
-    elements_in_polygon = sorted(
-        set(itertools.chain.from_iterable(elements_in_balls.values()))
-    )
+    # Node -> elements and unique element ids
+    elements_in_balls = {n: list(mesh.get_elems_i_from_node(n)) for n in nodes_sel}
+    elements_all = sorted(set(itertools.chain.from_iterable(elements_in_balls.values())))
 
-    # Step 3: Build element polygons
+    # Build element polygons (XY only)
     element_polygons = {}
-    for elem_i in elements_in_polygon:
+    for elem_i in elements_all:
         nodes_idx = mesh.elem(elem_i)
-        coords = mesh.nodes[nodes_idx, :2]  # x,y only
+        coords = mesh.nodes[nodes_idx, :2]
         element_polygons[elem_i] = Polygon(coords)
 
-    # Step 4: Open raster
+    # Defaults per path
+    elem_default_cont  = (fill_value if fill_value is not None else 0.0)
+    elem_default_class = 0.0
+
     with rasterio.open(path_raster) as src:
-        nodata = src.nodatavals[band - 1] if src.nodatavals else None
-
-        # Read full raster band if binning is needed
-        if classify_raster:
-            raster_array = src.read(band)
-            if mask_value is not None:
-                raster_array = np.where(
-                    raster_array == mask_value, fill_value, raster_array
-                )
-            digitized = np.digitize(raster_array, bins, right=False)
-            classified_array = np.array(mapped_values)[digitized]
-        else:
-            classified_array = None  # Will read on demand
-
-        # Step 5: Calculate mean per element
-        sav_in_elements = {}
-        for elem_i, polygon in element_polygons.items():
-            geom = [mapping(polygon)]
-            if classify_raster is False:
-                out_image, out_transform = mask(
-                    src, geom, crop=True, indexes=band, nodata=nodata
-                )
-                data = out_image[0]
+        # Minimal window that covers mesh ∩ raster
+        bbox = _elements_bounds_in_raster(list(element_polygons.values()), src.bounds)
+        # If nothing overlaps: return defaults
+        if bbox is None:
+            default_val = elem_default_class if classify else elem_default_cont
+            return np.full(len(nodes_sel), default_val, dtype=float)
+
+        win = from_bounds(*bbox, transform=src.transform)
+        full_win = from_bounds(*src.bounds, transform=src.transform)
+        win = win.intersection(full_win)
+
+        # Read data as masked array over that window
+        data = src.read(band, window=win, masked=True)  # (H,W) masked
+        if data.ndim == 3:  # (1,H,W) -> (H,W) if env returns 3D
+            data = data[0]
+        data = _union_missing_mask(data, mask_value)
+
+        # Rasterize element IDs on same window grid
+        window_transform = rasterio.windows.transform(win, src.transform)
+        shapes = [(poly, int(eid)) for eid, poly in element_polygons.items()]
+        labels = rasterize(
+            shapes,
+            out_shape=data.shape,
+            transform=window_transform,
+            fill=0,           # background label
+            dtype="int32",
+            all_touched=rasterize_all_touched,
+        )
+
+        # Prepare result per-element
+        elem_values = {eid: (elem_default_class if classify else elem_default_cont)
+                       for eid in element_polygons.keys()}
+
+        # VALID pixels only
+        valid = ~data.mask
+        labs_valid = labels[valid]
+        vals_valid = data.data[valid]
+        sel_valid = labs_valid != 0
+
+        if classify:
+            # If there are valid class pixels, aggregate them
+            if sel_valid.any():
+                bins_arr = np.asarray(bins)
+                mapped = np.asarray(mapped_values)
+                idx = np.digitize(vals_valid[sel_valid], bins_arr, right=False)  # 0..len(bins)
+                class_vals = mapped[idx]
+
+                max_id = int(labels.max()) if labels.size else 0
+                sums_valid = np.bincount(labs_valid[sel_valid], weights=class_vals, minlength=max_id + 1)
+                cnts_valid = np.bincount(labs_valid[sel_valid], minlength=max_id + 1)
             else:
-                # Clip manually from classified_array
-                bounds = polygon.bounds  # (minx, miny, maxx, maxy)
-                row_min, col_min = src.index(bounds[0], bounds[3])  # (xmin, ymax)
-                row_max, col_max = src.index(bounds[2], bounds[1])  # (xmax, ymin)
-                rows = slice(min(row_min, row_max), max(row_min, row_max) + 1)
-                cols = slice(min(col_min, col_max), max(col_min, col_max) + 1)
-                if rows.start == rows.stop or cols.start == cols.stop:
-                    sav_in_elements[elem_i] = 0.0
-                    continue
-                window = classified_array[rows, cols]
-                # Verify that window has nonzero size
-                if window.shape[0] == 0 or window.shape[1] == 0:
-                    sav_in_elements[elem_i] = 0.0
-                    continue
-                # Build temporary profile
-                transform = src.window_transform(
-                    ((rows.start, rows.stop), (cols.start, cols.stop))
-                )
-                with rasterio.io.MemoryFile() as memfile:
-                    with memfile.open(
-                        driver="GTiff",
-                        height=window.shape[0],
-                        width=window.shape[1],
-                        count=1,
-                        dtype=window.dtype,
-                        transform=transform,
-                        crs=src.crs,
-                        nodata=-999,
-                    ) as dataset:
-                        dataset.write(window, 1)
-                        out_image, out_transform = mask(
-                            dataset, geom, crop=True, indexes=1, nodata=-999
-                        )
-                        data = out_image[0]
-
-            masked = np.ma.masked_array(data, mask=(data == nodata))
-            mean_val = masked.mean() if masked.count() > 0 else 0.0
-            sav_in_elements[elem_i] = mean_val
-
-    # Step 6: Assemble node values
+                max_id = int(labels.max()) if labels.size else 0
+                sums_valid = np.zeros(max_id + 1, dtype=float)
+                cnts_valid = np.zeros(max_id + 1, dtype=float)
+
+            if missing_policy == "as_zero":
+                # Denominator: all pixels of element (valid+masked), excluding background
+                labs_all = labels.ravel()
+                sel_all = labs_all != 0
+                cnts_total = np.bincount(labs_all[sel_all], minlength=max_id + 1)
+                present = np.nonzero(cnts_total)[0]
+                for eid in present:
+                    denom = cnts_total[eid]
+                    num = sums_valid[eid]  # masked contribute 0
+                    elem_values[eid] = float(num / denom) if denom > 0 else elem_default_class
+            else:
+                # "exclude": denominator is count of valid only
+                present = np.nonzero(cnts_valid)[0]
+                for eid in present:
+                    denom = cnts_valid[eid]
+                    elem_values[eid] = float(sums_valid[eid] / denom) if denom > 0 else elem_default_class
+
+        else:
+            # Continuous: mean over valid pixels per element
+            if sel_valid.any():
+                max_id = int(labels.max()) if labels.size else 0
+                sums = np.bincount(labs_valid[sel_valid], weights=vals_valid[sel_valid], minlength=max_id + 1)
+                cnts = np.bincount(labs_valid[sel_valid], minlength=max_id + 1)
+                present = np.nonzero(cnts)[0]
+                for eid in present:
+                    elem_values[eid] = float(sums[eid] / cnts[eid]) if cnts[eid] > 0 else elem_default_cont
+            # else: keep defaults (no valid pixels anywhere)
+
+    # Area-weighted average from elements to nodes
     elem_areas = mesh.areas()
-    sav_at_nodes = np.empty((len(nodes_sel),), dtype=float)
+    out_vals = np.empty((len(nodes_sel),), dtype=float)
     for i, node_i in enumerate(nodes_sel):
         ball = elements_in_balls[node_i]
-        values = [sav_in_elements[e] for e in ball]
-        weights = elem_areas[ball]
-        sav_at_nodes[i] = np.average(values, weights=weights)
+        if not ball:
+            out_vals[i] = 0.0
+            continue
+        values = np.array([elem_values[e] for e in ball], dtype=float)
+        weights = np.array(elem_areas[ball], dtype=float)
+        bad = ~np.isfinite(values)
+        if bad.any():
+            weights = weights.copy(); weights[bad] = 0.0
+            values = np.where(bad, 0.0, values)
+        out_vals[i] = np.average(values, weights=weights) if weights.sum() > 0 else 0.0
 
-    return sav_at_nodes
+    return out_vals
diff --git a/schimpy/stacked_dem_fill.py b/schimpy/stacked_dem_fill.py
@@ -64,10 +64,11 @@ def _stacked_dem_fill(
         f_out = ensure_outdir(out_dir, "dem_misses.txt")
         miss_idx = np.isnan(values)
         miss_points = points[miss_idx]
+        nmiss = miss_points.shape[0]
         f = open(f_out, "w")
         for p in miss_points:
             f.write("%s,%s\n" % tuple(p))
-            message = "DEMs provided do not cover all the points. See file dem_misses.txt for locations"
+            message = f"DEMs provided do not cover {nmiss} points. See file dem_misses.txt for locations and to check for significance"
         f.close()
 
         gdf = gpd.GeoDataFrame(
