A Python script to create SFINCS model files using HydorMT-SFINCS

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
+```

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))
+
+
+

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