prevent high orthogonality after coastline cutting

veenstrajelmer commented 8 months ago

Is your feature request related to a problem? Please describe. In some rare cases, removing grid with coastline polygons deletes only one edge which result in high orthogonality. This is a known issue in workflows from the past, but this is the first time I encounter it in a grid generated with meshkernel.

The below example generates a grid with an orthogonality of max 0.0262, but after cutting it the max orthogonality is 0.418

Example (can probably be simplified to speed it up, but is also sensitive to change):

import matplotlib.pyplot as plt
plt.close('all')
import dfm_tools as dfmt
import xarray as xr
import contextily as ctx

# input params
# lon_min, lon_max, lat_min, lat_max = 142, 159, -46, -31
lon_min, lon_max, lat_min, lat_max = 142, 150, -42, -39
dxy = 0.5 #0.05
crs = 'EPSG:4326'

# grid generation and refinement with GEBCO bathymetry
# file_nc_bathy = r'p:\metocean-data\open\GEBCO\2021\GEBCO_2021.nc'
file_nc_bathy = r'http://opendap.deltares.nl/thredds/dodsC/opendap/deltares/Delft3D/netcdf_example_files/GEBCO_2022/GEBCO_2022_coarsefac08.nc'
data_bathy = xr.open_dataset(file_nc_bathy)
data_bathy_sel = data_bathy.sel(lon=slice(lon_min,lon_max),lat=slice(lat_min,lat_max)).elevation

#TODO: grid generation and bathy-refinement is still to be improved in meshkernel (https://github.com/Deltares/dfm_tools/issues/234)
mk_object = dfmt.make_basegrid(lon_min, lon_max, lat_min, lat_max, dx=dxy, dy=dxy, crs=crs)

fig, ax = plt.subplots()
mk_object.mesh2d_get().plot_edges(ax=ax,linewidth=1)
ctx.add_basemap(ax=ax, crs=crs, attribution=False)
dfmt.plot_coastlines(ax=ax, crs=crs)

#refine
min_edge_size = 500 #in meters
print('Refining basegrid')
dfmt.refine_basegrid(mk=mk_object, data_bathy_sel=data_bathy_sel, min_edge_size=min_edge_size)
print("max ortho after refining", mk_object.mesh2d_get_orthogonality().values.max()) #TODO: couple back to uds, currently ordering mismatch: https://github.com/Deltares/MeshKernelPy/issues/72

fig, ax = plt.subplots()
mk_object.mesh2d_get().plot_edges(ax=ax,linewidth=1)
ctx.add_basemap(ax=ax, crs=crs, attribution=False)
dfmt.plot_coastlines(ax=ax, crs=crs)

#cutcells
print('Cutting cells')
dfmt.meshkernel_delete_withcoastlines(mk=mk_object, res='i')
#TODO: illegalcells.pol necessary?

#TODO: cleanup grid necessary?
# print('mk_object.mesh2d_get_obtuse_triangles_mass_centers()')
# print(mk_object.mesh2d_get_obtuse_triangles_mass_centers().values)
# print('mk_object.mesh2d_get_orthogonality()')
print("max ortho after cutting",mk_object.mesh2d_get_orthogonality().values.max()) #TODO: couple back to uds, currently ordering mismatch: https://github.com/Deltares/MeshKernelPy/issues/72
# print('mk_object.mesh2d_get_hanging_edges()')
# print(mk_object.mesh2d_get_hanging_edges())
# mk_object.mesh2d_delete_hanging_edges()

fig, ax = plt.subplots()
mk_object.mesh2d_get().plot_edges(ax=ax,linewidth=1)
ctx.add_basemap(ax=ax, crs=crs, attribution=False)
dfmt.plot_coastlines(ax=ax, crs=crs)

Describe the solution you'd like Prevent high orthogonality with cutting grid.

Describe alternatives you've considered

it could be considered not to delete these single edges, or not to create cells that have six edges to avoid this issue.
in the past we worked with illegalcells.pol (DryPointsFile), which handled these cells properly. File is written by delft3dfm. Discuss with AvD.
passing the landboundary that was used for cutting to the mdu could potentially solve the problem but is quite a heavy solution.

veenstrajelmer commented 5 months ago

Simplified testcase:

import matplotlib.pyplot as plt
plt.close('all')
# import dfm_tools as dfmt
# import contextily as ctx
import numpy as np
from meshkernel import (ProjectionType, MakeGridParameters, MeshKernel, 
                        MeshRefinementParameters, GriddedSamples, RefinementType,
                        GeometryList, DeleteMeshOption,
                        )

# input params
lon_min, lon_max, lat_min, lat_max = 147.75, 147.9, -40.4, -40.25
dxy = 0.05 #0.05
crs = 'EPSG:4326'

# grid generation and refinement with GEBCO bathymetry
# create base grid
projection = ProjectionType.SPHERICAL
make_grid_parameters = MakeGridParameters(angle=0,
                                            origin_x=lon_min,
                                            origin_y=lat_min,
                                            upper_right_x=lon_max,
                                            upper_right_y=lat_max,
                                            block_size_x=dxy,
                                            block_size_y=dxy)

