Add support for reading ISG files

[Manangka]

In GitLab by @JoerivanEngelen on Feb 6, 2023, 17:31

From what I know, the ISG file format follows a fairly complicated data model.

It is however used in the regional models of the Netherlands, for the SFR package, as far as I know. (Not for the LHM).

[JoerivanEngelen]

Notes Part 1

I sat with @Huite to discuss possible approaches a while ago. He worked on a branch, but never got to finishing it.

My notes are as follows:

• Investigate what ISG functionality is actually required.
• The iMOD5 manual describes all ISG options
• ASFR = 0 is standard ISG, ASFR = 1 requires graph. Don't support ASFR=1.

Segment definition

ISP   ISD
|     |
+-----x-----x-----x-----+

ISD: Calculation points
ISP: Segment definition points

Files

Basically:

• ISP: segment definitions
• ISD: calculation points, contains waterlevels
• ISC: cross-sections

More specifically:

• ISP: contains xy coordinates segment
• ISD1: calculation point locations
• ISD2: timeseries
• ISC1: cross-sections settings
• ISC2: definition cross-sections

ISC1

N > 0 : cross-section
N < 0 : bathymetry

DIST: distance along segments

ISC2

DIST: Distance of the cross-section measured from the centre of the riverbed (minus to the left, and positive to the right)

N > 0

      DIST       DIST
   |<------>|<---------->|
   |        |            x
  -x-       |  x        /
 x   \      | / \      /
/     \     |x   \    /
       x   /|     x--x
        --x |      

N < 0

    +---+---+---+
    |   |   |   |
+---+---+---+---+
|   |   |   |   |
+---+---+---+---+
    |   |   |   |
    +---+---+---+
            | O |
            +---+

Data which defines cell O:

if dx > 0 and dy > 0:

• (x, y, z)
• where: z = 32-bit float

if dx < 0 and dy < 0

• (x, y, z_m, z_c, z_p),
• where: z_m = 16-bit int, z_c = 8-bit int, z_p = 8-bit int

IST

For weirs

w_level_up  ->
------------|
            |
            | w_level_down
            |---------------

ISQ

For Q-H relations

Tabulated rating curve (just like RIBASIM)

Q, width, depth, factor

[JoerivanEngelen]

Notes part 2

Gridding bathymetry (N<0)

Say for the following 2D bathymetry grid:

    +---+---+---+
    |   |   |   |
+---+---+---+---+
|   |   |   |   |
+---+---+---+---+
    |   |   |   |
    +---+---+---+
            |   |
            +---+

We have the following line orthogonal to the segment at calculation point X:

    +---+---+---+
    | 4\|   |   |
+---+-------+---+
|   |   |\3 |   |
+---+---+-X-+---+
    |   | 2\|   |
    +---+---+---+
            |\1 |
            +---+

With the following cross-sections and water-level h at the calculation point:

|                       |
|                       |
|                       |
x-----x                 |
   1  |~~~~~~~~~~~~~~~~~|<-- h
      |                 |
      x-----x     x-----x
         2  |     |  4
            x-----x
               3

The bathymetry can be gridded to a Ugrid2D. The orthogonal lines at calculation points can be used to map stages at calcution points to bathymetry grid cells. The grid can be fully filled with a Laplace interpolation. All cells active where h>bottom, else make nan

[JoerivanEngelen]

Notes part 3

Gridding cross-sections (N>0)

We discussed rasterrizing the cross-sections by using a buffer of their maximum width, but this idea was discarded because of the many (literal) corner cases.

For each profile, create a table with the wetted area/perimeter for each stage.
Conduct piecewise linear interpolation for inbetween values.

  h  | A
-----------
0.1  | 1.0
0.2  | 2.2
0.3  | 3.1
....

This can be used to compute the conductance. River bottom elevation is determined by simply taking the minimum elevation. Stage is defined by calculation points ISG.

In this way, we can keep things to rasterizing a 1D Ugrid network, as explained further.

For the following segment laying on a grid:

                       O
                      /
+------+------+------+
|      |      |     /|
|      |      |    / |
+------+------+---/--+
|      O------|--O   |
|     /|      |      |
+------+------+------+
|   /  |      |      |
|  O   |      |      |
+------+------+------+

Then we add extra points on the 1D network where it intersects the grid edges (the numbers indicate cellids for the first column):

                       O
                      /
+------+------+------X
|0     |      |     /|
|      |      |    / |
+------+------+---X--+
|3     O------X--O   |
|     /|      |      |
+----X-+------+------+
|6  /  |      |      |
|  O   |      |      |
+------+------+------+

The new intersection points can be filled with Laplace interpolation. With this, we can work towards the following table, which is used to compute conductances:

segment|cellid|length|h_int|wetted_fraction
-----------------------------------------
0      | 6    | 0.75 |  .  |     .
1      | 3    | 0.75 |  .  |     .
2      | 4    | 1.5  |  .  |     .
3      | 5    | 0.5  |  .  |     .
4      | 5    | 0.25 |  .  |     .
5      | 2    | 1.25 |  .  |     .

For this we require tools to refine networks, based on vertices. Care should be taken here with the indexes.

0               1    	    2
^               ^    	    ^
O---X---X---X---O---X---X---O
^   ^   ^   ^   ^   ^   ^   ^
0   1   2   3   4   5   6   7

Should become:

O---X---X---X---O---X---X---O
^   ^   ^   ^   ^   ^   ^   ^
0   4   5   6   2   7   8   3

This requires the following features:

• xugrid: Laplace interpolation. Done in Deltares/xugrid@90c03be
• xugrid: refine 1D networks, see Add a refine_with_vertices method to Ugrid1d xugrid#191
• numba_celltree: EdgeCellTree, see Add EdgeCellTree2d class numba_celltree#11
• iMOD Python: New Ugrid 1D
• iMOD Python: Create new empty dataset
• iMOD Python: Index with new edge nodes
• iMOD Python: Where index != -1 or where have data, put data
• iMOD Python: Laplace interpolation
• iMOD Python: Assign 1D network with tables to river grid
• iMOD Python: Read D-Hydro files

[JoerivanEngelen]

@Huite and I had a meeting with Wouter. It was mentioned there that:

• The Bathymetry grids are used mainly for big rivers in Moria & Ibrahym
• Big rivers with 1D cross-section profiles rarely occur

[Huite]

Laplace interpolation for 1D has already been fixed a little while ago: Deltares/xugrid@90c03be

[JoerivanEngelen]

Thanks, updated the list with issues to track in the other repositories.

[Huite]

Xugrid additions to do this cleanly:

• Ugrid1d.refine_with_vertices to insert intersection vertices into the river/stream network; at those sites we can assign properties.
• Laplace interpolation to fill the entire refined topology (should work)
• Translate data from nodes to edges; associate differents facets: Associate data from one aspect (node, face, edge) to another aspect (node, face, edge) xugrid#173
  (not very important because converting nodes to edges is so simple)
• Translation of the filled network to the actual 2D model grid. This can already be done via Ugrid2d.intersect_edges + custom pandas groupby aggregation, but nicer is a "Gridder" which translates LineStrings, Polygons, Ugrid1d to Ugrid2d; comparable to burn_polygons or burn_linestrings, but uses actual overlap and supports the different overlap methods. PolygonGridder, LineStringGridder, NetworkGridder.