Mesh

[hboisgon]

Issue addressed

Fixes #30

Explanation

Update self.mesh to include both 1D2D. The main changes is that now we only use self.dfmmodel when reading and writing.
This means that conversion functions from xugrid mesh to hydrolib network were added in a new mesh_utils.py script when we
want to use functions from hydrolib-core (eg add branches or create link1d2d etc).

Hopefully this makes it easier in the future to add new parts or data to the delft3dfm mesh and use mesh functions from xugrid.

For now, this PR was co-developed with hydromt core PR 412

Checklist

• Updated tests or added new tests
• Branch is up to date with main
• Tests & pre-commit hooks pass - will be fixed in another PR
• Updated documentation if needed - partly, can be finalized in the docs PR
• Updated changelog.rst if needed - will be fixed in another PR
• update config files in the examples
• optional: check if we can directly read/write mesh using xugrid rather than hydrolib core (would be handy if we start to store actual data like 2D DEM in self.mesh)
• optional: finalize clean up of some functions in dflowfm.py

Additional Notes

I had to put a pin on meshkernel version 2.0.2 as 2.1.0 did some breaking changes on a link1d2d function that hydrolib-core still needs to fix.

As region is now the total bounds of 1D2D mesh, each setup method preparing a 1D or 2D part of the mesh comes with its own region argument so that they can be different from each other. For ease of use, we did not want to set defaults here so region is mandatory for all the setup 1D or 2D mesh methods.

I added two optional tasks to the checklist but they could also be done in another PR.

For the cleanup one, while updating the API docs, I wonder if some of our functions should not be moved into workflows or included in setup_functions ? These are:

• - add_branches (maybe can be combined with set_branches?)
• - add_crosssections
• - get_boundaries (there is not a lot of get methods so maybe can be merged with set_boundaries?)

I think the rest is okay. @xldeltares what do you think?

[hboisgon]

Hi @xldeltares ! I made some progress on adapting the mesh object to 1D2D. Took a bit longer than what I thought... Do you mind reviewing what I did so far?

I updated the tests and they run but fail because the built models are apparently different than the examples ones and I think you're the best person to know if this is expected or not.

Happy to discuss details if you need and thanks in advance :)

[xldeltares]

Hi @hboisgon, nice work! look forward to review it!

Before our discussion, here is a look at the tests (I was able to install the working environment and ran the tests. ):

The example and test model differ in the following files:

• - fm_net.nc
• - nodefile.ini, i.e. manholes

I checked the fm_net.nc, The differences lie in the following variables (first dataset being the example, second dataset being the test, i.e. new results):
mesh2d_edge_nodes: Differences start from the 0th index, with the values in the first dataset being 0 and the corresponding values in the second dataset being 1.
mesh2d_face_nodes: Differences start from the 0th index, with the values in the first dataset being 0 and the corresponding values in the second dataset being 1.
link1d2d: Differences start from the 0th index, with the values in the first dataset being 0 and the corresponding values in the second dataset being 1.
link1d2d_ids: Differences start from the 0th index. The values in the first dataset are NaN, while the corresponding values in the second dataset are 1.
link1d2d_long_names: Differences start from the 0th index. The values in the first dataset are NaN, while the corresponding values in the second dataset are "1".
--> It can be concluded that the differences are mainly due to different means of starting index locations. I don't think this differences matter. :)

When checking nodefile.ini, I noticed that in the tests, i.e. new results, MH_14 did not get an street level sampled correctly from DEM. This is a bit strange because I see nothing changed about this part of the code [line 1415 in delft3dfm.py] (https://github.com/Deltares/hydromt_delft3dfm/pull/38/files#diff-33ac0e158093b8ff9226df60fd435fdbd07769738588ba1cc1400c4385beecd9R1415).
--> I will remark this and conduct debug another time.

[arthurvd]

Hey @xldeltares and @hboisgon,

I did not look into the details of this work and these files, but maybe the following remark helps.
Regarding:

I checked the fm_net.nc, The differences lie in the following variables (first dataset being the example, second dataset being the test, i.e. new results):
mesh2d_edge_nodes: Differences start from the 0th index, with the values in the first dataset being 0 and the corresponding values in the second dataset being 1.

D-Flow FM accepts UGRID-compliant files, and supports both 0- and 1-based indexing: both are valid, as long as the file has consistent data. This concerns the mesh topology variables *face_nodes, *edge_nodes, *edge_faces, *links, etc.
For example, when attribute :start_index=0 in a variable mesh2d_edge_nodes, then the values in that table can start at 0, and a value of "0" would correspond to the first mesh node in all node arrays, for example in the mesh2d_node_x array. I intentionally write "first", because the index of that first one simply depends on the programming language you're using. In Python it would be index 0, in Fortran index 1.
The same file could also be written using 1-based indexing: if :start_index=1 and all values in the topology tables start at 1, then the mesh is entirely equivalent.

[xldeltares]

As a follow up, findings about manholes seems to indicate a difference in behavior of function dem.raster.sample(manholes). At MH_14, the function returns a value of nan, which then passed to fillna, where a 0 value is filled as default.

However, it is unexpected, because all manholes are within the domain of dem, and dem does not have missing values.

As discussed with @hboisgon , it was found that the manholes reprojection in the dem.raster.sample method caused a shift in the location of the manhole, which exceeds the bounds of the dem.

How to fix? Make the region slightly larger than the current total mesh bounds. Will work on that during the review.

[xldeltares]

Hi @arthurvd, thanks for the quick feedback! Good to know that both are valid and how to make the entirely equivalent. :)