mk = MeshKernel(projection=projection)
mk.curvilinear_compute_rectangular_grid_on_extension(make_grid_parameters)
mk.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d

#refine
min_edge_size = 500 #in meters
print('Refining basegrid')
gebco_elev = np.array([[-37, -37, -36, -36],
        [-41, -39, -38, -34],
        [-42, -32, -27,  -6],
        [-38,  -6,  -2,   2],
        [-43, -42, -31, -20]],dtype=np.int16)
lat_np = np.array([-40.397917, -40.364583, -40.33125 , -40.297917, -40.264583])
lon_np = np.array([147.76875 , 147.802083, 147.835417, 147.86875 ])
values_np = gebco_elev.flatten()
gridded_samples = GriddedSamples(x_coordinates=lon_np,y_coordinates=lat_np,values=values_np)

#refinement
mesh_refinement_parameters = MeshRefinementParameters(min_edge_size=min_edge_size, #always in meters
                                                                 refinement_type=RefinementType(1), #Wavecourant/1,
                                                                 connect_hanging_nodes=True, #set to False to do multiple refinement steps (e.g. for multiple regions)
                                                                 smoothing_iterations=2,
                                                                 max_courant_time=120)

mk.mesh2d_refine_based_on_gridded_samples(gridded_samples=gridded_samples,
                                           mesh_refinement_params=mesh_refinement_parameters,
                                           use_nodal_refinement=True) #TODO: what does this do?

print("max ortho after refining", mk.mesh2d_get_orthogonality().values.max()) #TODO: couple back to uds, currently ordering mismatch: https://github.com/Deltares/MeshKernelPy/issues/72

#cutcells
print('Cutting cells')
geometry_separator = -999
xx = np.array([147.83625 , 147.839556, 147.855833, 147.877528, 147.904139,
       147.911222, 147.885861, 147.849111, 147.85625 , 147.83625 ])
yy = np.array([-40.305028, -40.301667, -40.293778, -40.30625 , -40.291222,
       -40.301639, -40.326278, -40.324611, -40.309222, -40.305028])
xx = np.concatenate([xx-0.075, [geometry_separator], xx])
yy = np.concatenate([yy, [geometry_separator], yy])
delete_pol_geom = GeometryList(x_coordinates=xx, y_coordinates=yy, geometry_separator=geometry_separator) #TODO: .copy()/to_numpy() makes the array contiguous in memory, which is necessary for meshkernel.mesh2d_delete()
mk.mesh2d_delete(geometry_list=delete_pol_geom, 
                 delete_option=DeleteMeshOption.INSIDE_NOT_INTERSECTED,
                 invert_deletion=False)

ortho = mk.mesh2d_get_orthogonality().values
# print('mk_object.mesh2d_get_obtuse_triangles_mass_centers()')
# print(mk_object.mesh2d_get_obtuse_triangles_mass_centers().values)
# print('mk_object.mesh2d_get_orthogonality()')
print("max ortho after cutting",ortho.max()) #TODO: couple back to uds, currently ordering mismatch: https://github.com/Deltares/MeshKernelPy/issues/72
# print('mk_object.mesh2d_get_hanging_edges()')
# print(mk_object.mesh2d_get_hanging_edges())
# mk_object.mesh2d_delete_hanging_edges()

m2d = mk.mesh2d_get()
bool_ortho_high = ortho>0.4
sel_x = m2d.edge_x[bool_ortho_high]
sel_y = m2d.edge_y[bool_ortho_high]
sel_o = ortho[bool_ortho_high]

fig, ax = plt.subplots()
m2d.plot_edges(ax=ax,linewidth=1)
ax.scatter(sel_x,sel_y,5,c=sel_o, zorder=5)
# ctx.add_basemap(ax=ax, crs=crs, attribution=False)
# dfmt.plot_coastlines(ax=ax, res='i', min_area=5, crs=crs)

Results in:

max ortho after refining 0.0023339926563505273
max ortho after cutting 0.5525717808005531

lucacarniato commented 4 months ago

closed with latest meshkernelpy pr

veenstrajelmer commented 4 months ago

Adding to the code in https://github.com/Deltares/MeshKernelPy/issues/150#issuecomment-2076921746, one can now use the new mesh2d_get_face_polygons() API to get a geometrylist around the "illegal" cell to use as a drypointsfile in the FM mdu:

illegalcells_geom = mk.mesh2d_get_face_polygons(num_edges=6)
illegalcells_np = np.c_[illegalcells_geom.x_coordinates, illegalcells_geom.y_coordinates]
illegalcells_np[illegalcells_np==illegalcells_geom.geometry_separator] = np.nan

fig, ax = plt.subplots()
m2d.plot_edges(ax=ax,linewidth=1)
# ax.scatter(sel_x,sel_y,5,c=sel_o, zorder=5)
ax.plot(illegalcells_np[:,0], illegalcells_np[:,1], "r-")

Results in:

Deltares / MeshKernelPy

prevent high orthogonality after coastline cutting #150