** ERROR : network is not orthogonal

[thopri]

Hi,

Thanks for your help previously, I now have a new issue that I could use some guidance for. I am trying a new model domain and I am getting the following error when running dflowfm

** ERROR : network is not orthogonal

I assume this means my grid needs more refinement but I am not sure if there is a way for dfm_tools to do that? Some pointers would be appreciated.

Many Thanks

[thopri]

Here are some screenshots of the grid, I am not that familar with unstructured grids so don't know if this acceptable?

[veenstrajelmer]

Hi @thopri, this is probably because of this issue: Deltares/MeshKernelPy#150

A simple fix would be to use a larger value for area_min for dfmt.meshkernel_delete_withcoastlines() (the unit is square kilometers). This excludes all small islands to be deleted and it often solves this issue. Does it work for you?

[thopri]

Hello! Yes that is the issue, I have applied the mesh2d_get_orthogonality().values.max() method before and after meshkernel_delete_withcoastlines() returning values of 0.00015 before and 0.53695 after. However setting min_area at various values from 0.3 to 10km2 has not fixed it. I can't see any holes like in the issue you linked so am slightly stumped. Does 10km2 sound like a reasonable amount? the cell size is between 1/20th degree at the coarse side and 300m at the fine side.

[thopri]

I am using the high res coastline, maybe using a lower res option would be helpful? The coastline is quite complicated
Update: that made it worse!

[thopri]

Ok, increasing the min area parameter doesn't appear to help, if anything it makes the orthogonality worse. (~0.75). Different resolution of coastlines doesn't seem to have much impact. So any more help would be appreciated!

[thopri]

Just to update, the model runs fine if I don't clip the coastline. I have been measuring the size of some of the holes in the coastline due to the esturaries/rivers and am increasing the min size to remove them. Currently at 250 km2 and 0.4 orthogonality so I think this will work if I get to a certain level. (currently trying 750km2)

Update: an min area of 1500km2 worked :)

[veenstrajelmer]

That is a bit much, but if it works it works.. We are working on this in the issue that I linked, hopefully this will be a bit more charming in the future.

You can plot the edges like this:

from meshkernel import MeshKernel, MakeGridParameters, ProjectionType
import matplotlib.pyplot as plt
plt.close('all')
# import numpy as np

# get dummy mk variable
lon_min, lon_max, lat_min, lat_max = -68.55, -67.9, 11.8, 12.6
dxy = 0.05
make_grid_parameters = MakeGridParameters(
                                        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=ProjectionType.SPHERICAL)
mk.curvilinear_compute_rectangular_grid_on_extension(make_grid_parameters)
mk.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d

# plot
mesh2d = mk.mesh2d_get()
ortho = mk.mesh2d_get_orthogonality().values
# ortho[ortho==-999] = np.nan
fig,ax = plt.subplots()
mesh2d.plot_edges(ax)

It should be fairly simple to get the orthogonality values as line colors, but I have no time to do so now. This would make it easier to see where the cells with issues are.

[thopri]

Oooh helpful thanks! I actually found that 250km2 was sufficent as I thought I would need lower orthogonal values than 0.4 but it turns out the model will at least start with values like that.

[thopri]

Thought I would add some extra experience here, creating a new domain resulted in the issue discussed above, increasing min area didn't help at all in this case. But I found the solution was one I mentioned above, changing the coastline shapefile resolution. The options are f (full), h, (high), i, (intermediate), l, (low) and c, (coarse).

By default I use "h" but I found setting res to "i" resulted in an value low enough for the model to run. "l" was even better but then the coastline was pretty inaccurate. Hope that helps someone (probably me) in the future.

[veenstrajelmer]

There will always be edge cases for which one of these suggestion works or will not work and there are also downsides to using the lower resolution coastlines. I expect that there will be a new meshkernel release next week that can generate illegalcells.pol to use as drypoints in the model, with which this issue will be solved properly. I will try to add an example to dfm_tools asap since I think many users run into this.

[veenstrajelmer]

From dfm_tools 0.24.0 (not released yet), there will be a new function called dfmt.meshkernel_get_illegalcells() that derives polygons about the 6-edged-cells that cause this high orthogonality. The example notebook already includes how to work with this (including the relevant mdu setting). dfm_tools 0.24.0 was not released yet, but you can try it out after the release is done.