high orthogonality with refinement for higher latitudes

[veenstrajelmer]

Describe the bug
Originally reported in Deltares/dfm_tools#1078. The refinement algorithm result in decent grids for arbitrary regions most of the time. However, an edge case was discovered where the base grid has a max orthogonality of 0, but the refined grid as a max orthogonality of >0.99. This is unexpected, I guess (and hope) some inconvenient bug in the algorithm?

To Reproduce

import xarray as xr
import meshkernel

# domain and resolution
model_name = 'Famke'
if model_name=='Bonaire':
    lon_min, lon_max, lat_min, lat_max = -68.55, -67.9, 11.8, 12.6
    # max ortho after refinement: 0.0003921491771577449
elif model_name=='Vietnam':
    lon_min, lon_max, lat_min, lat_max = 105.8, 106.85, 17.75, 18.5
    # max ortho after refinement: 0.000606091522901473
elif model_name=='Famke':
    lon_min, lon_max, lat_min, lat_max =  -1, 2.5, 52.5, 55
    # max ortho after refinement: 0.9999866054173161
dxy = 0.05

# generate spherical regular grid
projection = meshkernel.ProjectionType.SPHERICAL
make_grid_parameters = meshkernel.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_object = meshkernel.MeshKernel(projection=projection)
mk_object.curvilinear_compute_rectangular_grid_on_extension(make_grid_parameters)
mk_object.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d
print('max ortho before refinement:', mk_object.mesh2d_get_orthogonality().values.max())

# get bathymetry dataset
file_nc_bathy = "https://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).elevation
data_bathy_sel = data_bathy.sel(lon=slice(lon_min-1, lon_max+1), lat=slice(lat_min-1, lat_max+1))
lon_np = data_bathy_sel.lon.to_numpy()
lat_np = data_bathy_sel.lat.to_numpy()
values_np = data_bathy_sel.to_numpy().flatten()

# refine grid
min_edge_size = 300 # in meters
gridded_samples = meshkernel.GriddedSamples(x_coordinates=lon_np,y_coordinates=lat_np,values=values_np)
mesh_refinement_parameters = meshkernel.MeshRefinementParameters(
    min_edge_size=min_edge_size, #always in meters
    refinement_type=meshkernel.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_object.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 refinement:', mk_object.mesh2d_get_orthogonality().values.max())

Gives:

max ortho before refinement: 0.0
max ortho after refinement: 0.9999866054173161

Expected behavior
A max orthogonality <0.1, just like for the other cases included.

Screenshots
Figure created in the original issue example code shows where the high orthogonality occurs:

Version info (please complete the following information):

• OS: Windows
• Version: 5.0.2
• Python 3.11

Additional information
The cause is a high aspect ratio of resulting dx/dy (in meters) because of the relatively high latitude of the grid. When supplying block_size_y=dxy * 1.4, the resulting orthogonality is 0.001505. This was actually already reported in #153. Also for DCSM applications, the base grid was originally generated by forcing a factor between dx and dy (dy = 1.5 * dx), so it kind of makes sense. It would be convenient to supply the dx/dy in meters or to be able to quickly compute a preferred aspect ratio between dx and dy. Could this land in meshkernel or should it be external?

Update: if i understood correctly there is auto scaling of dy over latitude. It would just require starting off with a different value, ignoring user input