Add support for reading ISG files #326

Open Manangka opened 1 year ago

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).

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

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

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:

[x] xugrid: Laplace interpolation. Done in https://github.com/Deltares/xugrid/commit/90c03be9d9ac6edd4a6045c3ca50a86714dfd4de
[ ] xugrid: refine 1D networks, see https://github.com/Deltares/xugrid/issues/191
[ ] numba_celltree: EdgeCellTree, see https://github.com/Deltares/numba_celltree/issues/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

@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

Laplace interpolation for 1D has already been fixed a little while ago: https://github.com/Deltares/xugrid/commit/90c03be9d9ac6edd4a6045c3ca50a86714dfd4de

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

Deltares / imod-python