[xldeltares]

Hi @hboisgon , I have finished the review and made changes accordingly - the testing of building/updating 1D model, 2D model and 1D2D model are completed. Would you do a quick check with the TODOs at the bottom in mind? Afterwards I think we are done!

Checklist

• Updated tests and examples
• Updated documentation if needed

Notebook

• The example local model is not fully correct yet, and will be further investigated in issue Investigate why 1d mesh node cannot be correctly generated at MH_2 in local test #56 .
• Also the pre-commit is not passed, and will be done in a sperate issue pre-commit fails on ruff check #28

TODOs:

• check the crs property (now we expect the user to always specify the crs as global property) --> @hboisgon
• update the environment to use git installation for hydromt main --> @hboisgon I think you might know better when mesh PR is merged into main
• ci tests pass (needs the correct env for the tests, we could also exclude the linting and doc) --> @hboisgon
• branch is up to date with main (need to merge again when Santacruz #58 is merged to main). --> @xldeltares
• approve the PR (if everything is ok) --> @rhutten I think authors automatically disqualified from being a reviewer

[hboisgon]

Hi @xldeltares. I am now looking into your changes :)
I'll write more later but while reviewing I saw this:

FIXME: question: Do we always need to read and write the mesh? there are updates that are related to geometry changes and not related geometry changes. The latter might not need a read/write.

That's a good question. Can we actually do anything if we don't read the mesh? Right now for self.crs, we try to read if from mesh. Branches and staticgeoms can also only be derived by reading the mesh file (eg for crosssections or boundaries but I am guessing also the structures). The only thing so far we can change without reading the mesh file is changing the config or dimr file. Do you see other cases?

For writting in update mode, if you do not want to re-write the exact same mesh file we could think of adding extra checks functionnality to only write the mesh file if something changed. Note that via the hydromt config file, users can already decide to only write the objects they want (adding any [write_*] to the config like write_config will then tell hydromt it should only write the listed conponents instead of all).

What do you think?

[hboisgon]

Other question in mesh_utils hydrolib_network_from_mesh:

dfm_network._mesh2d._process(
            grids["mesh2d"].mesh
)  # FIXME: test what if the mesh had bedlevel variable

Can the hydrolib network object store variables like bedlevel ? If so we can decide to include it or not. If we are using hydrolib core writers for mesh than I would say yes. If we only use it for network operations then it's not needed per say and avoids large memory usage

[hboisgon]

I reviewed your changes @xldeltares. Very nice and very detailed thanks a lot!! I think we are almost there.
Can you check my small review commit and the two questions above?

The mesh branch of hydromt was merged to main now so we should be able to update the env file temporarily to install from git to get the ci to pass before we merge. Can you have a look at that or shall I?

An extra one question: I am not sure if we want to commit the batch files to run delft3d-fm or if we keep them local?

[xldeltares]

TODO:

• remove the bat file
• open issue about data variables in mesh file: Support data variables in mesh file #71
• open issue about (nice to have feature) allow non-geometry updates components to be completely independent form self properties: Improve Efficiency in Update Mode: Reduce Computational Effort in Reading Model Geometry #72
• check the latest commit
• update env, install hydromt main and hydrolib-core 0.5.2, --> needed to fix llvmlite version for hydrolib-core
• check ci test, if differs, then update examples --> small fix in reading maps.

[xldeltares]

The updated crs funcs look good!

[sonarcloud]

Kudos, SonarCloud Quality Gate passed!   

0 Bugs
0 Vulnerabilities
0 Security Hotspots
12 Code Smells

No Coverage information
1.6% Duplication