Smoother refinement along coastlines with bathy-based refinement

[veenstrajelmer]

Is your feature request related to a problem? Please describe.
When using bathymetry based refinement, the resolution along the coastline is not always smooth. This might be due to the lack of a smoothiterations parameter, but maybe an additional algorithm should be implemented to find coastlines. The latter was implemented in the dflowfm kernel in UNST-5236.

Updated to smaller testcase:

Describe the solution you'd like
Smooth refinmenent along coastlines (less triangles)

Additional context
MWE:

import meshkernel
import xarray as xr
import matplotlib.pyplot as plt
plt.close('all')
import numpy as np
import contextily as ctx

#general settings
lon_min,lon_max = -1,-0.2
lat_min,lat_max = 49.1,49.6
lon_res,lat_res = 0.1,0.1
num_x = int(np.ceil((lon_max-lon_min)/lon_res))
num_y = int(np.ceil((lat_max-lat_min)/lat_res))

figsize = (10,4)
crs = 'EPSG:4326'

"""
Make a regular (potentially rotated) rectilinear grid. First generate a curvilinear grid than convert the curvilinear grid into unstructured grid. The steps are the following:
- curvilinear_make_uniform, see the following notebook: https://github.com/Deltares/MeshKernelPy/blob/AddCurvilinearGridSupport/docs/examples/04_curvilineargrid_basics.ipynb
- curvilinear_convert_to_mesh2d: https://github.com/Deltares/MeshKernelPy/blob/118cb4953c4e95d5b18ed283bb37f391134b2bb2/meshkernel/meshkernel.py#L1399 
"""

# Create an instance of MakeGridParameters and set the values
make_grid_parameters = meshkernel.MakeGridParameters()
make_grid_parameters.num_columns = num_x
make_grid_parameters.num_rows = num_y
make_grid_parameters.angle = 0.0
#make_grid_parameters.block_size = 0.0
make_grid_parameters.origin_x = lon_min
make_grid_parameters.origin_y = lat_min
make_grid_parameters.block_size_x = lon_res
make_grid_parameters.block_size_y = lat_res

pol_x = np.empty(0, dtype=np.double)
pol_y = np.empty(0, dtype=np.double)
geometry_list = meshkernel.GeometryList(pol_x, pol_y)

mk2 = meshkernel.MeshKernel()
mk2.curvilinear_make_uniform(make_grid_parameters, geometry_list)
mk2.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d

"""
Mesh refinement in MeshKernelPy with bathymetry samples and plot result
"""
#select and plot bathy
file_nc_bathy = r'p:\metocean-data\open\GEBCO\2021\GEBCO_2021.nc'
data_bathy = xr.open_dataset(file_nc_bathy)
data_bathy_sel = data_bathy.sel(lon=slice(lon_min-1/4,lon_max+1/4),lat=slice(lat_min-1/4,lat_max+1/4))

if 0: #write coarse grid and cutout of gebco asc
    import dfm_tools as dfmt
    uds_mk2 = dfmt.meshkernel_to_UgridDataset(mk2)
    uds_mk2.ugrid.to_netcdf('base_net.nc')
    # data_bathy_sel.to_netcdf('GEBCO_2021_cutout.nc')
    # lon_grid,lat_grid = np.meshgrid(data_bathy_sel.lon,data_bathy_sel.lat)
    # gebco_xyz = np.c_[lon_grid.ravel(),lat_grid.ravel(),data_bathy_sel.elevation.to_numpy().ravel()]
    # np.savetxt(fname='GEBCO_2021_cutout.xyz',X=gebco_xyz,fmt='%16.12f %16.12f %5.1f')
    dfmt.write_bathy_toasc(filename_asc='GEBCO_2021_cutout.asc',lon_sel_ext=data_bathy_sel.lon,lat_sel_ext=data_bathy_sel.lat,elev_sel_ext=data_bathy_sel.elevation)#,asc_fmt='%9.2f',nodata_val=32767)

