diff --git a/CHANGELOG.md b/CHANGELOG.md index f0271ae9..519cd8e0 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/coupling.py b/src/coupling.py index fbbd4ce9..0a2314a7 100644 --- a/src/coupling.py +++ b/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 diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index 0f5c7b74..49d68f24 100644 --- a/src/ebfm/INIT.py +++ b/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. diff --git a/src/ebfm/LOOP_write_to_file.py b/src/ebfm/LOOP_write_to_file.py index 1587eb7d..f7e35df7 100644 --- a/src/ebfm/LOOP_write_to_file.py +++ b/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 diff --git a/src/ebfm/config.py b/src/ebfm/config.py index 8bca5ab6..4a9aa45f 100644 --- a/src/ebfm/config.py +++ b/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 diff --git a/src/ebfm/grid.py b/src/ebfm/grid.py index 5428b33c..9f39dfd6 100644 --- a/src/ebfm/grid.py +++ b/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" diff --git a/src/main.py b/src/main.py index 74345ad0..379ae278 100644 --- a/src/main.py +++ b/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: diff --git a/src/reader.py b/src/reader.py index e29db52f..bdae6c42 100644 --- a/src/reader.py +++ b/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: