Unexpectedly high orthogonality for relatively high latitudes

[famkekovacs-deltares]

A new issue has arrived with refining the grid due to a relative high latitude. The initial grid has a strong aspect ratio. The following code creates the issue:

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 =  -1, 2.5, 52.5, 55

dxy = 0.05
min_edge_size = 300 # 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.01)
bnd_gdf_interp = dfmt.interpolate_bndpli(bnd_gdf, res=0.05)

file_nc_bathy_sel = 'https://opendap.deltares.nl/thredds/dodsC/opendap/deltares/Delft3D/netcdf_example_files/GEBCO_2022/GEBCO_2022_coarsefac08.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='h')

# 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())
print('mk.mesh2d_get_orthogonality()')
print(mk_object.mesh2d_get_orthogonality().values.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)
mk_object.mesh2d_get().plot_edges(ax,zorder=1)

In the picture the issue arrives not at the coast as with previous issues but more seaward or completely at sea.

[veenstrajelmer]

This happens because the aspect ratio (in meters) in the base grid is far from 1 already because of the high latitude. A workaround is to use dy=dxy*1.5. Would be great if dy would be automatically scaled for the latitude so this would be avoided, I will see if that is possible in meshkernel or dfm_tools.