#convert bathy data to geomlist
samp_x,samp_y = np.meshgrid(data_bathy_sel.lon.to_numpy(),data_bathy_sel.lat.to_numpy())
samp_z = data_bathy_sel.elevation.to_numpy().astype(float) 
samp_x = samp_x.ravel()
samp_y = samp_y.ravel()
samp_z = samp_z.ravel()
geomlist = meshkernel.GeometryList(x_coordinates=samp_x, y_coordinates=samp_y, values=samp_z)

#refinement
mesh_refinement_parameters = meshkernel.MeshRefinementParameters(refine_intersected=False,
                                                                 use_mass_center_when_refining=False, 
                                                                 min_face_size=0.01, 
                                                                 refinement_type=meshkernel.RefinementType(1), #Wavecourant/1,
                                                                 connect_hanging_nodes=True, #set to False to do multiple refinement steps (e.g. for multiple regions)
                                                                 account_for_samples_outside_face=True, #outsidecell argument for --refine?
                                                                 max_refinement_iterations=5,
                                                                 )

mk2.mesh2d_refine_based_on_samples(samples=geomlist,
                                   relative_search_radius=0.5, 
                                   minimum_num_samples=3,
                                   mesh_refinement_params=mesh_refinement_parameters,
                                   )

mesh2d_grid2 = mk2.mesh2d_get()

#zoomed in plot to focus on patchy coastlines
fig, ax = plt.subplots(figsize=figsize)
mesh2d_grid2.plot_edges(ax,linewidth=1)
ctx.add_basemap(ax=ax, crs=crs, attribution=False)

[veenstrajelmer]

When using the dflowfm executable in the command line, a command like this is uses:
/p/d-hydro/dimrset/weekly/2.23.03.78224/lnx64/bin/run_dflowfm.sh --refine:hmin=0.01:directional=0:outsidecell=1:connect=1:smoothiters=5:jasfer3d=0 base_net.nc GEBCO_2021_cutout.asc

Resulting in this:

Some differences:

