From c472074f968c65745d6f0dd2b3e217b0cbfdbc98 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 12 Dec 2025 16:52:35 +0100 Subject: [PATCH 1/7] Remove unnecessary read_elmer_xios_grid. GridInputType.ELMER_XIOS does not allow for independent grid reading of the XIOS NetCDF. The XIOS NetCDF is always the DEM and comes with an additional Elmer grid. --- src/ebfm/INIT.py | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index 49d68f24..67de7342 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -309,33 +309,12 @@ 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. From 4ef27ffcc5abf37f8bb994f72d42209ec0978341 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Mon, 15 Dec 2025 15:57:53 +0100 Subject: [PATCH 2/7] Move relevant code up into DEM part. --- src/ebfm/INIT.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index 67de7342..c4091dca 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -171,14 +171,23 @@ 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["x"] = np.zeros_like(grid["z"]) # remove? 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 + + if grid["input_type"] is GridInputType.ELMERXIOS: + min_thickness_glacier = 1.0 # minimum ice thickness to consider grid cell as glacier (m) + + # treats grid cells as glacier where ice thickness exceeds threshold + grid["mask"] = (grid["h"] > min_thickness_glacier).astype(int) + else: + grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier + + grid["gpsum"] = np.sum(grid["mask"]) # number of glacier grid cells grid["slope_x"] = np.zeros_like(grid["x"]) # test values! grid["slope_y"] = np.zeros_like(grid["x"]) # test values! grid["slope_beta"] = np.zeros_like(grid["x"]) # test values! From 8abcc347214eba0ca78d228fb11d975a2ba0e668 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Mon, 15 Dec 2025 18:25:08 +0100 Subject: [PATCH 3/7] Introduce function for glacier cells. --- src/ebfm/INIT.py | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index c4091dca..caf76770 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -151,6 +151,19 @@ def init_constants(): return C +def compute_number_of_glacier_cells(grid): + """ + Computes the number of glacier cells in the grid based on the mask. + + Parameters: + grid (dict): Dictionary containing grid-related parameters + + Returns: + int: Number of glacier cells (gpsum) + """ + return np.sum(grid["mask"] == 1) + + def init_grid(grid, io, config: GridConfig): grid["is_partitioned"] = config.is_partitioned grid["is_unstructured"] = config.is_unstructured @@ -187,7 +200,7 @@ def init_grid(grid, io, config: GridConfig): else: grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier - grid["gpsum"] = np.sum(grid["mask"]) # number of glacier grid cells + grid["gpsum"] = compute_number_of_glacier_cells(grid) grid["slope_x"] = np.zeros_like(grid["x"]) # test values! grid["slope_y"] = np.zeros_like(grid["x"]) # test values! grid["slope_beta"] = np.zeros_like(grid["x"]) # test values! @@ -214,7 +227,7 @@ def init_grid(grid, io, config: GridConfig): 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["gpsum"] = compute_number_of_glacier_cells(grid) # TODO later add slope # grid["slope_x"], grid["slope_y"] = mesh.dzdy, mesh.dzdy @@ -254,8 +267,7 @@ def init_grid(grid, io, config: GridConfig): # Calculate grid spacing grid["dx"] = grid["x_2D"][0][1] - grid["x_2D"][0][0] - # Calculate number of modeled grid cells (gpsum) - grid["gpsum"] = np.sum(grid["mask_2D"] == 1) + grid["gpsum"] = compute_number_of_glacier_cells(grid) # Create 1-D mask grid["mask"] = grid["mask_2D"][grid["mask_2D"] == 1] From 4a2a308716f849bc6d621cc874dc9c772160b368 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Mon, 15 Dec 2025 18:33:55 +0100 Subject: [PATCH 4/7] Fix --- src/ebfm/INIT.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index caf76770..d34fa2d2 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -267,10 +267,9 @@ def init_grid(grid, io, config: GridConfig): # Calculate grid spacing grid["dx"] = grid["x_2D"][0][1] - grid["x_2D"][0][0] - grid["gpsum"] = compute_number_of_glacier_cells(grid) - # Create 1-D mask grid["mask"] = grid["mask_2D"][grid["mask_2D"] == 1] + grid["gpsum"] = compute_number_of_glacier_cells(grid) # Calculate latitude & longitude fields (from the original UTM coordinates) utmzone = grid["utmzone"] # Assume this is already part of the grid From 9de84de1a1c4f6aeb721d79a393bb0fb7aef4463 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Mon, 15 Dec 2025 18:42:27 +0100 Subject: [PATCH 5/7] Only use mask_2D inside matlab grid reader. --- src/ebfm/INIT.py | 2 +- src/ebfm/LOOP_EBM_insolation.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index d34fa2d2..c8a34e95 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -287,7 +287,7 @@ def init_grid(grid, io, config: GridConfig): grid["x"] = grid["x_2D"].flatten()[mask_flat == 1] grid["y"] = grid["y_2D"].flatten()[mask_flat == 1] grid["z"] = grid["z_2D"].flatten()[mask_flat == 1] - grid["ind"] = np.where(grid["mask_2D"].flatten() == 1) + grid["ind"] = np.where(mask_flat == 1) grid["xind"], grid["yind"] = np.where(grid["mask_2D"] == 1) # --------------------------------------------------------------------- diff --git a/src/ebfm/LOOP_EBM_insolation.py b/src/ebfm/LOOP_EBM_insolation.py index 27d3a6c4..601f80d2 100644 --- a/src/ebfm/LOOP_EBM_insolation.py +++ b/src/ebfm/LOOP_EBM_insolation.py @@ -140,7 +140,7 @@ def main(grid, time2, OUT): # Build the final shade 2D array shade_2D = np.ones((xl, yl), dtype=int) shade_2D.flat[grid["ind"]] = shade - OUT["shade"] = shade_2D.flatten()[grid["mask_2D"].flatten() == 1] + OUT["shade"] = shade_2D.flatten()[grid["ind"]] else: # TODO: Shading based on lookup tables from Elmer/Ice # For now: assume flat surface and no shading by surrounding terrain From 57fc84631e6230db65a3c5bb62e2f4443b297226 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Mon, 15 Dec 2025 18:43:51 +0100 Subject: [PATCH 6/7] Use local variable to not pollute grid dict. --- src/ebfm/INIT.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index c8a34e95..73cb4204 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -242,7 +242,7 @@ def init_grid(grid, io, config: GridConfig): grid["x_2D"] = input_data["x"][0][0] grid["y_2D"] = input_data["y"][0][0] grid["z_2D"] = input_data["z"][0][0] - grid["mask_2D"] = input_data["mask"][0][0] + mask_2D = input_data["mask"][0][0] # Determine domain extent grid["has_shading"] = True @@ -255,20 +255,20 @@ def init_grid(grid, io, config: GridConfig): grid["x_2D"] = np.flipud(grid["x_2D"]) grid["y_2D"] = np.flipud(grid["y_2D"]) grid["z_2D"] = np.flipud(grid["z_2D"]) - grid["mask_2D"] = np.flipud(grid["mask_2D"]) + mask_2D = np.flipud(mask_2D) _, fx = np.gradient(grid["x_2D"]) if fx[0, 0] < 0: grid["x_2D"] = np.fliplr(grid["x_2D"]) grid["y_2D"] = np.fliplr(grid["y_2D"]) grid["z_2D"] = np.fliplr(grid["z_2D"]) - grid["mask_2D"] = np.fliplr(grid["mask_2D"]) + mask_2D = np.fliplr(mask_2D) # Calculate grid spacing grid["dx"] = grid["x_2D"][0][1] - grid["x_2D"][0][0] # Create 1-D mask - grid["mask"] = grid["mask_2D"][grid["mask_2D"] == 1] + grid["mask"] = mask_2D[mask_2D == 1] grid["gpsum"] = compute_number_of_glacier_cells(grid) # Calculate latitude & longitude fields (from the original UTM coordinates) @@ -283,12 +283,12 @@ def init_grid(grid, io, config: GridConfig): grid["lat_2D"] = lat.reshape(grid["y_2D"].shape) # Store 1-D (vectorized) grid information - mask_flat = grid["mask_2D"].flatten() + mask_flat = mask_2D.flatten() grid["x"] = grid["x_2D"].flatten()[mask_flat == 1] grid["y"] = grid["y_2D"].flatten()[mask_flat == 1] grid["z"] = grid["z_2D"].flatten()[mask_flat == 1] grid["ind"] = np.where(mask_flat == 1) - grid["xind"], grid["yind"] = np.where(grid["mask_2D"] == 1) + grid["xind"], grid["yind"] = np.where(mask_2D == 1) # --------------------------------------------------------------------- # Grid slope and aspect From e53c62cdeec4ec366de1bc4f66286e42456d3ac1 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Tue, 16 Dec 2025 07:59:34 +0100 Subject: [PATCH 7/7] Update src/ebfm/INIT.py --- src/ebfm/INIT.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index 73cb4204..c9139ceb 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -184,7 +184,6 @@ 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["x"] = np.zeros_like(grid["z"]) # remove? 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!