Non-orthogonal grid and illegal cells function

[n-aleksandrova]

Hi,

I am generating a grid using the latest release of dfm_tools (0.24.0), and using the dfmt.meshkernel_get_illegalcells method to identify non-orthogonal cells. However, I encountered cases where this function was not able to identify illegal cells, which is the case for a specific cell type, see snapshot below (where the blue circle shows the problematic cell that is not picked up by the identification method, and the red rectangle shows another type of problematic cell that is marked automatically as an illegal cell).
I was able to run the model after manually adding a polygon around the cells that caused the issue (my grid had 2 such cells), but would be great if this can be done automatically.

[veenstrajelmer]

Thanks for creating this issue. Would it be possible to share a minimal example that generates this part of the grid? It does make sense that this cell is not found, since it searches for cells with 6 edges, yours has only 4. However, I understand that this cell shape causes an issue, it should have an extra edge dividing it in two cells just like the similar cell a bit above right of this one.

[n-aleksandrova]

Indeed the missing edge is the problem, and I think this happens because of the rather complex coastline in that area. I tried to use the coastline with lower resolution, but that does not change the situation. In this particular grid I found two such cells with missing edges.

Below is the code used to generate the grid, let me know if you are missing any information.

import dfm_tools as dfmt
import xarray as xr

lon_min, lon_max, lat_min, lat_max =  32.55, 47.22, -27.84, -9.6
dxy = 0.1
min_edge_size = 400 # in meters

# generate spherical regular grid
mk_object = dfmt.make_basegrid(lon_min, lon_max, lat_min, lat_max, dx=dxy, dy=dxy, crs=crs)

# generate plifile from grid extent and coastlines
bnd_gdf = dfmt.generate_bndpli_cutland(mk=mk_object, res='h', buffer=0.2)
bnd_gdf_interp = dfmt.interpolate_bndpli(bnd_gdf, res=0.03)

# filter out boundary sections that are very short (<5km in this case)
ids = []
for ii in range(len(bnd_gdf_interp.index)):
    if bnd_gdf_interp.iloc[ii].values.to_crs('EPSG:3857').length < 5000: # approx. 111 km in 1 degree
        ids.append(ii)

bnd_gdf_interp = bnd_gdf_interp.drop(ids)

file_nc_bathy_sel = r'p:\11210471-001-compass\02_Models\Delft3DFM\mozambique_model\boundary_conditions\GEBCO_MZB\gebco_2023_n-8.0_s-29.0_w31.5_e48.5.nc'
data_bathy_sel = xr.open_dataset(file_nc_bathy_sel).elevation
data_bathy_sel.load()

# refine grid
dfmt.refine_basegrid(mk=mk_object, data_bathy_sel=data_bathy_sel, min_edge_size=min_edge_size)

# remove land with GSHHS coastlines
dfmt.meshkernel_delete_withcoastlines(mk=mk_object,res='i',min_area=50,crs=crs)

# derive illegalcells geodataframe
illegalcells_gdf = dfmt.meshkernel_get_illegalcells(mk=mk_object)

# convert to xugrid
xu_grid_uds = dfmt.meshkernel_to_UgridDataset(mk=mk_object, crs=crs)

# interpolate bathymetry onto the grid
data_bathy_interp = data_bathy_sel.interp(lon=xu_grid_uds.obj.mesh2d_node_x, lat=xu_grid_uds.obj.mesh2d_node_y)
xu_grid_uds['mesh2d_node_z'] = data_bathy_interp.clip(max=10)

# write xugrid grid to netcdf
netfile = os.path.join(dir_output_run, f'{model_name}_net.nc')
xu_grid_uds.ugrid.to_netcdf(netfile)

[veenstrajelmer]

Thanks, I have minimized your code so the issue can be easily reproduced:

import dfm_tools as dfmt
import xarray as xr
import matplotlib.pyplot as plt
plt.close("all")

