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..4baadf6 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,14 @@ 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 @@ -413,23 +417,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) @@ -442,23 +473,22 @@ def from_stations(region, nx, ny, datum_in, datum_out, shapefiles=None): 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) 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") + 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..a517214 100644 --- a/src/transformez/transform.py +++ b/src/transformez/transform.py @@ -468,41 +468,74 @@ 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", + 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 + # Fetch the global FES MSS baseline to tie the envelope to the Ellipsoid + global_shift, d_global = self._get_global_chain( + "mss", model="fes2014" ) - 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") + + 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 - ) hydro_shift = GridEngine.fill_nans( hydro_shift, decay_pixels=self.decay_pixels, buffer_pixels=10, - land_mask=land_mask, + land_mask=GridEngine.create_land_mask( + self.region, self.nx, self.ny, coast_shapefiles + ), ) desc.append("Inland Hydro Decay") @@ -610,6 +643,21 @@ 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)"