From 42264410da0eece05167aea5a0eab0b1149d573b Mon Sep 17 00:00:00 2001 From: "Zhengtao.Cui" Date: Tue, 10 Feb 2026 22:12:18 +0000 Subject: [PATCH] added the codes to create SFINCS model form a given polygon. --- coastal/SFINCS/model_creation/README.md | 60 +++++ coastal/SFINCS/model_creation/SCHISMNcGrid.py | 96 ++++++++ .../SFINCS/model_creation/SCHISMnwmReaches.py | 59 +++++ .../model_creation/sfincs_model_setup.py | 223 ++++++++++++++++++ 4 files changed, 438 insertions(+) create mode 100644 coastal/SFINCS/model_creation/README.md create mode 100644 coastal/SFINCS/model_creation/SCHISMNcGrid.py create mode 100644 coastal/SFINCS/model_creation/SCHISMnwmReaches.py create mode 100755 coastal/SFINCS/model_creation/sfincs_model_setup.py diff --git a/coastal/SFINCS/model_creation/README.md b/coastal/SFINCS/model_creation/README.md new file mode 100644 index 0000000..6b6fbfe --- /dev/null +++ b/coastal/SFINCS/model_creation/README.md @@ -0,0 +1,60 @@ +# Overview +This directory contains Python scripts to create SFINCS model files for a given polygon in one of the 4 NWM v3 coastal domains, that is, hawaii, prvi, pacific and atlgulf. + + +# Prerequisite +First, the bathymetry data for the NWM v3 domains are needed. The bathmetry data is specified in the HydroMT datacatlog file in YAML format. Here is an example of a datacatalog file. + +```yaml +>meta: +> version: '1' +> root: ../../hydromt_sfincs/examples/tmpdir +>hawaii: +> crs: 4326 +> data_type: RasterDataset +> driver: raster +> meta: +> category: topography +> source_url: s3://edfs-data/surface/nws-topobathy  +> unit: m+EGM2008 +> version: 2010 +> path: hawaii_30m.tif +``` + +Second, a geojson polygon file that defines the modeling area is needed. +Third, a geojson polygon file that defines the seaward boundary is needed to define the waterlevel boundary. +Third, two SCHISM parameter files, hgrid.nc and nwmReaches.csv, are needed to create the sfincs.src file. + +# Script description +makeOcean_fvcom.py : +sfincs_model_setup.py : The main driver of the Python scripts. Includes the top level logic for creating a SFINCS model files for a given area. +SCHISMNcGrid.py : The class to manage a SCHISM grid, hgrid.nc, file. +SCHISMnwmReaches : The class to manage a SCHISM to NWM domain crosswalk file. + + +# Script Usage +The driver script has a help option (-h) for printing usage information. + +```python +>usage: sfincs_model_setup.py [-h] +> data_catalog model_root_dir polygon_file nwm_domain resolution +> waterlevel_boundary_file hgrid_file nwmreach_file +> +>positional arguments: +> data_catalog the hydromt data catalog filename +> model_root_dir the directory where model files will be saved +> polygon_file the jeojson file defines the simulation domain +> nwm_domain the nwm domain name, one of atlgulf, pacific, prvi and hawaii +> resolution defind grid resolution in meters +> waterlevel_boundary_file +> the jeoson file defines the waterlevel boundary +> hgrid_file SCHISM hgrid.nc file +> nwmreach_file SCHISM nwmReaches.csv file +> +>options: +> -h, --help show this help message and exit +``` +#### Examples #### +```bash +./sfincs_model_setup.py ./data_catalog_atlgulf.yml sfincs_atlgulf /shareds3/zhengtao/domain_polygon_atl3.geojson "atlgulf" 200 /shareds3/zhengtao/boundary_atl2.geojson /efs/ngwpc-coastal/parm/coastal/atlgulf/hgrid.nc /efs/ngwpc-coastal/parm/coastal/atlgulf/nwmReaches.csv +``` diff --git a/coastal/SFINCS/model_creation/SCHISMNcGrid.py b/coastal/SFINCS/model_creation/SCHISMNcGrid.py new file mode 100644 index 0000000..5b95eca --- /dev/null +++ b/coastal/SFINCS/model_creation/SCHISMNcGrid.py @@ -0,0 +1,96 @@ +############################################################################### +# Module name: SCHISNcMGrid # +# # +# Author : Zhengtao Cui (Zhengtao.Cui@rtx.com) # +# # +# Initial version date: 12/12/2025 # +# # +# Last modification date: # +# # +# Description: Manages a SCHISM hgrid.nc file # +# # +############################################################################### +import numpy as np +import xarray as xr +from pyproj import Transformer + +class SCHISMNcGrid: + """ + SCHISM hgrid.nc file + """ + + def __init__(self, hgridncfile ): + self.source = hgridncfile + + #self.elems = dict() + #self.centerCoords = dict() + #self.nodeCoords = dict() + #self.numElementCounts = [] + grid = xr.open_dataset(hgridncfile) + + self._elems = grid['elementConn'].data + self._centerCoords = grid['centerCoords'].data + self._nodeCoords = grid['nodeCoords'].data + self._numElementCounts = grid['numElementConn'] + + grid.close() + + + @property + def source(self): + return self._source + + @source.setter + def source(self, s): + self._source = s + + @property + def elems(self): + return self._elems + + @elems.setter + def elems(self, e): + self._elems = e + + @property + def centerCoords(self): + return self._centerCoords + + @centerCoords.setter + def centerCoords(self, c): + self._centerCoords = c + + @property + def nodeCoords(self): + return self._nodeCoords + + @nodeCoords.setter + def nodeCoords(self, n): + self._nodeCoords = n + + @property + def numElementCounts(self): + return self._numElementCounts + + @numElementCounts.setter + def numElementCounts(self, n): + self._numElementCounts = n + + def getElemLat(self, e ): + return self.centerCoords[e, 1] + + def getElemLon(self, e ): + return self.centerCoords[e, 0] + + def getElemCenterCoordsInCrs(self, e, crs_to ): + # EPSG:4326 is WGS84 (latitude/longitude) + transformer = Transformer.from_crs("EPSG:4326", crs_to, always_xy=True) + + # Perform the transformation on a specific coordinate + # Input order is typically longitude, latitude (x, y) for geographic CRS + # Return order is the same as the input order + + return transformer.transform(self.getElemLon(e), self.getElemLat(e)) + + + diff --git a/coastal/SFINCS/model_creation/SCHISMnwmReaches.py b/coastal/SFINCS/model_creation/SCHISMnwmReaches.py new file mode 100644 index 0000000..d257d8a --- /dev/null +++ b/coastal/SFINCS/model_creation/SCHISMnwmReaches.py @@ -0,0 +1,59 @@ +############################################################################### +# Module name: SCHISMnwmReaches # +# # +# Author : Zhengtao Cui (Zhengtao.Cui@rtx.com) # +# # +# Initial version date: 12/12/2025 # +# # +# Last modification date: # +# # +# Description: Manages a SCHISM corsswalk file, nwmReaches.csv # +# # +############################################################################### +import os +from os import path + +class SCHISMnwmReaches: + """ + SCHISM nwmReaches.csv file + """ + + def __init__(self, reachfile ): + self.source = reachfile + + self._soelem_ids = dict() + self._sielem_ids = dict() + with open(path.join(reachfile)) as f: + nso = int(f.readline()) + for i in range(nso): + line = f.readline() + self._soelem_ids[int(line.split()[0]) -1] = line.split()[1] + next(f) + nsi = int(f.readline()) + for i in range(nsi): + line = f.readline() + self._sielem_ids[int(line.split()[0]) -1 ] = line.split()[1] + + @property + def source(self): + return self._source + + @source.setter + def source(self, s): + self._source = s + + @property + def soelem_ids(self): + return self._soelem_ids + + @soelem_ids.setter + def soelem_ids(self, e): + self._soelem_ids = e + + @property + def sielem_ids(self): + return self._sielem_ids + + @sielem_ids.setter + def sielem_ids(self, c): + self._sielem_ids = c diff --git a/coastal/SFINCS/model_creation/sfincs_model_setup.py b/coastal/SFINCS/model_creation/sfincs_model_setup.py new file mode 100755 index 0000000..7e18d5c --- /dev/null +++ b/coastal/SFINCS/model_creation/sfincs_model_setup.py @@ -0,0 +1,223 @@ +#!/usr/bin/env python + +# Step 0: Imports +import os +import logging +import argparse +import requests +import numpy as np +import pandas as pd +import geopandas as gpd +from shapely.geometry import Point +from pyproj import CRS + +from hydromt_sfincs import SfincsModel +from hydromt_sfincs import utils + +from SCHISMNcGrid import SCHISMNcGrid +from SCHISMnwmReaches import SCHISMnwmReaches + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument('data_catalog', type=str, help='the hydromt data catalog filename') + parser.add_argument('model_root_dir', type=str, help='the directory where model files will be saved') + parser.add_argument('polygon_file', type=str, help='the jeojson file defines the simulation domain') + parser.add_argument('nwm_domain', type=str, help='the nwm domain name, one of atlgulf, pacific, prvi and hawaii') + parser.add_argument('resolution', type=str, help='defind grid resolution in meters', default=200) + parser.add_argument('waterlevel_boundary_file', type=str, help='the jeoson file defines the waterlevel boundary') + parser.add_argument('hgrid_file', type=str, help='SCHISM hgrid.nc file') + parser.add_argument('nwmreach_file', type=str, help='SCHISM nwmReaches.csv file') + + args = parser.parse_args() + + # Step 1: Initialize the SfincsModel class with the data catalog, the "w+" option overwrites an existing model in the root directory if it exists + sf = SfincsModel(data_libs=[args.data_catalog], root=args.model_root_dir, mode="w+") + + # Step 2: Use the pre-defined polygon geojson to define the grid + + #HydorMT doesn't support crs from a string, only support EPSG codes + #lcc = CRS.from_proj4(proj4_str) + #lcc = CRS.from_user_input(esri_pe_string) + + epsg_code = None + if ( args.nwm_domain == 'hawaii' ): + epsg_code = 32604 #UTM + elif ( args.nwm_domain == 'prvi' ): + epsg_code = 3920 #UTM + elif ( args.nwm_domain == 'atlgulf' or args.nwm_domain == 'pacific' ): + epsg_code = 5070 + else: + raise ValueError(f"Unknown nwm_domain name - {args.nwm_domain}") + + sf.setup_grid_from_region( + region={'geom': args.polygon_file}, # we will use the open_boundary from SCHISM, which outlines this particular domain + res=float(args.resolution), # resolution in meters - set to be fairly coarse for initial testing + rotated=False, # do not rotate + crs=epsg_code, # crs/ epsg code + ) + + # Display the automatically generated input file + print(sf.config) + + da = sf.data_catalog.get_rasterdataset(args.nwm_domain, geom=sf.region, buffer=5) + + # show the model grid outline + _ = sf.plot_basemap(fn_out=f"model_grid.png", plot_region=True, bmap="sat", zoomlevel = 12) + + # Step 3: Load elevation dataset(s) and map them to the model grid + + datasets_dep = [ + #{"elevtn":""}, # optional zmin, zmax let's you determine what part of data you want + #{"elevtn":""}, # the dataset order is important + {"elevtn":args.nwm_domain}, + #{} + ] + + dep = sf.setup_dep(datasets_dep=datasets_dep) + + # plot + _ = sf.plot_basemap(fn_out="f{args.nwm_domain}_elev.png", variable='dep',plot_region=True, bmap="sat", zoomlevel=12) + + # Step 4a: Make the mask by setting inactive cells (=0) and active cells (=1) + + sf.setup_mask_active(include_mask=args.polygon_file, reset_mask=True) + + # plot + _ = sf.plot_basemap(fn_out=f"{args.nwm_domain}_mask.png", variable="msk",plot_region=False,bmap="sat",zoomlevel=12) + + # Step 4b (optional): exclude SCHISM land boundary areas (set inactive (=0)) + + # These are messing up the water level boundary part - skipping this step + + # These are the land boundaries from SCHISM - the user may wish to change the inland extents to something different + # exclude_files = [ + # 'hi_boundaries\land_boundary_0.geojson', + # 'hi_boundaries\land_boundary_1.geojson', + # 'hi_boundaries\land_boundary_2.geojson', + # 'hi_boundaries\land_boundary_3.geojson', + # 'hi_boundaries\land_boundary_4.geojson', + # ] + + # gdf_exclude = gdf_exclude = gpd.GeoDataFrame( + # pd.concat([gpd.read_file(f) for f in exclude_files], ignore_index=True), + # crs=gpd.read_file(exclude_files[0]).crs, + # ) + + # sf.setup_mask_active(exclude_mask=gdf_exclude,reset_mask=False) + + # _ = sf.plot_basemap(variable="msk",plot_region=False,bmap="sat",zoomlevel=12) + + # Step 5: Set the mask open water level boundary (=2) + + sf.setup_mask_bounds(btype="waterlevel",include_mask=args.waterlevel_boundary_file, reset_bounds=True) + + # plot + _ = sf.plot_basemap(fn_out=f"{args.nwm_domain}_make_wl_bnd.png", variable="msk",plot_region=False,bmap="sat",zoomlevel=12) + + # Step 6: Add river inflow points + + # we will need to connect to NWM/NGen here, will skip this step for now + # Step 7: Spatially varying roughness data + + # do we have access to a roughness dataset? or a land use dataset that we can convert to roughness? Will skip this step for now and set Manning in step 8. + # Step 8: Make Subgrid Tables + + sf.setup_subgrid( + datasets_dep = datasets_dep, + nr_subgrid_pixels = 20, + manning_land = 0.04, + manning_sea = 0.02, + write_dep_tif = True, + ) + + # Step 9: Add spatially varying infiltration data + + # Will skip this step, but option to include it is available + # Step 10: Set water level boundary points + + # Loops around the water level boundary selecting points every 25,000m + #sf.setup_waterlevel_bnd_from_mask(distance=25000,merge=False) + sf.setup_waterlevel_bnd_from_mask(distance=20000,merge=False) + + _ = sf.plot_basemap(fn_out=f"{args.nwm_domain}_wl_bnd_pts.png", plot_region=False,bmap="sat",zoomlevel=12) + + # Step 11: Add observation points + + # Query lat lon points of NOAA stations within the open water boundary geojson + REGION_GEOJSON = args.polygon_file + BASE_URL = "https://api.tidesandcurrents.noaa.gov/mdapi/prod/webapi/stations.json" + + def fetch_waterlevel_stations(): + """Fetch NOAA CO-OPS stations with water-level data.""" + r = requests.get(BASE_URL, params={"type": "waterlevels"}, timeout=30) + r.raise_for_status() + data = r.json().get("stations", []) + return pd.DataFrame( + [ + { + "station_id": s["id"], + "raw_name": s.get("name", ""), + "lat": float(s["lat"]), + "lon": float(s["lng"]), + } + for s in data + ] + ) + + def get_noaa_waterlevel_stations(region_file): + # 1. Load region polygon + region = gpd.read_file(region_file).to_crs(4326) + #region = gpd.read_file(region_file).to_crs(lcc) + geom = region.geometry.unary_union + + # 2. Fetch water level stations + df = fetch_waterlevel_stations() + + # 3. Format name as "Name (ID)" + df["name"] = df["raw_name"] + " (" + df["station_id"] + ")" + + # 4. Convert to GeoDataFrame + gdf = gpd.GeoDataFrame( + df, + geometry=[Point(lon, lat) for lon, lat in zip(df["lon"], df["lat"])], + crs="EPSG:4326", + #crs=lcc, + ) + + # 5. Spatial filter + return gdf[gdf.within(geom)].copy()[["station_id", "name", "lon", "lat", "geometry"]] + + # ---------------------------- + # Usage + # ---------------------------- + gdf_obs = get_noaa_waterlevel_stations(REGION_GEOJSON) + + sf.setup_observation_points(locations=gdf_obs) + + # Check that the observation points are stored in the sf.geoms dictionary + sf.geoms.keys() + + # Plot the model setup to check it + + _ = sf.plot_basemap(fn_out=f"{args.nwm_domain}_mode_setup.png", bmap="sat",zoomlevel=12) + + # And save it + + sf.write() + + #create the sfincs.src file + + reaches = SCHISMnwmReaches(os.path.join(args.nwmreach_file)) + + schm_grid = SCHISMNcGrid( args.hgrid_file ) + + with open(f"{args.model_root_dir}/sfincs_nwm.src", "w") as file: + for e in reaches.soelem_ids: + x,y = schm_grid.getElemCenterCoordsInCrs(e, f"EPSG:{epsg_code}" ) + file.write(f"{x:.6f} {y:.6f} {reaches.soelem_ids[e]}\n") + +if __name__ == "__main__": + try: + main() + except Exception as e: + logging.error("Failed to get program options.", exc_info=True)