# lon_min, lon_max, lat_min, lat_max =  32.55, 47.22, -27.84, -9.6
lon_min, lon_max, lat_min, lat_max =  46.2, 46.8, -18, -15.85
dxy = 0.1
min_edge_size = 400 # in meters
crs = "EPSG:4326"

# generate spherical regular grid
mk_object = dfmt.make_basegrid(lon_min, lon_max, lat_min, lat_max, dx=dxy, dy=dxy, crs=crs)

# generate plifile from grid extent and coastlines
bnd_gdf = dfmt.generate_bndpli_cutland(mk=mk_object, res='h', buffer=0.2)
bnd_gdf_interp = dfmt.interpolate_bndpli(bnd_gdf, res=0.03)

file_nc_bathy_sel = r'p:\11210471-001-compass\02_Models\Delft3DFM\mozambique_model\boundary_conditions\GEBCO_MZB\gebco_2023_n-8.0_s-29.0_w31.5_e48.5.nc'
data_bathy_sel = xr.open_dataset(file_nc_bathy_sel).elevation
data_bathy_sel = data_bathy_sel.sel(lon=slice(lon_min-0.1,lon_max+0.1),
                                    lat=slice(lat_min-0.1,lat_max+0.1))
data_bathy_sel.load()

# refine grid
dfmt.refine_basegrid(mk=mk_object, data_bathy_sel=data_bathy_sel, min_edge_size=min_edge_size)

# remove land with GSHHS coastlines
# dfmt.meshkernel_delete_withcoastlines(mk=mk_object,res='i',min_area=50,crs=crs)
mesh_bnds = mk_object.mesh2d_get_mesh_boundaries_as_polygons()
bbox = (mesh_bnds.x_coordinates.min(), mesh_bnds.y_coordinates.min(), mesh_bnds.x_coordinates.max(), mesh_bnds.y_coordinates.max())
coastlines_gdf = dfmt.get_coastlines_gdb(bbox=bbox, res='i', min_area=50, crs=crs)
dfmt.meshkernel_delete_withgdf(mk=mk_object, coastlines_gdf=coastlines_gdf)

# derive illegalcells geodataframe
illegalcells_gdf = dfmt.meshkernel_get_illegalcells(mk=mk_object)

# convert to xugrid
xu_grid_uds = dfmt.meshkernel_to_UgridDataset(mk=mk_object, crs=crs)

# add orthogonality
ortho = mk_object.mesh2d_get_orthogonality()
ortho_vals = ortho.values
ortho_vals[ortho_vals==ortho.geometry_separator] = 0
print("ortho max:", ortho_vals.max())
xu_grid_uds['ortho'] = xr.DataArray(ortho_vals, dims=xu_grid_uds.grid.edge_dimension)

# plot NOTE: the order of edges/orthogonality is correctly aligned from xugrid 0.11.0
fig, ax = plt.subplots()
xu_grid_uds['ortho'].ugrid.plot(ax=ax)
coastlines_gdf.plot(ax=ax, color="none", edgecolor='r')

Gives:

Zooming in to the area with the issue shows that the coastline slighly touches the triangle cell below, if this would not have happened, this issue would not have occurred since that cell would also have been deleted:

The issue is that only the two cells are marked as to be deleted by the meshkernel deletion algorithm. This is an edge case that I hope we can avoid in the future. We already worked around it for square cells since they result in cells with six edges (Deltares/MeshKernelPy#150, for which the illegalcells algorithm works properly). However, it seems we also need a fix for triangular cells.

For now, you can manually add illegalcells as you did, but you could also try using res='h' in dfmt.meshkernel_delete_withcoastlines(). It solves this specific cell, but it might still result in other similar cells.

[lucacarniato]

In the upcoming release, a new API will be introduced to retrieve face polygons with all edges having orthogonality values falling within a specified range. For an example, please take a look at the image below.

[veenstrajelmer]

Implemented in Deltares/MeshKernel#367, also call this function in dfm_tools