Skip to content

carve vdatum pixels out of the land_mask #57

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
matth-love merged 1 commit into from
Apr 15, 2026
Merged

carve vdatum pixels out of the land_mask #57

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 10 additions & 4 deletions src/transformez/grid_engine.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -211,7 +211,8 @@ def coastal_aware_composite(
region,
nx,
ny,
shapefiles=None,
shapefiles=None, # Deprecated
land_mask=None,
decay_pixels=100,
buffer_pixels=10,
blend_pixels=50,
Expand All @@ -221,11 +222,16 @@ def coastal_aware_composite(
"""

final_grid = vdatum_grid.copy()
land_mask = None
if shapefiles:

# Fallback just in case shapefiles are passed directly
if shapefiles and land_mask is None:
land_mask = GridEngine.create_land_mask(region, nx, ny, shapefiles)
if land_mask is not None:
global_grid[~land_mask] = np.nan
# Carve out the VDatum rivers!
land_mask[~np.isnan(vdatum_grid)] = False

if land_mask is not None:
global_grid[~land_mask] = np.nan

is_vdatum = ~np.isnan(vdatum_grid)
is_ocean = ~np.isnan(global_grid)
Expand Down
44 changes: 15 additions & 29 deletions src/transformez/transform.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -468,11 +468,17 @@ 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 <---
land_mask = None
if coast_shapefiles:
land_mask = GridEngine.create_land_mask(
self.region, self.nx, self.ny, coast_shapefiles
)
if land_mask is not None:
valid_vdatum = ~np.isnan(hydro_shift)
land_mask[valid_vdatum] = True

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,
Expand All @@ -483,7 +489,6 @@ def _get_vdatum_chain(self, datum_name, geoid_name):
)

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"
)
Expand All @@ -494,16 +499,11 @@ def _get_vdatum_chain(self, datum_name, geoid_name):
)
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,
Expand All @@ -512,12 +512,7 @@ def _get_vdatum_chain(self, datum_name, geoid_name):
)
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
)
logger.warning(" [Override] FES baseline missing.")
hydro_shift = GridEngine.fill_nans(
hydro_shift,
decay_pixels=self.decay_pixels,
Expand All @@ -526,16 +521,12 @@ def _get_vdatum_chain(self, datum_name, geoid_name):
)
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.")
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
),
land_mask=land_mask,
)
desc.append("Inland Hydro Decay")

Expand All @@ -548,7 +539,6 @@ def _get_vdatum_chain(self, datum_name, geoid_name):
)

if global_shift is not None and np.any(global_shift):
# FES is Proxy -> WGS84. We convert to NAD83.
htdp_wgs_to_nad = self._get_htdp_shift(
WGS84_EPSG, NAD83_EPSG, self.epoch_in, 2010.0
)
Expand All @@ -560,7 +550,7 @@ def _get_vdatum_chain(self, datum_name, geoid_name):
region=self.region,
nx=self.nx,
ny=self.ny,
shapefiles=coast_shapefiles,
land_mask=land_mask,
decay_pixels=self.decay_pixels,
buffer_pixels=10,
)
Expand All @@ -570,19 +560,15 @@ def _get_vdatum_chain(self, datum_name, geoid_name):
hydro_shift,
decay_pixels=self.decay_pixels,
buffer_pixels=10,
land_mask=GridEngine.create_land_mask(
self.region, self.nx, self.ny, coast_shapefiles
),
land_mask=land_mask,
)
desc.append("Inland Hydro Decay")
else:
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
),
land_mask=land_mask,
)
desc.append("Inland Hydro Decay")

Expand Down
Loading