MDT and decay (#56)

* use MDT (with EGM2008) for FES shifts

* decay_pixels = 100 is default

* ruff format

Author: matth-love

src/transformez/api.py

@@ -190,7 +190,7 @@ def transform_raster(
     input_raster: str,
     datum_in: str,
     datum_out: str,
-    decay_pixels: int = 0,
+    decay_pixels: int = 100,
     output_raster: Optional[str] = None,
     cache_dir: Optional[str] = None,
     z_unit_in: Optional[str] = "auto",

src/transformez/cli.py

@@ -174,7 +174,7 @@ def transform_run(
 )
 @click.option("--out", "-o", help="Output filename (default: auto-named).")
 @click.option(
-    "--decay-pixels", type=int, default=0, help="Pixels to decay tidal shifts inland."
+    "--decay-pixels", type=int, default=100, help="Pixels to decay tidal shifts inland."
 )
 @click.option(
     "--use-stations",
@@ -252,7 +252,7 @@ def transform_grid(
 )
 @click.option("--out", "-o", help="Output filename (default: auto-named).")
 @click.option(
-    "--decay-pixels", type=int, default=0, help="Pixels to decay tidal shifts inland."
+    "--decay-pixels", type=int, default=100, help="Pixels to decay tidal shifts inland."
 )
 @click.option(
     "--use-stations",

src/transformez/definitions.py

@@ -402,11 +402,13 @@ def get_unit_factor(cls, unit_str):
             "name": "EGM2008 height",
             "vdatum_id": "egm2008:m:height",
             "default_geoid": "egm2008",
+            "ellipsoid": 4979,
         },
         5773: {
             "name": "EGM96 height",
             "vdatum_id": "egm96:m:height",
             "default_geoid": "egm96",
+            "ellipsoid": 4979,
         },
         # # Ellipsoidal (Hubs) - No Geoid needed
         # 6319: {'name': 'NAD83(2011)', 'vdatum_id': 'nad83_2011:m:height'},

src/transformez/grid_engine.py

@@ -251,20 +251,25 @@ def coastal_aware_composite(
         return final_grid
 
     @staticmethod
-    def fill_nans(data, decay_pixels=0, buffer_pixels=10, land_mask=None):
+    def fill_nans(data, decay_pixels=100, buffer_pixels=10, land_mask=None):
         """Fills NaNs by extrapolating nearest valid coastal values.
         Melted Voronoi ridges ensure C1 continuity deep inland.
         """
 
-        mask = np.isnan(data)
+        out_data = data.copy()
+
+        if land_mask is not None:
+            out_data[~land_mask] = np.nan
+
+        mask = np.isnan(out_data)
         if not mask.any() or mask.all():
-            return data
+            return out_data
 
         dist, indices = ndimage.distance_transform_edt(
             mask, return_distances=True, return_indices=True
         )
 
-        raw_extrapolation = data[tuple(indices)]
+        raw_extrapolation = out_data[tuple(indices)]
         # Blur the "Voronoi Ridges" deep inland
         blurred_extrapolation = ndimage.gaussian_filter(raw_extrapolation, sigma=25)
         # Crossfade! Beach = Raw Data, Inland = Blurred Data
@@ -273,8 +278,6 @@ def fill_nans(data, decay_pixels=0, buffer_pixels=10, land_mask=None):
             blurred_extrapolation * blur_blend
         )
 
-        out_data = data.copy()
-
         if decay_pixels and decay_pixels > 0:
             # --- Inland Decay ---
             effective_dist = np.clip(dist - buffer_pixels, 0, None)

src/transformez/transform.py

@@ -594,12 +594,13 @@ def _get_vdatum_chain(self, datum_name, geoid_name):
     def _get_global_chain(self, datum_name, model="fes2014"):
         """Builds shift: Global Tidal -> WGS84 Native."""
 
-        total_shift = np.zeros((self.ny, self.nx))
+        tidal_shift = np.zeros((self.ny, self.nx))
+        mss_shift = np.zeros((self.ny, self.nx))
         desc = []
 
         model_def = Datums.MODELS.get(model)
         if not model_def:
-            return total_shift, "Error"
+            return tidal_shift, "Error"
 
         provider = model_def["provider"]
 
@@ -608,26 +609,18 @@ def _get_global_chain(self, datum_name, model="fes2014"):
             grid_name = model_def["grids"].get("lat")
             if grid_name:
                 grid = self._get_grid(provider, grid_name)
-                # Sign Correction (LAT should be negative relative to MSS)
                 if np.nanmean(grid) > 0:
                     grid *= -1.0
-
-                total_shift += grid
+                tidal_shift += grid
                 desc.append("LAT->MSS")
 
         elif datum_name == "hat":
             lat_name = model_def["grids"].get("lat")
             if lat_name:
-                logger.info("HAT requested. Computing FES deep-water symmetry...")
                 grid = self._get_grid(provider, lat_name)
-
                 if np.nanmean(grid) > 0:
                     grid *= -1.0
-
-                # Symmetry Math: HAT offset is the exact inverse of LAT offset
-                hat_grid = grid * -1.0
-
-                total_shift += hat_grid
+                tidal_shift += grid * -1.0
                 desc.append("HAT->MSS(Symmetry)")
 
         # MSS -> WGS84
@@ -640,23 +633,36 @@ def _get_global_chain(self, datum_name, model="fes2014"):
             else:
                 desc.append("MSS->Ellipsoid")
 
-            total_shift += mss_grid
-
-        if np.isnan(total_shift).any():
-            coast_shapefiles = self._fetch_coastline_shapefiles()
-            land_mask = None
-            if coast_shapefiles:
-                land_mask = GridEngine.create_land_mask(
-                    self.region, self.nx, self.ny, coast_shapefiles
-                )
+            # --- MDT ---
+            # Isolate the Mean Dynamic Topography (MDT) by subtracting EGM2008
+            egm_grid, _ = self._fetch_geoid_with_fallback("egm2008")
+            mdt_grid = mss_grid - egm_grid
+            desc.append("MDT Decay -> EGM2008")
+
+        coast_shapefiles = self._fetch_coastline_shapefiles()
+        land_mask = None
+        if coast_shapefiles:
+            land_mask = GridEngine.create_land_mask(
+                self.region, self.nx, self.ny, coast_shapefiles
+            )
 
-            total_shift = GridEngine.fill_nans(
-                total_shift,
+        if mss_name:
+            mdt_decayed = GridEngine.fill_nans(
+                mdt_grid,
                 decay_pixels=self.decay_pixels,
                 buffer_pixels=10,
                 land_mask=land_mask,
             )
+            mss_shift = mdt_decayed + egm_grid
 
+        tidal_shift = GridEngine.fill_nans(
+            tidal_shift,
+            decay_pixels=self.decay_pixels,
+            buffer_pixels=10,
+            land_mask=land_mask,
+        )
+
+        total_shift = tidal_shift + mss_shift
         if not desc:
             return total_shift, "Global Chain (Empty)"
 
@@ -764,7 +770,8 @@ def _step_from_hub(self, epsg, ref_type, geoid=None, epoch=None):
                 s, d = self._get_global_chain(datum_name)
                 if s is not None:
                     total_out -= s
-                    desc_parts.append(f"Native -> Global({datum_name})")
+                    # desc_parts.append(f"Native -> Global({datum_name})")
+                    desc_parts.append(d)
 
         elif ref_type == "cdn":
             target_geoid = geoid if geoid else "g2018"