• the parameter dtmax seems to be missing from meshkernel.MeshRefinementParameters.
• the parameter hmin is in meters, while min_face_size is in degrees in case of a spherical model. This could be because is_geographic is not implemented yet and therefore not used here (Generating a spherical grid raises error #39).
• the parameter hmin is an edge length, while min_face_size is a surface. This does not matter too much, but it is different.
• the dflowfm executable performs bilinear interpolation on the gridded dataset, while meshkernel currently only accepts xyz formatted data (Add support for gridded datasets (and bilinear interpolation option) #34).
• the dflowfm executable searches a landsea interface implemented in [UNST-5236](https://issuetracker.deltares.nl/browse/UNST-5236

base_net.nc.txt
GEBCO_2021_cutout.asc.txt
base_spherical_net.nc.txt

[veenstrajelmer]

When using the latest whl, is_geographic is implemented (#39), but when doing mk2 = meshkernel.MeshKernel(is_geographic=True), this example script raises OSError: exception: access violation reading 0xFFFFFFFFFFFFFFFF upon refinement

[veenstrajelmer]

Using the latest whl from the release branch, the base_grid does not extend to 49.6 latitude, the block_size_y seems to be 0.05 instead of 0.1, also visible in block [4] in the example notebook:

This is intented behaviour to keep the aspect ratio close to squares in spherical coordinates, but might be better to ask different user input in this case.

Also, when adjusting the num_y to num_y*2, the refinement does not match the coastline, it is shifted southwards:

There are other issues when using varying settings, which are documented in TODO comments in the below example code. I propose to focus on them in a next iteration.

import meshkernel
import xarray as xr
import matplotlib.pyplot as plt
plt.close('all')
import numpy as np
import contextily as ctx

is_geographic = True #TODO: is_geographic implemented in release branch
gridded_samples = True #TODO: normal samples take ages and give completely different result, gridded samples are fast but refinement is south of actual coastline

#general settings
lon_min,lon_max = -1,-0.2
lat_min,lat_max = 49.1,49.6
lon_res,lat_res = 0.1,0.1
num_x = int(np.ceil((lon_max-lon_min)/lon_res))
num_y = int(np.ceil((lat_max-lat_min)/lat_res))
if is_geographic:
    num_y = num_y*2 #TODO: remove *2, necessary to get correct lat grid extent with is_geographic=True

figsize = (10,4)
crs = 'EPSG:4326'

"""
Make a regular (potentially rotated) rectilinear grid. First generate a curvilinear grid than convert the curvilinear grid into unstructured grid. The steps are the following:
- curvilinear_make_uniform, see the following notebook: https://github.com/Deltares/MeshKernelPy/blob/AddCurvilinearGridSupport/docs/examples/04_curvilineargrid_basics.ipynb
- curvilinear_convert_to_mesh2d: https://github.com/Deltares/MeshKernelPy/blob/118cb4953c4e95d5b18ed283bb37f391134b2bb2/meshkernel/meshkernel.py#L1399 
"""

# Create an instance of MakeGridParameters and set the values
make_grid_parameters = meshkernel.MakeGridParameters()
make_grid_parameters.num_columns = num_x
make_grid_parameters.num_rows = num_y
make_grid_parameters.angle = 0.0
#make_grid_parameters.block_size = 0.0
make_grid_parameters.origin_x = lon_min
make_grid_parameters.origin_y = lat_min
make_grid_parameters.block_size_x = lon_res
make_grid_parameters.block_size_y = lat_res

pol_x = np.empty(0, dtype=np.double)
pol_y = np.empty(0, dtype=np.double)
geometry_list = meshkernel.GeometryList(pol_x, pol_y)

mk2 = meshkernel.MeshKernel(is_geographic=is_geographic)
mk2.curvilinear_make_uniform(make_grid_parameters, geometry_list)
mk2.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d

mesh2d = mk2.mesh2d_get()
fig, ax = plt.subplots()
mesh2d.plot_edges(ax)
ctx.add_basemap(ax=ax, crs=crs, attribution=False)

"""
Mesh refinement in MeshKernelPy with bathymetry samples and plot result
"""
#select and plot bathy
file_nc_bathy = r'p:\metocean-data\open\GEBCO\2021\GEBCO_2021.nc'
data_bathy = xr.open_dataset(file_nc_bathy)
data_bathy_sel = data_bathy.sel(lon=slice(lon_min-1/4,lon_max+1/4),lat=slice(lat_min-1/4,lat_max+1/4))

if gridded_samples:
    lon_np = data_bathy_sel.lon.to_numpy()
    lat_np = data_bathy_sel.lat.to_numpy()
    values_np = data_bathy_sel.elevation.to_numpy().flatten().astype('float') #TODO: astype to avoid "TypeError: incompatible types, c_short_Array_74880 instance instead of LP_c_double instance"
    gridded_samples = meshkernel.GriddedSamples(n_cols=data_bathy_sel.elevation.shape[1],#-1,
                                                n_rows=data_bathy_sel.elevation.shape[0],#-1,
                                                x_origin=lon_np.min(),
                                                y_origin=lat_np.min(),
                                                cell_size=0.0041666666666, #TODO: make more precise, prefferably by providing all lat/lon values,
                                                values=values_np)
    #gridded_samples = meshkernel.GriddedSamples(x_coordinates=lon_np,y_coordinates=lat_np,values=values_np) #TODO: does not result in refinement
else:
    #convert bathy data to geomlist
    samp_x,samp_y = np.meshgrid(data_bathy_sel.lon.to_numpy(),data_bathy_sel.lat.to_numpy())
    samp_z = data_bathy_sel.elevation.to_numpy().astype(float) 
    samp_x = samp_x.ravel()
    samp_y = samp_y.ravel()
    samp_z = samp_z.ravel()
    geomlist = meshkernel.GeometryList(x_coordinates=samp_x, y_coordinates=samp_y, values=samp_z)

#refinement
if is_geographic: #min_edge_size in meters, so is different for is_geographic=True/False
    min_edge_size = 500
else:
    min_edge_size=0.01
mesh_refinement_parameters = meshkernel.MeshRefinementParameters(refine_intersected=False,
                                                                 use_mass_center_when_refining=False,
                                                                 min_edge_size=min_edge_size, 
                                                                 refinement_type=meshkernel.RefinementType(1), #Wavecourant/1,
                                                                 connect_hanging_nodes=True, #set to False to do multiple refinement steps (e.g. for multiple regions)
                                                                 account_for_samples_outside_face=False, #outsidecell argument for --refine?
                                                                 max_refinement_iterations=5,
                                                                 smoothing_iterations=5,
                                                                 max_courant_time=120.0, #TODO: important to add?
                                                                 directional_refinement=0,
                                                                 )

if gridded_samples:
    mk2.mesh2d_refine_based_on_gridded_samples(gridded_samples=gridded_samples,
                                               mesh_refinement_params=mesh_refinement_parameters,
                                               use_nodal_refinement=True)
else:
    mk2.mesh2d_refine_based_on_samples(samples=geomlist,
                                       relative_search_radius=0.5, 
                                       minimum_num_samples=3,
                                       mesh_refinement_params=mesh_refinement_parameters,
                                       )

mesh2d_grid2 = mk2.mesh2d_get()

#zoomed in plot to focus on patchy coastlines
fig, ax = plt.subplots(figsize=figsize)
mesh2d_grid2.plot_edges(ax,linewidth=1)
ctx.add_basemap(ax=ax, crs=crs, attribution=False)

[lucacarniato]

The number of columns in a gridded sample set is data_bathy_sel.elevation.shape[1] - 1. The number of rows is data_bathy_sel.elevation.shape[0] - 1. Once n_cols and n_cols are corrected, the refinement occurs at the land/sea interface as expected.

[veenstrajelmer]

Tested with whl from 29-06-2023.

• reference is all three booleans from below MWE set to True
• when setting gridded_samples=False, refining takes ages (known issue)
• when setting is_geographic=False and gridded_samples=False, the refinement looks odd. (image below)
• when setting numrowcols=False (using upper_right_x and upper_right_y instead), the resulting grid is of a too small extent (image below) >> setting num_columns=0 and num_rows=0 solves this, preferrably make these the meshkernel default values.

Odd refinement:

Wrong extent:

Updated and simplified MWE:

import meshkernel
import xarray as xr
import matplotlib.pyplot as plt
plt.close('all')
import numpy as np
import contextily as ctx

is_geographic = True #TODO: False results in weird refinement pattern
gridded_samples = True #TODO: False (so normal samples) takes ages and for is_geographic=True and give completely different result for is_geographic=False
numrowcols = True #TODO: False results in grid of too small extent if num_columns/num_rows are not set to 0

#general settings
lon_min,lon_max = -1,-0.2
lat_min,lat_max = 49.1,49.6
lon_res,lat_res = 0.1,0.1

figsize = (10,4)
crs = 'EPSG:4326'

# Create an instance of MakeGridParameters and set the values
make_grid_parameters = meshkernel.MakeGridParameters()
make_grid_parameters.angle = 0.0
make_grid_parameters.origin_x = lon_min
make_grid_parameters.origin_y = lat_min
if numrowcols:
    num_x = int(np.ceil((lon_max-lon_min)/lon_res))
    num_y = int(np.ceil((lat_max-lat_min)/lat_res))
    if is_geographic:
        num_y = num_y*2 #TODO: remove *2, necessary to get correct lat grid extent with is_geographic=True
    make_grid_parameters.num_columns = num_x
    make_grid_parameters.num_rows = num_y
else:
    make_grid_parameters.num_columns = 0 #TODO: this is necessary or the grid extent will be wrong
    make_grid_parameters.num_rows = 0
    make_grid_parameters.upper_right_x = lon_max
    make_grid_parameters.upper_right_y = lat_max
make_grid_parameters.block_size_x = lon_res
make_grid_parameters.block_size_y = lat_res

mk2 = meshkernel.MeshKernel(is_geographic=is_geographic)
mk2.curvilinear_make_uniform(make_grid_parameters)
mk2.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d

mesh2d = mk2.mesh2d_get()
fig, ax = plt.subplots()
mesh2d.plot_edges(ax)
ctx.add_basemap(ax=ax, crs=crs, attribution=False)

#select and plot bathy
file_nc_bathy = r'p:\metocean-data\open\GEBCO\2021\GEBCO_2021.nc'
data_bathy = xr.open_dataset(file_nc_bathy)
data_bathy_sel = data_bathy.sel(lon=slice(lon_min-1/4,lon_max+1/4),lat=slice(lat_min-1/4,lat_max+1/4))

if gridded_samples:
    lon_np = data_bathy_sel.lon.to_numpy()
    lat_np = data_bathy_sel.lat.to_numpy()
    values_np = data_bathy_sel.elevation.to_numpy().flatten().astype('float') #TODO: astype to avoid "TypeError: incompatible types, c_short_Array_74880 instance instead of LP_c_double instance"
    gridded_samples = meshkernel.GriddedSamples(x_coordinates=lon_np,y_coordinates=lat_np,values=values_np) #TODO: does not result in refinement
else:
    #convert bathy data to geomlist
    samp_x,samp_y = np.meshgrid(data_bathy_sel.lon.to_numpy(),data_bathy_sel.lat.to_numpy())
    samp_z = data_bathy_sel.elevation.to_numpy().astype(float) 
    samp_x = samp_x.ravel()
    samp_y = samp_y.ravel()
    samp_z = samp_z.ravel()
    geomlist = meshkernel.GeometryList(x_coordinates=samp_x, y_coordinates=samp_y, values=samp_z)

#refinement
if is_geographic: #min_edge_size in meters, so is different for is_geographic=True/False
    min_edge_size = 500
else:
    min_edge_size = 0.005
mesh_refinement_parameters = meshkernel.MeshRefinementParameters(refine_intersected=False,
                                                                 use_mass_center_when_refining=False,
                                                                 min_edge_size=min_edge_size, 
                                                                 refinement_type=meshkernel.RefinementType(1), #Wavecourant/1,
                                                                 connect_hanging_nodes=True, #set to False to do multiple refinement steps (e.g. for multiple regions)
                                                                 account_for_samples_outside_face=False, #outsidecell argument for --refine?
                                                                 max_refinement_iterations=5,
                                                                 smoothing_iterations=5,
                                                                 max_courant_time=120.0, #TODO: important to add?
                                                                 directional_refinement=0,
                                                                 )

if gridded_samples:
    mk2.mesh2d_refine_based_on_gridded_samples(gridded_samples=gridded_samples,
                                               mesh_refinement_params=mesh_refinement_parameters,
                                               use_nodal_refinement=True)
else:
    mk2.mesh2d_refine_based_on_samples(samples=geomlist,
                                       relative_search_radius=0.5, 
                                       minimum_num_samples=3,
                                       mesh_refinement_params=mesh_refinement_parameters,
                                       )

mesh2d_grid2 = mk2.mesh2d_get()

#plot
fig, ax = plt.subplots(figsize=figsize)
mesh2d_grid2.plot_edges(ax,linewidth=1)
ctx.add_basemap(ax=ax, crs=crs, attribution=False)

[veenstrajelmer]

Tested with wheel of 4-7-2023.

Findings:

• refinement based on gridded samples did not stop when min_edge_size was reached >> solved
• difficulties in creating grid for lat/lon extent, solved by introducing curvilinear_make_uniform_on_extension
• user feedback to be improved, follow up issue: Improve curvilinear_make_uniform user feedback #70
• refinement with gridded or non-gridded samples results in different results since the landsea refinement is only implemented for gridded samples, follow up issue: Align refinement on gridded/nongridded samples #71