Bugfix 1dconnect

[hboisgon]

Issue addressed

Fixes #251

Explanation

There were several issues with the method.

• Bugfixes to connect to a river: subcatch did not always have unique ids + problem with adjacent tributaries cells + unwanted behavior if the area_max threshold is less than the wflow river upa threshold (tributaries gauges would then not end up on the river)
• Connect to several rivers at the same time (see figure below): no more changes than the above bugfixes seemed to be needed. @shartgring maybe you could also still test with your previous example?
• Allow to pass snapping options for the river edges/nodes. I started to add the snap uparea method but as we start from the river (so the whole river link uparea rather than the edge) I think you may get unexpected results... Maybe we should improve the method to allow the user to give either the river geodataframe or the river nodes/edges geodataframe directly. Then uparea snapping would make sense. But as this is a bigger issue, I prefer to leave it until we get use cases or demand from users.

Checklist

• Updated tests or added new tests
• Branch is up to date with main
• Tests & pre-commit hooks pass
• Updated documentation if needed -> method not yet released so I did not update the changelog

Additional Notes (optional)

Pirai case with two rivers

[shartgring]

Hi @hboisgon ! As a heads up, I was reviewing the code today, I expect to have the review finished next week. I was testing the method with the Vechte basin and noticed a few things which are worth mentioning/asking. In the images, the white dotted line is the river gdf, the blue features are the subcatch_1dmodel, red are riv_1dmodel and green dots are gauges_1dmodel

When the provided geometry lies outside of the model domain, the error ValueError("Coordinates outside domain.") gets raised. Maybe this can be prevented by clipping the gdf to the bounds of the model in the workflow? Or checking the bounds of the gdf to the model bounds, and clipping when necessary?

Also related to clipping (as I had to clip my gdf's to get them working), using the WflowModel.bounds or WflowModel.grid.raster.bounds to clip results in incorrect behaviour when using include_river_boundaries=True. The entire model domain gets in subcatch_1dmodelI think this is caused by the bounds depending on the high-res basin and not the basin at model resolution? Using the basin geom to clip seemed to solve it for me. Clipping the gdf by the basin geom might therefore prevent such behaviour?

When area_max is very large (1e6), we would not expect any point inflows right? As all cells will contribute to the 1d_model through subcatch_1dmodel as all upstream areas are smaller than area_max. That is also confirmed by the log. Still, I get some tributary inflows when using such a large value (such as the green dot below). Is that expected behaviour?

I also compared the difference between a gdf without small tributaries, and one with small tributaries (using include_river_boundaries=False to prevent that they get turned into point inflows). The differences are quite larger, especially for a large area_max. Shouldn't the total area of the subcatch_1dmodel be similar in those cases?

The Santa Cruz case works indeed nicely with the changes when using a single geodataframe!

[shartgring]

Nice work Hélène! I did some tests (see also the comments in the PR) and provided a suggestion to include clipping of the rivers, as I did encounter some problems with it in my case. There might be a better way to include the clipping, so I look forward to your thoughts on this

[hboisgon]

Thanks for the review @shartgring ! I fixed the clipping of the river. Could you send me your example model so that I can check what is going wrong if you use a large area max?

[shartgring]

testmodel_vecht.zip
@hboisgon here it is!

[shartgring]

Looks good! Let me know if any changes are added when you have addressed the example model. Looking forward to your comments on that

[hboisgon]

I looked partly at the Vecht and some things I notice is that the river shapefile contains very small inlet for tributaries (see picture) and these are not long enough to be snapped to the right wflow river (without upstream area propertyor a higher res wflow model). Because it's more a snapping issue, I don't think I want to solve it in hydromt-wflow, it's too complex and it's because of wrong snapping (so maybe impossible).

If you remove these small rivers, then things work nicely again apart for the most downstream part of the wflow basin for which no subbasin is built. Maybe this I can still improve in this PR but for the rest, I think it looks good for now for a first version.