From 264bb8664c639d86c417af94f930ea4ce721800b Mon Sep 17 00:00:00 2001 From: Matthew Love Date: Tue, 14 Apr 2026 14:47:46 -0700 Subject: [PATCH 1/2] overhaul; inland-decay is optional, extrapolate is the default; use tide station navd88 if exists otherwise msl->fes --- src/transformez/api.py | 2 +- src/transformez/cli.py | 4 +- src/transformez/grid_engine.py | 98 +++++++++++++++++++++------------- src/transformez/transform.py | 75 ++++++++++++++++---------- 4 files changed, 111 insertions(+), 68 deletions(-) diff --git a/src/transformez/api.py b/src/transformez/api.py index 47063f2..a48e07e 100644 --- a/src/transformez/api.py +++ b/src/transformez/api.py @@ -190,7 +190,7 @@ def transform_raster( input_raster: str, datum_in: str, datum_out: str, - decay_pixels: int = 100, + decay_pixels: int = 0, output_raster: Optional[str] = None, cache_dir: Optional[str] = None, z_unit_in: Optional[str] = "auto", diff --git a/src/transformez/cli.py b/src/transformez/cli.py index a45be57..c6172ed 100644 --- a/src/transformez/cli.py +++ b/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=100, help="Pixels to decay tidal shifts inland." + "--decay-pixels", type=int, default=0, 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=100, help="Pixels to decay tidal shifts inland." + "--decay-pixels", type=int, default=0, help="Pixels to decay tidal shifts inland." ) @click.option( "--use-stations", diff --git a/src/transformez/grid_engine.py b/src/transformez/grid_engine.py index b81b7f8..428ead4 100644 --- a/src/transformez/grid_engine.py +++ b/src/transformez/grid_engine.py @@ -261,11 +261,10 @@ def coastal_aware_composite( return final_grid @staticmethod - def fill_nans(data, decay_pixels=100, buffer_pixels=50, land_mask=None): - """Extrapolates nearest valid value for 'buffer_pixels', - then decays to zero over 'decay_pixels'. - - If land_mask is provided, prevents decay from covering land areas. + def fill_nans(data, decay_pixels=0, buffer_pixels=10, land_mask=None): + """Fills NaNs by extrapolating nearest valid coastal values. + If decay_pixels > 0, it will smoothly ramp the values to 0.0 inland. + If decay_pixels == 0, it extrapolates the coastal values infinitely inland. """ mask = np.isnan(data) @@ -276,18 +275,20 @@ def fill_nans(data, decay_pixels=100, buffer_pixels=50, land_mask=None): mask, return_distances=True, return_indices=True ) coast_values = data[tuple(indices)] - - # 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] - if land_mask is not None: - # We only mask out areas that were originally nan. - # If VDatum contained valid data over "land", we keep it! - out_data[mask & land_mask] = 0.0 + if decay_pixels and decay_pixels > 0: + # --- Inland Decay --- + effective_dist = np.clip(dist - buffer_pixels, 0, None) + decay_factor = np.clip((decay_pixels - effective_dist) / decay_pixels, 0, 1) + out_data[mask] = coast_values[mask] * decay_factor[mask] + + if land_mask is not None: + # Force areas deep inland (and covered by the land mask) to exactly 0.0 + out_data[mask & ~land_mask] = 0.0 + else: + # --- Infinite Nearest-Neighbor Extrapolation (Default) --- + out_data[mask] = coast_values[mask] return out_data @@ -378,11 +379,12 @@ def write(filename, data, region): class GridGen: @staticmethod - def from_stations(region, nx, ny, datum_in, datum_out, shapefiles=None): - """Generates a tidal shift grid (datum_in -> datum_out) - using live NOAA CO-OPS tide station data and RBF interpolation. + def from_stations(region, nx, ny, datum_in, datum_out, shapefiles=None, baseline_grid=None): + """Dynamically generates a tidal shift grid using live tide stations. + If a station lacks the target datum, it falls back to MSL and uses the + baseline_grid (FES) to bridge the gap to the geodetic frame. """ - + import os import json from fetchez.modules.tides import Tides import fetchez @@ -401,7 +403,7 @@ def from_stations(region, nx, ny, datum_in, datum_out, shapefiles=None): logger.error(f"GeoJSON file not found: {geojson_path}") return None - with open(geojson_path, "r") as f: + with open(geojson_path, 'r') as f: data = json.load(f) features = data.get("features", []) @@ -413,23 +415,50 @@ def from_stations(region, nx, ny, datum_in, datum_out, shapefiles=None): d_in = datum_in.lower() d_out = datum_out.lower() + # Calculate pixel resolution for baseline grid sampling + res_x = (region.xmax - region.xmin) / nx + res_y = (region.ymax - region.ymin) / ny + for feat in features: props = feat.get("properties", {}) geom = feat.get("geometry", {}) - if d_in in props and d_out in props: - val_in = props.get(d_in) - val_out = props.get(d_out) + val_in = props.get(d_in) + if val_in is None or val_in < -90000: + continue - if val_in is None or val_out is None: - continue + val_out = props.get(d_out) + shift = None + units = props.get("units", "meters").lower() + # --- Perfect Data (Station has NAVD88) --- + if val_out is not None and val_out > -90000: shift = val_in - val_out - - units = props.get("units", "meters").lower() if units == "feet": shift *= 0.3048 + # --- Floating Station (Lacks NAVD88, but has MSL) --- + elif baseline_grid is not None and "msl" in props: + val_msl = props.get("msl") + if val_msl is not None and val_msl > -90000: + lon = geom["coordinates"][0] + lat = geom["coordinates"][1] + + x_idx = int((lon - region.xmin) / res_x) + y_idx = int((region.ymax - lat) / res_y) + + if 0 <= x_idx < nx and 0 <= y_idx < ny: + fes_offset = baseline_grid[y_idx, x_idx] + if not np.isnan(fes_offset): + # Get the local tidal envelope to MSL + shift_to_msl = val_in - val_msl + if units == "feet": + shift_to_msl *= 0.3048 + + # Add the FES baseline offset to mathematically tie it to NAVD88 + shift = shift_to_msl + fes_offset + + if shift is not None: x.append(geom["coordinates"][0]) y.append(geom["coordinates"][1]) z.append(shift) @@ -439,26 +468,21 @@ def from_stations(region, nx, ny, datum_in, datum_out, shapefiles=None): return None if len(z) < 3: - logger.warning( - f"Only {len(z)} station(s) found. Applying a constant average offset instead of RBF." - ) - # Average the offsets (works for 1 or 2 stations) + logger.warning(f"Only {len(z)} station(s) found. Applying a constant average offset instead of RBF.") constant_shift = sum(z) / len(z) - # Create a flat sheet rbf_grid = np.full((ny, nx), constant_shift, dtype=np.float32) else: - logger.info( - f"Interpolating surface using {len(z)} coastal tide stations..." - ) - rbf = Rbf(x, y, z, function="thin_plate") + logger.info(f"Interpolating surface using {len(z)} coastal tide stations...") + rbf = Rbf(x, y, z, function='linear') xi = np.linspace(region.xmin, region.xmax, nx) yi = np.linspace(region.ymax, region.ymin, ny) XI, YI = np.meshgrid(xi, yi) rbf_grid = rbf(XI, YI).astype(np.float32) + rbf_grid = np.clip(rbf_grid, min(z), max(z)) if shapefiles: - logger.info("Applying vector coastline mask to RBF surface...") + logger.info("Applying vector coastline mask to station surface...") land_mask = GridEngine.create_land_mask(region, nx, ny, shapefiles) if land_mask is not None: rbf_grid[~land_mask] = np.nan diff --git a/src/transformez/transform.py b/src/transformez/transform.py index 1824d08..de8fea7 100644 --- a/src/transformez/transform.py +++ b/src/transformez/transform.py @@ -468,41 +468,47 @@ def _get_vdatum_chain(self, datum_name, geoid_name): coast_shapefiles = self._fetch_coastline_shapefiles() proxy_name = Datums.get_global_proxy(datum_name) + # ---> CHECK FLAG FIRST <--- if self.use_stations: logger.info(" [Override] Forcing Tide Station RBF interpolation...") + + # Ask the RBF to solve the local tidal envelope (to MSL instead of NAVD88) rbf_grid = GridGen.from_stations( - self.region, - self.nx, - self.ny, - datum_in=datum_name, - datum_out="navd88", - shapefiles=coast_shapefiles, + self.region, self.nx, self.ny, datum_in=datum_name, datum_out="msl", shapefiles=coast_shapefiles ) + if rbf_grid is not None: - vdatum_empty = np.isnan(hydro_shift) - hydro_shift[vdatum_empty] = rbf_grid[vdatum_empty] - land_mask = GridEngine.create_land_mask( - self.region, self.nx, self.ny, coast_shapefiles - ) - hydro_shift = GridEngine.fill_nans( - hydro_shift, - decay_pixels=self.decay_pixels, - buffer_pixels=10, - land_mask=land_mask, - ) - desc.append("Station RBF + Inland Decay") + # Fetch the global FES MSS baseline to tie the envelope to the Ellipsoid + global_shift, d_global = self._get_global_chain("mss", model="fes2014") + + if global_shift is not None and np.any(global_shift): + htdp_wgs_to_nad = self._get_htdp_shift(WGS84_EPSG, NAD83_EPSG, self.epoch_in, 2010.0) + fes_nad83 = global_shift + htdp_wgs_to_nad + fes_navd88 = fes_nad83 - geoid_grid + + # Stack the local tide on top of the global surface! + combined_shift = rbf_grid + fes_navd88 + + vdatum_empty = np.isnan(hydro_shift) + hydro_shift[vdatum_empty] = combined_shift[vdatum_empty] + + land_mask = GridEngine.create_land_mask(self.region, self.nx, self.ny, coast_shapefiles) + hydro_shift = GridEngine.fill_nans( + hydro_shift, + decay_pixels=self.decay_pixels, + buffer_pixels=10, + land_mask=land_mask + ) + desc.append("Station RBF (Tidal) + FES (MSS) + Inland Decay") + else: + logger.warning(" [Override] FES baseline missing. Cannot tie RBF to Ellipsoid.") + land_mask = GridEngine.create_land_mask(self.region, self.nx, self.ny, coast_shapefiles) + hydro_shift = GridEngine.fill_nans(hydro_shift, decay_pixels=self.decay_pixels, buffer_pixels=10, land_mask=land_mask) + desc.append("Inland Hydro Decay") else: - logger.warning( - " [Override] Tide Station RBF failed. Falling back to inland decay." - ) - land_mask = GridEngine.create_land_mask( - self.region, self.nx, self.ny, coast_shapefiles - ) + logger.warning(" [Override] Tide Station RBF failed. Falling back to inland decay.") hydro_shift = GridEngine.fill_nans( - hydro_shift, - decay_pixels=self.decay_pixels, - buffer_pixels=10, - land_mask=land_mask, + hydro_shift, decay_pixels=self.decay_pixels, buffer_pixels=10, land_mask=GridEngine.create_land_mask(self.region, self.nx, self.ny, coast_shapefiles) ) desc.append("Inland Hydro Decay") @@ -610,6 +616,19 @@ def _get_global_chain(self, datum_name, model="fes2014"): 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) + + total_shift = GridEngine.fill_nans( + total_shift, + decay_pixels=self.decay_pixels, + buffer_pixels=10, + land_mask=land_mask + ) + if not desc: return total_shift, "Global Chain (Empty)" From 4f36e7ce0f00e652773329d6e5b60d41a9850d3b Mon Sep 17 00:00:00 2001 From: Matthew Love Date: Tue, 14 Apr 2026 14:50:49 -0700 Subject: [PATCH 2/2] ruff format --- src/transformez/grid_engine.py | 16 ++++-- src/transformez/transform.py | 93 ++++++++++++++++++++++------------ 2 files changed, 72 insertions(+), 37 deletions(-) diff --git a/src/transformez/grid_engine.py b/src/transformez/grid_engine.py index 428ead4..4baadf6 100644 --- a/src/transformez/grid_engine.py +++ b/src/transformez/grid_engine.py @@ -379,7 +379,9 @@ def write(filename, data, region): class GridGen: @staticmethod - def from_stations(region, nx, ny, datum_in, datum_out, shapefiles=None, baseline_grid=None): + def from_stations( + region, nx, ny, datum_in, datum_out, shapefiles=None, baseline_grid=None + ): """Dynamically generates a tidal shift grid using live tide stations. If a station lacks the target datum, it falls back to MSL and uses the baseline_grid (FES) to bridge the gap to the geodetic frame. @@ -403,7 +405,7 @@ def from_stations(region, nx, ny, datum_in, datum_out, shapefiles=None, baseline logger.error(f"GeoJSON file not found: {geojson_path}") return None - with open(geojson_path, 'r') as f: + with open(geojson_path, "r") as f: data = json.load(f) features = data.get("features", []) @@ -468,13 +470,17 @@ def from_stations(region, nx, ny, datum_in, datum_out, shapefiles=None, baseline return None if len(z) < 3: - logger.warning(f"Only {len(z)} station(s) found. Applying a constant average offset instead of RBF.") + logger.warning( + f"Only {len(z)} station(s) found. Applying a constant average offset instead of RBF." + ) constant_shift = sum(z) / len(z) rbf_grid = np.full((ny, nx), constant_shift, dtype=np.float32) else: - logger.info(f"Interpolating surface using {len(z)} coastal tide stations...") - rbf = Rbf(x, y, z, function='linear') + logger.info( + f"Interpolating surface using {len(z)} coastal tide stations..." + ) + rbf = Rbf(x, y, z, function="linear") xi = np.linspace(region.xmin, region.xmax, nx) yi = np.linspace(region.ymax, region.ymin, ny) XI, YI = np.meshgrid(xi, yi) diff --git a/src/transformez/transform.py b/src/transformez/transform.py index de8fea7..a517214 100644 --- a/src/transformez/transform.py +++ b/src/transformez/transform.py @@ -474,41 +474,68 @@ def _get_vdatum_chain(self, datum_name, geoid_name): # Ask the RBF to solve the local tidal envelope (to MSL instead of NAVD88) rbf_grid = GridGen.from_stations( - self.region, self.nx, self.ny, datum_in=datum_name, datum_out="msl", shapefiles=coast_shapefiles + self.region, + self.nx, + self.ny, + datum_in=datum_name, + datum_out="msl", + shapefiles=coast_shapefiles, ) if rbf_grid is not None: - # Fetch the global FES MSS baseline to tie the envelope to the Ellipsoid - global_shift, d_global = self._get_global_chain("mss", model="fes2014") - - if global_shift is not None and np.any(global_shift): - htdp_wgs_to_nad = self._get_htdp_shift(WGS84_EPSG, NAD83_EPSG, self.epoch_in, 2010.0) - fes_nad83 = global_shift + htdp_wgs_to_nad - fes_navd88 = fes_nad83 - geoid_grid - - # Stack the local tide on top of the global surface! - combined_shift = rbf_grid + fes_navd88 - - vdatum_empty = np.isnan(hydro_shift) - hydro_shift[vdatum_empty] = combined_shift[vdatum_empty] - - land_mask = GridEngine.create_land_mask(self.region, self.nx, self.ny, coast_shapefiles) - hydro_shift = GridEngine.fill_nans( - hydro_shift, - decay_pixels=self.decay_pixels, - buffer_pixels=10, - land_mask=land_mask - ) - desc.append("Station RBF (Tidal) + FES (MSS) + Inland Decay") - else: - logger.warning(" [Override] FES baseline missing. Cannot tie RBF to Ellipsoid.") - land_mask = GridEngine.create_land_mask(self.region, self.nx, self.ny, coast_shapefiles) - hydro_shift = GridEngine.fill_nans(hydro_shift, decay_pixels=self.decay_pixels, buffer_pixels=10, land_mask=land_mask) - desc.append("Inland Hydro Decay") + # Fetch the global FES MSS baseline to tie the envelope to the Ellipsoid + global_shift, d_global = self._get_global_chain( + "mss", model="fes2014" + ) + + if global_shift is not None and np.any(global_shift): + htdp_wgs_to_nad = self._get_htdp_shift( + WGS84_EPSG, NAD83_EPSG, self.epoch_in, 2010.0 + ) + fes_nad83 = global_shift + htdp_wgs_to_nad + fes_navd88 = fes_nad83 - geoid_grid + + # Stack the local tide on top of the global surface! + combined_shift = rbf_grid + fes_navd88 + + vdatum_empty = np.isnan(hydro_shift) + hydro_shift[vdatum_empty] = combined_shift[vdatum_empty] + + land_mask = GridEngine.create_land_mask( + self.region, self.nx, self.ny, coast_shapefiles + ) + hydro_shift = GridEngine.fill_nans( + hydro_shift, + decay_pixels=self.decay_pixels, + buffer_pixels=10, + land_mask=land_mask, + ) + desc.append("Station RBF (Tidal) + FES (MSS) + Inland Decay") + else: + logger.warning( + " [Override] FES baseline missing. Cannot tie RBF to Ellipsoid." + ) + land_mask = GridEngine.create_land_mask( + self.region, self.nx, self.ny, coast_shapefiles + ) + hydro_shift = GridEngine.fill_nans( + hydro_shift, + decay_pixels=self.decay_pixels, + buffer_pixels=10, + land_mask=land_mask, + ) + desc.append("Inland Hydro Decay") else: - logger.warning(" [Override] Tide Station RBF failed. Falling back to inland decay.") + logger.warning( + " [Override] Tide Station RBF failed. Falling back to inland decay." + ) hydro_shift = GridEngine.fill_nans( - hydro_shift, decay_pixels=self.decay_pixels, buffer_pixels=10, land_mask=GridEngine.create_land_mask(self.region, self.nx, self.ny, coast_shapefiles) + hydro_shift, + decay_pixels=self.decay_pixels, + buffer_pixels=10, + land_mask=GridEngine.create_land_mask( + self.region, self.nx, self.ny, coast_shapefiles + ), ) desc.append("Inland Hydro Decay") @@ -620,13 +647,15 @@ def _get_global_chain(self, datum_name, model="fes2014"): 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) + land_mask = GridEngine.create_land_mask( + self.region, self.nx, self.ny, coast_shapefiles + ) total_shift = GridEngine.fill_nans( total_shift, decay_pixels=self.decay_pixels, buffer_pixels=10, - land_mask=land_mask + land_mask=land_mask, ) if not desc: