Improve orthogonality of rotated spherical grid

Feature Request
To be discussed. Is it possible to generate more orthogonal spherical rotated grids? When we generate a sperical rotated grid like the example below, it results in a max orthogonality of approximately 0.18. When we start with a cartesian grid (orthogonality of approximately 0) and set it on a spherical meshkernel instance, the orthogonality is 0.05

To Reproduce

import meshkernel
import numpy as np
import matplotlib.pyplot as plt
plt.close("all")
import dfm_tools as dfmt

lon_min, lon_max, lat_min, lat_max = 105.7, 106.59, 18.2, 18.46
dxy = 0.02
angle = -43

make_grid_parameters = meshkernel.MakeGridParameters(angle=angle,
                                                     origin_x=lon_min,
                                                     origin_y=lat_min,
                                                     num_columns=45,
                                                     num_rows=45,
                                                     block_size_x=dxy,
                                                     block_size_y=dxy)

fig,(ax1,ax2) = plt.subplots(1,2, figsize=(9,4))

# generate spherical rotated grid: max ortho 0.1642
mk0 = meshkernel.MeshKernel(projection=meshkernel.ProjectionType.SPHERICAL)
mk0.curvilinear_compute_rectangular_grid(make_grid_parameters)
mk0.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d
ortho_values = mk0.mesh2d_get_orthogonality().values
ortho_values[ortho_values == -999] = np.nan
print(np.nanmin(ortho_values), np.nanmax(ortho_values))
mk0.mesh2d_get().plot_edges(ax=ax1)
dfmt.plot_coastlines(ax=ax1, aspect=None, min_area=1000)

# generate cartesian rotated grid: max ortho 1.272853350497083e-12
mk1 = meshkernel.MeshKernel(projection=meshkernel.ProjectionType.CARTESIAN)
mk1.curvilinear_compute_rectangular_grid(make_grid_parameters)
mk1.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d
ortho_values = mk1.mesh2d_get_orthogonality().values
ortho_values[ortho_values == -999] = np.nan
print(np.nanmin(ortho_values), np.nanmax(ortho_values))

# convert to spherical: max ortho 0.0548
mk2 = meshkernel.MeshKernel(projection=meshkernel.ProjectionType.SPHERICAL)
mk2.mesh2d_set(mk1.mesh2d_get())
ortho_values = mk2.mesh2d_get_orthogonality().values
ortho_values[ortho_values == -999] = np.nan
print(np.nanmin(ortho_values), np.nanmax(ortho_values))
mk2.mesh2d_get().plot_edges(ax=ax2)
dfmt.plot_coastlines(ax=ax2, aspect=None, min_area=1000)

Resulting figure:

Reference RGFGRID grid: n:\Projects\11209500\11209671\B. Measurements and calculations\01_Updated_course_materials\02_exercises\Ex01_rectilinear_grd
Main difference is that in these course materials we start with cartesian, then convert to spherical. This also results in an orthogonality of 0.05. Would this also be possible if we start with spherical directly?

Version info (please complete the following information):

• OS: Windows
• Version 3.0.0

Start with the reproduction script and the MeshKernel entry points curvilinear_compute_rectangular_grid, curvilinear_convert_to_mesh2d, mesh2d_set, and mesh2d_get_orthogonality. Compare direct spherical generation with Cartesian generation followed by spherical conversion using the reported parameters. Done would require an agreed approach or documented result that addresses whether the target orthogonality is achievable.