add decay_pixels arg to api/cli; coastal_aware_blend to merge FES and VDATUM, only decay inland

Author: matth-love

src/transformez/api.py

@@ -61,9 +61,10 @@ def generate_grid(
     increment: Union[str, float],
     datum_in: str,
     datum_out: str,
+    decay_pixels: int = 100,
     out_fn: Optional[str] = None,
     cache_dir: Optional[str] = None,
-    verbose: bool = False
+    verbose: bool = False,
 ) -> Optional[np.ndarray]:
     """Generate a vertical shift grid for a specific region.
 
@@ -72,6 +73,7 @@ def generate_grid(
         increment: Resolution (e.g., '3s' or 0.0008333).
         datum_in: Source datum (e.g., 'mllw', '5703').
         datum_out: Target datum (e.g., '4979', '6319').
+        decay_pixels: Set the pixel decay in case extrapolation is required.
         out_fn: If provided, saves the grid to this file (.tif or .gtx).
         cache_dir: Path to store downloaded grids.
         verbose: Enable debug logging.
@@ -109,6 +111,7 @@ def generate_grid(
         nx=nx, ny=ny,
         epsg_in=epsg_in, epsg_out=epsg_out,
         geoid_in=geoid_in, geoid_out=geoid_out,
+        decay_pixels=decay_pixels,
         cache_dir=cache_dir,
         verbose=verbose
     )
@@ -130,6 +133,7 @@ def transform_raster(
     input_raster: str,
     datum_in: str,
     datum_out: str,
+    decay_pixels: int = 100,
     output_raster: Optional[str] = None,
     cache_dir: Optional[str] = None,
     verbose: bool = False
@@ -141,6 +145,7 @@ def transform_raster(
         datum_in: Source datum of the DEM.
         datum_out: Target datum for the output DEM.
         output_raster: Path to save the transformed DEM. If None, auto-generates a name.
+        decay_pixels: Set the pixel decay in case extrapolation is required.
         cache_dir: Path to store downloaded grids.
         verbose: Enable debug logging.
 
@@ -177,6 +182,7 @@ def transform_raster(
         nx=nx, ny=ny,
         epsg_in=epsg_in, epsg_out=epsg_out,
         geoid_in=geoid_in, geoid_out=geoid_out,
+        decay_pixels=decay_pixels,
         cache_dir=cache_dir,
         verbose=verbose
     )

src/transformez/cli.py

@@ -100,6 +100,10 @@ def transformez_cli():
     grp_out.add_argument('-o', '--output', help='Output filename (default: auto-named).')
     grp_out.add_argument('--preview', action='store_true', help='Show plot of shift grid before saving.')
 
+    grp_proc = parser.add_argument_group('Process')
+    grp_proc.add_argument('--decay-pixels', type=int, default=100,
+                          help='Number of pixels to decay tidal shifts inland (default: 100). Set to 0 for strict VDatum boundaries.')
+
     grp_sys = parser.add_argument_group('System')
     grp_sys.add_argument('--list-datums', action='store_true', help='List supported datums.')
     grp_sys.add_argument('--download-htdp', action='store_true', help='Download the NGS HTDP software to the current directory.')
@@ -202,6 +206,7 @@ def transformez_cli():
         nx=nx, ny=ny,
         epsg_in=epsg_in, epsg_out=epsg_out,
         geoid_in=geoid_in, geoid_out=geoid_out,
+        decay_pixels=args.decay_pixels,
         cache_dir=cache_dir,
         verbose=args.verbose
     )

src/transformez/grid_engine.py

@@ -45,7 +45,8 @@ def plot_grid(grid_array, region, title="Vertical Shift Preview"):
     plt.figure(figsize=(10, 6))
     plot_region = [region.xmin, region.xmax, region.ymin, region.ymax]
 
-    im = plt.imshow(masked_data, extent=plot_region, cmap="RdBu_r", origin="upper")
+    # im = plt.imshow(masked_data, extent=plot_region, cmap="RdBu_r", origin="upper")
+    im = plt.imshow(masked_data, extent=plot_region, cmap="viridis", origin="upper")
     cbar = plt.colorbar(im)
     cbar.set_label("Vertical Shift (meters)")
     plt.title(title)
@@ -110,14 +111,71 @@ def load_and_interpolate(source_files, target_region, nx, ny, decay_pixels=100):
                 continue
 
         # Fill inland areas (decaying to 0) before we clear the remaining NaNs
-        mosaic = GridEngine.fill_nans(mosaic, decay_pixels=decay_pixels)
-        mosaic[np.isnan(mosaic)] = 0.0
+        # mosaic = GridEngine.fill_nans(mosaic, decay_pixels=decay_pixels)
+        # mosaic[np.isnan(mosaic)] = 0.0
 
         return mosaic
 
     @staticmethod
-    def fill_nans(data, decay_pixels=100):
-        """Fill NaNs with the nearest valid value, decayed to zero over distance."""
+    def smart_blend(in_grid, background_grid, blend_pixels=50):
+        """Smoothly blends the grid into a background grid."""
+
+        mask = np.isnan(in_grid)
+
+        if not mask.any():
+            return in_grid
+
+        if mask.all():
+            return background_grid
+
+        dist = ndimage.distance_transform_edt(mask)
+        alpha = np.clip(dist / blend_pixels, 0.0, 1.0)
+
+        nearest_indices = ndimage.distance_transform_edt(mask, return_distances=False, return_indices=True)
+        extended_vdatum = in_grid.copy()
+        extended_vdatum[mask] = in_grid[tuple(nearest_indices)][mask]
+
+        blended_data = (extended_vdatum * (1.0 - alpha)) + (background_grid * alpha)
+
+        return blended_data
+
+    @staticmethod
+    def coastal_aware_composite(vdatum_grid, global_grid, decay_pixels=100, buffer_pixels=10, max_discontinuity=0.5):
+        """Intelligently handles inland decay vs. offshore blending, while
+        filtering out low-resolution global artifacts.
+        """
+
+        final_grid = vdatum_grid.copy()
+        vdatum_mask = np.isnan(vdatum_grid)
+        if not vdatum_mask.all():
+            nearest_idx = ndimage.distance_transform_edt(vdatum_mask, return_distances=False, return_indices=True)
+            nearest_vdatum_vals = vdatum_grid[tuple(nearest_idx)]
+
+            fes_anomaly_mask = np.abs(global_grid - nearest_vdatum_vals) > max_discontinuity
+
+            global_grid[fes_anomaly_mask] = np.nan
+
+        is_vdatum = ~np.isnan(vdatum_grid)
+        is_ocean = ~np.isnan(global_grid)
+
+        is_inland = ~is_vdatum & ~is_ocean
+        is_offshore = ~is_vdatum & is_ocean
+
+        if is_offshore.any():
+            blended_ocean = GridEngine.smart_blend(vdatum_grid, global_grid, blend_pixels=50)
+            final_grid[is_offshore] = blended_ocean[is_offshore]
+
+        if is_inland.any():
+            decayed_inland = GridEngine.fill_nans(vdatum_grid, decay_pixels=decay_pixels, buffer_pixels=buffer_pixels)
+            final_grid[is_inland] = decayed_inland[is_inland]
+
+        return final_grid
+
+    @staticmethod
+    def fill_nans(data, decay_pixels=100, buffer_pixels=50):
+        """Extrapolates nearest valid value for 'buffer_pixels',
+        then decays to zero over 'decay_pixels'.
+        """
 
         mask = np.isnan(data)
         if not mask.any() or mask.all():
@@ -127,7 +185,10 @@ def fill_nans(data, decay_pixels=100):
             mask, return_distances=True, return_indices=True
         )
         coast_values = data[tuple(indices)]
-        decay_factor = np.clip((decay_pixels - dist) / decay_pixels, 0, 1)
+
+        # Subtract the buffer from the distance. Anything <= 0 is in the buffer zone.
+        effective_dist = np.clip(dist - buffer_pixels, 0, None)
+        decay_factor = np.clip((decay_pixels - effective_dist) / decay_pixels, 0, 1)
 
         out_data = data.copy()
         out_data[mask] = coast_values[mask] * decay_factor[mask]

src/transformez/transform.py

@@ -33,7 +33,7 @@ class VerticalTransform:
 
     def __init__(self, region, nx, ny, epsg_in, epsg_out,
                  geoid_in=None, geoid_out=None, epoch_in=2010.0, epoch_out=2010.0,
-                 cache_dir=None, verbose=True):
+                 decay_pixels=100, cache_dir=None, verbose=True):
 
         self.region = region
         self.nx = nx
@@ -55,6 +55,8 @@ def __init__(self, region, nx, ny, epsg_in, epsg_out,
         self.ref_in = Datums.get_frame_type(self.epsg_in)
         self.ref_out = Datums.get_frame_type(self.epsg_out)
 
+        self.decay_pixels = decay_pixels
+
         # --- HUB SELECTION ---
         # Determine the Native Ellipsoid of Input and Output
         native_in = self._get_native_ellipsoid(self.epsg_in, self.ref_in)
@@ -178,9 +180,9 @@ def _get_grid(self, provider, name):
         if provider == 'seanoe' or provider == 'fes':
             var_name = "lat_elevation" if "lat" in name.lower() else "msl_elevation"
             nc_path = f"netcdf:{files[0]}:{var_name}"
-            return GridEngine.load_and_interpolate([nc_path], self.region, self.nx, self.ny)
+            return GridEngine.load_and_interpolate([nc_path], self.region, self.nx, self.ny, decay_pixels=self.decay_pixels)
 
-        return GridEngine.load_and_interpolate(files, self.region, self.nx, self.ny)
+        return GridEngine.load_and_interpolate(files, self.region, self.nx, self.ny, decay_pixels=self.decay_pixels)
 
     def _get_htdp_shift(self, epsg_from, epsg_to, epoch_from, epoch_to):
         """Calculate Frame Shift via HTDP."""
@@ -232,33 +234,65 @@ def _fetch_geoid_with_fallback(self, target_geoid):
     # =========================================================================
     def _get_vdatum_chain(self, datum_name, geoid_name):
         """Builds shift: Tidal -> [NAD83 Native]."""
-        total_shift = np.zeros((self.ny, self.nx))
+        hydro_shift = np.zeros((self.ny, self.nx))
         desc = []
 
         # Tidal -> LMSL
         if datum_name not in ['msl', '5714', 'lmsl']:
             grid = self._get_grid('vdatum', datum_name)
-            if not np.any(grid):
+            if np.isnan(grid).all() or (grid == 0).all():
                 return None, f"Missing Tidal Grid: {datum_name}"
-            total_shift += grid
+            hydro_shift += grid
             desc.append(f"({datum_name}->LMSL)")
 
         # LMSL -> Ortho (TSS)
         tss = self._get_grid('vdatum', 'tss')
-
-        if not np.any(tss):
+        if np.isnan(tss).all() or (tss == 0).all():
             return None, "Outside VDatum coverage (Missing TSS)"
 
-        total_shift += tss
+        hydro_shift += tss
         desc.append("TSS(LMSL->NAVD88)")
 
-        # Ortho -> NAD83 (Smart Geoid Fallback)
+        # Ortho -> NAD83 (Geoid)
+        # We fetch the geoid, but DO NOT add it to the shift yet!
         actual_geoid = geoid_name if geoid_name else 'g2018'
         geoid_grid, used_geoid = self._fetch_geoid_with_fallback(actual_geoid)
-
-        total_shift += geoid_grid
         desc.append(f"Geoid({used_geoid}->NAD83)")
 
+        # =======================================================
+        # Coastal Blend
+        # =======================================================
+        total_shift = np.zeros((self.ny, self.nx))
+
+        if np.isnan(hydro_shift).any():
+            proxy_name = Datums.get_global_proxy(datum_name)
+            if proxy_name:
+                logger.info(f"Partial VDatum coverage detected. Fetching {proxy_name.upper()} (FES) for offshore blending...")
+                global_shift, d_global = self._get_global_chain(proxy_name, model='fes2014')
+
+                if global_shift is not None and np.any(global_shift):
+                    # We have valid FES data. We must align it to NAD83.
+                    htdp_wgs_to_nad = self._get_htdp_shift(WGS84_EPSG, NAD83_EPSG, self.epoch_in, 2010.0)
+                    fes_full = global_shift + htdp_wgs_to_nad
+
+                    hydro_shift = GridEngine.coastal_aware_composite(
+                        vdatum_grid=hydro_shift,
+                        global_grid=fes_full,
+                        decay_pixels=self.decay_pixels,
+                        buffer_pixels=10,
+                        max_discontinuity=0.5
+                    )
+                    desc.append(f"Blended w/ Global({proxy_name.upper()})")
+                else:
+                    hydro_shift = GridEngine.fill_nans(hydro_shift, decay_pixels=self.decay_pixels, buffer_pixels=10)
+                    desc.append("Inland Hydro Decay")
+            else:
+                hydro_shift = GridEngine.fill_nans(hydro_shift, decay_pixels=self.decay_pixels, buffer_pixels=10)
+                desc.append("Inland Hydro Decay")
+
+        total_shift = hydro_shift + geoid_grid
+        total_shift[np.isnan(total_shift)] = 0.0
+
         return total_shift, " + ".join(desc)
 
     def _get_global_chain(self, datum_name, model="fes2014"):