Read surface topography from Elmer-Xios output

CHANGELOG.md

@@ -6,8 +6,9 @@ SPDX-License-Identifier: CC-BY-4.0
 
 # develop
 
+* Support new input mesh format. EBFM now accepts Elmer/Ice mesh file for xy-coordinates and separate unstructured NetCDF elevation file obtained from XIOS. To use this feature please provide `--elmer-mesh` together with the new option `--netcdf-mesh-unstructured`. See https://github.com/EBFMorg/EBFM/pull/12.
 * Use `setuptools_scm` as backend for `--version` information. See https://github.com/EBFMorg/EBFM/pull/46.
-* Remove `pathlib` from requirements, because this can lead to a bug. https://github.com/EBFMorg/EBFM/pull/48
+* Remove `pathlib` from requirements, because this can lead to a bug. https://github.com/EBFMorg/EBFM/pull/48.
 
 # v0.1.0
 

src/coupling.py

@@ -360,6 +360,7 @@ def exchange_elmer_ice(self, put_data: Dict[str, np.array]) -> Dict[str, np.arra
         received_data = {}
         for field_name, field in self.target_fields.items():
             if field_name in elmer_fields:
+                logger.debug(f"EBFM: {field.component_name}, {field.grid_name}, {field.name}")
                 role = self.interface.get_field_role(field.component_name, field.grid_name, field.name)
                 assert (
                     role == yac.ExchangeType.TARGET

src/ebfm/INIT.py

@@ -11,7 +11,7 @@
 from netCDF4 import Dataset, num2date
 
 from pathlib import Path
-from reader import read_elmer_mesh, read_dem
+from reader import read_elmer_mesh, read_dem, read_dem_xios
 
 from elmer.mesh import Mesh
 from ebfm.config import GridConfig
@@ -153,8 +153,14 @@ def init_constants():
 
 def init_grid(grid, io, config: GridConfig):
     grid["is_partitioned"] = config.is_partitioned
-
-    if config.grid_type is GridInputType.CUSTOM:  # Read grid from Elmer, elevations from BedMachine
+    grid["is_unstructured"] = config.is_unstructured
+
+    # Read grid from Elmer, elevations from BedMachine
+    if config.dem_file:
+        grid_input_type_supporting_dem = [GridInputType.CUSTOM, GridInputType.ELMERXIOS]
+        assert (
+            config.grid_type in grid_input_type_supporting_dem
+        ), f"DEM file can only be specified for {grid_input_type_supporting_dem}."
         if config.is_partitioned:
             mesh: Mesh = read_elmer_mesh(
                 mesh_root=config.mesh_file,
@@ -165,15 +171,18 @@ def init_grid(grid, io, config: GridConfig):
             mesh: Mesh = read_elmer_mesh(mesh_root=config.mesh_file)
 
         grid["x"], grid["y"] = mesh.x_vertices, mesh.y_vertices
-        grid["z"] = read_dem(config.dem_file, grid["x"], grid["y"])
+        if config.grid_type is GridInputType.CUSTOM:
+            grid["z"] = read_dem(config.dem_file, grid["x"], grid["y"])
+            grid["lat"] = np.zeros_like(grid["x"]) + 75  # test values!
+            grid["lon"] = np.zeros_like(grid["x"]) + 320  # test values!
+        if config.grid_type is GridInputType.ELMERXIOS:
+            grid = read_dem_xios(config.dem_file, grid)
+        grid["mask"] = np.ones_like(grid["x"])  # treats every grid cell as glacier
+        grid["gpsum"] = np.sum(grid["mask"] == 1)  # number of modelled grid cells
         grid["slope_x"] = np.zeros_like(grid["x"])  # test values!
         grid["slope_y"] = np.zeros_like(grid["x"])  # test values!
-        grid["lat"] = np.zeros_like(grid["x"]) + 75  # test values!
-        grid["lon"] = np.zeros_like(grid["x"]) + 320  # test values!
         grid["slope_beta"] = np.zeros_like(grid["x"])  # test values!
         grid["slope_gamma"] = np.zeros_like(grid["x"])  # test values!
-        grid["mask"] = np.ones_like(grid["x"])  # treats every grid cell as glacier
-        grid["gpsum"] = np.sum(grid["mask"] == 1)  # number of modelled grid cells
         grid["mesh"] = mesh
         grid["has_shading"] = False  # TODO: see https://github.com/EBFMorg/EBFM/issues/11
         # TODO later add slope
@@ -300,12 +309,33 @@ def init_grid(grid, io, config: GridConfig):
         grid["slope_gamma"][(grid["slope_x"] > 0) & (grid["slope_y"] == 0)] = np.pi / 2
         grid["slope_gamma"][(grid["slope_x"] < 0) & (grid["slope_y"] == 0)] = -np.pi / 2
         grid["slope_gamma"] = -grid["slope_gamma"]
+    elif grid["input_type"] is GridInputType.ELMER_XIOS:
+        grid = read_elmer_xios_grid(grid=grid, gridfile=config.dem_file)
+        grid["gpsum"] = grid["z"].shape[0]
+        grid["mask"] = (grid["h"] > 1.0) * 1.0
+        grid["x"] = np.zeros_like(grid["z"])
+        grid["slope_beta"] = np.zeros_like(grid["x"])  # test values!
+        grid["slope_gamma"] = np.zeros_like(grid["x"])  # test values!
     else:
         raise ValueError(f"Unsupported grid input type {config.grid_type} specified in configuration.")
 
     return grid
 
 
+def read_elmer_xios_grid(grid, gridfile: Path):
+    import netCDF4 as nc
+
+    print(gridfile)
+    with nc.Dataset(gridfile) as file:
+        grid["lat"] = file["mesh2D_node_x"][:]
+        grid["lon"] = file["mesh2D_node_y"][:]
+        grid["x"] = file["x"][:]
+        grid["y"] = file["y"][:]
+        grid["z"] = np.squeeze(file["zs"][:].data)
+        grid["h"] = np.squeeze(file["h"][:].data)
+    return grid
+
+
 def read_MATLAB_grid(gridfile: Path):
     """
     Provides grid information by reading from a .mat file or allowing user input.

src/ebfm/LOOP_write_to_file.py

@@ -141,11 +141,16 @@ def save_netCDF_file():
             # Create NetCDF file
             nc_filepath = os.path.join(io["outdir"], "model_output.nc")
             io["nc_file"] = Dataset(nc_filepath, "w", format="NETCDF4")
-
-            # Define dimensions
             io["nc_file"].createDimension("time", None)  # Unlimited time dimension
-            io["nc_file"].createDimension("y", grid["x_2D"].shape[0])  # 2D grid rows
-            io["nc_file"].createDimension("x", grid["x_2D"].shape[1])  # 2D grid columns
+
+            if grid["is_unstructured"]:
+                # Define dimensions
+                io["nc_file"].createDimension("y", grid["lat"].shape[0])  # 2D grid rows
+            else:  # structured grid
+                # Define dimensions
+                io["nc_file"].createDimension("y", grid["x_2D"].shape[0])  # 2D grid rows
+                io["nc_file"].createDimension("x", grid["x_2D"].shape[1])  # 2D grid columns
+
             io["nc_file"].createDimension("nl", grid["nl"])  # Vertical layers for `sub` variables
 
             # Define standard output variables
@@ -154,35 +159,39 @@ def save_netCDF_file():
                 var_units = entry[1]
                 var_desc = entry[3]
 
-                # Check if variable is a `sub` variable
-                if varname.startswith("sub"):
-                    # Define variable as 4D: (time, y, x, nl)
-                    nc_var = io["nc_file"].createVariable(
-                        varname,
-                        np.float32,
-                        ("time", "y", "x", "nl"),
-                        zlib=True,
-                        complevel=4,
-                        fill_value=-9999.0,  # Fill missing values
-                        chunksizes=(
-                            1,
-                            grid["x_2D"].shape[0],
-                            grid["x_2D"].shape[1],
-                            grid["nl"],
-                        ),
-                    )
+                dimensions: tuple
+                chunksizes: tuple
+
+                if grid["is_unstructured"]:
+                    # Check if variable is a `sub` variable
+                    if varname.startswith("sub"):
+                        # Define variable as 3D: (time, y, nl)
+                        dimensions = ("time", "y", "nl")
+                        chunksizes = (1, grid["lat"].shape[0], grid["nl"])
+                    else:
+                        # Define variable as 2D: (time, y)
+                        dimensions = ("time", "y")
+                        chunksizes = (1, grid["lat"].shape[0])
                 else:
-                    # Define variable as 3D: (time, y, x)
-                    nc_var = io["nc_file"].createVariable(
-                        varname,
-                        np.float32,
-                        ("time", "y", "x"),
-                        zlib=True,
-                        complevel=4,
-                        fill_value=-9999.0,
-                        chunksizes=(1, grid["x_2D"].shape[0], grid["x_2D"].shape[1]),
-                    )
-
+                    # Check if variable is a `sub` variable
+                    if varname.startswith("sub"):
+                        # Define variable as 4D: (time, y, x, nl)
+                        dimensions = ("time", "y", "x", "nl")
+                        chunksizes = (1, grid["x_2D"].shape[0], grid["x_2D"].shape[1], grid["nl"])
+                    else:
+                        # Define variable as 3D: (time, y, x)
+                        dimensions = ("time", "y", "x")
+                        chunksizes = (1, grid["x_2D"].shape[0], grid["x_2D"].shape[1])
+
+                nc_var = io["nc_file"].createVariable(
+                    varname=varname,
+                    datatype=np.float32,
+                    dimensions=dimensions,
+                    zlib=True,
+                    complevel=4,
+                    fill_value=-9999.0,  # Fill missing values
+                    chunksizes=chunksizes,
+                )
                 # Assign metadata
                 nc_var.units = var_units
                 nc_var.description = var_desc

src/ebfm/config.py

@@ -62,8 +62,9 @@ class GridConfig:
 
     grid_type: GridInputType  # Name of the grid used in coupling
     mesh_file: Path  # Path to the grid file
-    dem_file: Path  # Path to the DEM file (only relevant for CUSTOM grid type)
+    dem_file: Path = None  # Path to the DEM file (only relevant for CUSTOM grid type)
     is_partitioned: bool  # Whether the grid is partitioned
+    is_unstructured: bool = False  # Whether the grid is unstructured
     partition_id: int  # Partition ID (only relevant if is_partitioned is True)
 
     def __init__(self, args: Namespace):
@@ -103,6 +104,11 @@ def __init__(self, args: Namespace):
             self.grid_type = GridInputType.CUSTOM
             self.mesh_file = args.elmer_mesh
             self.dem_file = args.netcdf_mesh
+        elif args.netcdf_mesh_unstructured and args.elmer_mesh:
+            self.grid_type = GridInputType.ELMERXIOS
+            self.mesh_file = args.elmer_mesh
+            self.dem_file = args.netcdf_mesh_unstructured
+            self.is_unstructured = True
         elif args.elmer_mesh:
             self.grid_type = GridInputType.ELMER
             self.mesh_file = args.elmer_mesh

src/ebfm/grid.py

@@ -6,6 +6,14 @@
 
 
 class GridInputType(Enum):
-    MATLAB = "matlab"  # .mat grid file with elevation
-    CUSTOM = "custom"  # Elmer/Ice mesh file for xy-coordinates and separate NetCDF elevation file
-    ELMER = "elmer"  # Elmer/Ice mesh file with elevation in z-coordinate
+    # .mat grid file with elevation
+    MATLAB = "matlab"
+
+    # Elmer/Ice mesh file for xy-coordinates and separate NetCDF elevation file
+    CUSTOM = "custom"
+
+    # Elmer/Ice mesh file for xy-coordinates and separate unstructured NetCDF elevation file obtained from XIOS
+    ELMERXIOS = "elmerxios"
+
+    # Elmer/Ice mesh file with elevation in z-coordinate
+    ELMER = "elmer"

src/main.py

@@ -97,7 +97,7 @@ def extract_active_coupling_features(args: argparse.Namespace) -> List[str]:
 
 
 def main():
-    parser = argparse.ArgumentParser()
+    parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
     parser.add_argument(
         "--version",
@@ -126,20 +126,27 @@ def main():
         " If --netcdf-mesh is provided elevations will be read from the given NetCDF mesh file.",
     )
 
+    input_group.add_argument(
+        "--netcdf-mesh-unstructured",
+        type=Path,
+        help="Path to the unstructured NetCDF mesh file. Optional if using --elmer-mesh."
+        " If --netcdf-mesh is provided elevations will be read from the given NetCDF mesh file.",
+    )
+
     time_group = parser.add_argument_group("time configuration")
 
     time_group.add_argument(
         "--start-time",
         type=str,
-        help="Start time of the simulation in format 'DD-Mon-YYYY HH:MM', e.g., '01-Jan-1979 00:00'.",
+        help="Start time of the simulation in format 'DD-Mon-YYYY HH:MM'",
         default="1-Jan-1979 00:00",
     )
 
     time_group.add_argument(
         "--end-time",
         type=str,
-        help="End time of the simulation in format 'DD-Mon-YYYY HH:MM', e.g., '02-Jan-1979 09:00'.",
-        default="2-Jan-1979 09:00",
+        help="End time of the simulation in format 'DD-Mon-YYYY HH:MM'",
+        default="2-Jan-1979 00:00",
     )
 
     time_group.add_argument(
@@ -297,11 +304,10 @@ def main():
 
         # Write output to files (only in uncoupled run and for unpartitioned grid)
         if not grid["is_partitioned"] and not coupler.has_coupling:
-            assert (
-                grid_config.grid_type is GridInputType.MATLAB
-            ), "Output writing currently only implemented for MATLAB input grids."
-            io, OUTFILE = LOOP_write_to_file.main(OUTFILE, io, OUT, grid, t, time, C)
-            pass
+            if grid_config.grid_type is GridInputType.MATLAB:
+                io, OUTFILE = LOOP_write_to_file.main(OUTFILE, io, OUT, grid, t, time, C)
+            else:
+                logger.warning("Skipping writing output to file for Elmer input grids.")
         elif grid["is_partitioned"] or coupler.has_coupling:
             logger.warning("Skipping writing output to file for coupled or partitioned runs.")
         else:

src/reader.py

@@ -89,6 +89,28 @@ def read_elmer_mesh(mesh_root: Path, is_partitioned: bool = False, partition_id:
     )
 
 
+def read_dem_xios(dem_file: Path, grid: dict):
+    """Read digital elevation model (DEM) file.
+    Args:
+        dem_file (Path): Path to the DEM NetCDF file.
+        grid (dict): grid-related parameters
+    Returns:
+        grid (dict): dictionary containing grid-related parameters
+    """
+    assert dem_file.is_file(), f"DEM file {dem_file} does not exist."
+
+    import netCDF4
+
+    nc = netCDF4.Dataset(dem_file)
+    assert (
+        np.squeeze(nc["x"][:]).shape == grid["x"].shape
+    ), "Surface mesh and Elmer mesh do not have the same number of vertices"
+    grid["z"] = np.squeeze(nc["zs"][:]).data
+    grid["lat"] = np.squeeze(nc["mesh2D_node_x"][:]).data
+    grid["lon"] = np.squeeze(nc["mesh2D_node_y"][:]).data
+    return grid
+
+
 def read_dem(dem_file: Path, xs: NDArray[np.float64], ys: NDArray[np.float64]):
     """Read digital elevation model (DEM) file.
     Args: