Skip to content

DEM#

DEM#

hapi.dem.DEM #

Bases: Dataset

Digital Elevation Model dataset with flow-direction helpers.

DEM wraps a GDAL-backed raster dataset (via pyramids.dataset.Dataset) and adds methods that convert D8 flow-direction codes into downstream-cell indices and upstream-cell lookup tables. Three encodings are supported:

  • ESRI (default): powers-of-2 codes (1, 2, 4, 8, 16, 32, 64, 128).
  • SAGA: codes 0--7, starting East counter-clockwise.
  • GRASS: codes 1--8, starting North clockwise.

Use DEM.read_file(path) to open a raster from disk, exactly like pyramids.dataset.Dataset.

Source code in src/hapi/dem.py 
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
class DEM(Dataset):
    """Digital Elevation Model dataset with flow-direction helpers.

    ``DEM`` wraps a GDAL-backed raster dataset (via
    ``pyramids.dataset.Dataset``) and adds methods that convert
    D8 flow-direction codes into downstream-cell indices and
    upstream-cell lookup tables.  Three encodings are supported:

    - **ESRI** (default): powers-of-2 codes
      (1, 2, 4, 8, 16, 32, 64, 128).
    - **SAGA**: codes 0--7, starting East counter-clockwise.
    - **GRASS**: codes 1--8, starting North clockwise.

    Use ``DEM.read_file(path)`` to open a raster from disk, exactly
    like ``pyramids.dataset.Dataset``.
    """

    def flow_direction_index(self, encoding: str = "esri") -> np.ndarray:
        """Convert flow-direction codes into downstream-cell indices.

        Reads the flow-direction band from the underlying raster and
        maps each D8 direction code to the row/column index of the
        downstream neighbour cell.

        Args:
            encoding: The D8 flow-direction encoding used by the
                raster.  Supported values:

                - ``"esri"`` (default) -- ArcGIS / ESRI powers-of-2
                  codes (1, 2, 4, 8, 16, 32, 64, 128).
                - ``"saga"`` -- SAGA GIS codes (0--7, starting East
                  counter-clockwise).  Produced by QGIS Processing
                  SAGA tools such as *Fill Sinks* and *Channel
                  Network*.
                - ``"grass"`` -- GRASS GIS codes (1--8, starting
                  North clockwise).  Produced by QGIS Processing
                  GRASS tools such as *r.watershed*.

        Returns:
            numpy.ndarray: A 3-D array of shape ``(rows, cols, 2)``.
                The first layer (``[:, :, 0]``) holds the row index
                and the second layer (``[:, :, 1]``) holds the column
                index of the downstream cell.  Cells with no valid
                flow direction are set to ``NaN``.

        Raises:
            ValueError: If *encoding* is not one of the supported
                names, or if the raster contains direction values
                outside the expected set for the chosen encoding.
        """
        encoding = encoding.lower()
        if encoding not in D8_ENCODINGS:
            raise ValueError(
                f"Unsupported encoding {encoding!r}. Choose from {list(D8_ENCODINGS)}"
            )
        offsets = D8_ENCODINGS[encoding]

        no_val = self.no_data_value[0]
        cols = self.columns
        rows = self.rows

        fd = self.read_array(band=0)
        fd_val = np.unique(fd[~np.isclose(fd, no_val, rtol=0.00001)])
        valid_codes = set(offsets)
        if not all(int(v) in valid_codes for v in fd_val):
            raise ValueError(
                f"Flow direction raster should contain only "
                f"{sorted(valid_codes)} for encoding {encoding!r}"
            )

        fd_cell = np.full((rows, cols, 2), np.nan)

        row_idx, col_idx = np.meshgrid(np.arange(rows), np.arange(cols), indexing="ij")
        d_row = np.full((rows, cols), np.nan)
        d_col = np.full((rows, cols), np.nan)

        for code, (dr, dc) in offsets.items():
            mask = fd == code
            d_row[mask] = dr
            d_col[mask] = dc

        valid = ~np.isnan(d_row)
        fd_cell[valid, 0] = row_idx[valid] + d_row[valid]
        fd_cell[valid, 1] = col_idx[valid] + d_col[valid]

        return fd_cell

    def flow_direction_table(self, encoding: str = "esri") -> dict:
        """Build an upstream-cell lookup table from flow directions.

        Uses ``flow_direction_index`` to determine downstream
        neighbours, then inverts the relationship so that each cell
        maps to the list of cells that flow directly into it.

        Args:
            encoding: The D8 flow-direction encoding.  See
                ``flow_direction_index`` for supported values.

        Returns:
            dict[str, list[tuple[int, int]]]: A dictionary keyed by
                ``"row,col"`` strings.  Each value is a list of
                ``(row, col)`` tuples identifying the cells whose
                flow direction points directly into the key cell.
        """
        flow_direction_index = self.flow_direction_index(encoding=encoding)

        rows = self.rows
        cols = self.columns

        cell_i = []
        cell_j = []
        celli_content = []
        cellj_content = []
        for i in range(rows):
            for j in range(cols):
                if not np.isnan(flow_direction_index[i, j, 0]):
                    # store the indexes of not empty cells and the indexes stored inside these cells
                    cell_i.append(i)
                    cell_j.append(j)
                    # store the index of the receiving cells
                    celli_content.append(flow_direction_index[i, j, 0])
                    cellj_content.append(flow_direction_index[i, j, 1])

        flow_acc_table: dict[str, list[tuple[int, int]]] = {}
        # for each cell store the directly giving cells
        for i in range(rows):
            for j in range(cols):
                if not np.isnan(flow_direction_index[i, j, 0]):
                    # get the indexes of the cell and use it as a key in a dictionary
                    name = str(i) + "," + str(j)
                    flow_acc_table[name] = []
                    for k in range(len(celli_content)):
                        # search if any cell are giving this cell
                        if i == celli_content[k] and j == cellj_content[k]:
                            flow_acc_table[name].append((cell_i[k], cell_j[k]))

        return flow_acc_table

flow_direction_index(encoding: str = 'esri') -> np.ndarray #

Convert flow-direction codes into downstream-cell indices.

Reads the flow-direction band from the underlying raster and maps each D8 direction code to the row/column index of the downstream neighbour cell.

Parameters:

Name Type Description Default
encoding str

The D8 flow-direction encoding used by the raster. Supported values:

  • "esri" (default) -- ArcGIS / ESRI powers-of-2 codes (1, 2, 4, 8, 16, 32, 64, 128).
  • "saga" -- SAGA GIS codes (0--7, starting East counter-clockwise). Produced by QGIS Processing SAGA tools such as Fill Sinks and Channel Network.
  • "grass" -- GRASS GIS codes (1--8, starting North clockwise). Produced by QGIS Processing GRASS tools such as r.watershed.
 'esri'

Returns:

Type Description
ndarray

A 3-D array of shape (rows, cols, 2). The first layer ([:, :, 0]) holds the row index and the second layer ([:, :, 1]) holds the column index of the downstream cell. Cells with no valid flow direction are set to NaN.

Raises:

Type Description
ValueError

If encoding is not one of the supported names, or if the raster contains direction values outside the expected set for the chosen encoding.
Source code in src/hapi/dem.py 
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
def flow_direction_index(self, encoding: str = "esri") -> np.ndarray:
    """Convert flow-direction codes into downstream-cell indices.

    Reads the flow-direction band from the underlying raster and
    maps each D8 direction code to the row/column index of the
    downstream neighbour cell.

    Args:
        encoding: The D8 flow-direction encoding used by the
            raster.  Supported values:

            - ``"esri"`` (default) -- ArcGIS / ESRI powers-of-2
              codes (1, 2, 4, 8, 16, 32, 64, 128).
            - ``"saga"`` -- SAGA GIS codes (0--7, starting East
              counter-clockwise).  Produced by QGIS Processing
              SAGA tools such as *Fill Sinks* and *Channel
              Network*.
            - ``"grass"`` -- GRASS GIS codes (1--8, starting
              North clockwise).  Produced by QGIS Processing
              GRASS tools such as *r.watershed*.

    Returns:
        numpy.ndarray: A 3-D array of shape ``(rows, cols, 2)``.
            The first layer (``[:, :, 0]``) holds the row index
            and the second layer (``[:, :, 1]``) holds the column
            index of the downstream cell.  Cells with no valid
            flow direction are set to ``NaN``.

    Raises:
        ValueError: If *encoding* is not one of the supported
            names, or if the raster contains direction values
            outside the expected set for the chosen encoding.
    """
    encoding = encoding.lower()
    if encoding not in D8_ENCODINGS:
        raise ValueError(
            f"Unsupported encoding {encoding!r}. Choose from {list(D8_ENCODINGS)}"
        )
    offsets = D8_ENCODINGS[encoding]

    no_val = self.no_data_value[0]
    cols = self.columns
    rows = self.rows

    fd = self.read_array(band=0)
    fd_val = np.unique(fd[~np.isclose(fd, no_val, rtol=0.00001)])
    valid_codes = set(offsets)
    if not all(int(v) in valid_codes for v in fd_val):
        raise ValueError(
            f"Flow direction raster should contain only "
            f"{sorted(valid_codes)} for encoding {encoding!r}"
        )

    fd_cell = np.full((rows, cols, 2), np.nan)

    row_idx, col_idx = np.meshgrid(np.arange(rows), np.arange(cols), indexing="ij")
    d_row = np.full((rows, cols), np.nan)
    d_col = np.full((rows, cols), np.nan)

    for code, (dr, dc) in offsets.items():
        mask = fd == code
        d_row[mask] = dr
        d_col[mask] = dc

    valid = ~np.isnan(d_row)
    fd_cell[valid, 0] = row_idx[valid] + d_row[valid]
    fd_cell[valid, 1] = col_idx[valid] + d_col[valid]

    return fd_cell

flow_direction_table(encoding: str = 'esri') -> dict #

Build an upstream-cell lookup table from flow directions.

Uses flow_direction_index to determine downstream neighbours, then inverts the relationship so that each cell maps to the list of cells that flow directly into it.

Parameters:

Name Type Description Default
encoding str

The D8 flow-direction encoding. See flow_direction_index for supported values.

 'esri'

Returns:

Type Description
dict[str, list[tuple[int, int]]]

A dictionary keyed by "row,col" strings. Each value is a list of (row, col) tuples identifying the cells whose flow direction points directly into the key cell.
Source code in src/hapi/dem.py 
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
def flow_direction_table(self, encoding: str = "esri") -> dict:
    """Build an upstream-cell lookup table from flow directions.

    Uses ``flow_direction_index`` to determine downstream
    neighbours, then inverts the relationship so that each cell
    maps to the list of cells that flow directly into it.

    Args:
        encoding: The D8 flow-direction encoding.  See
            ``flow_direction_index`` for supported values.

    Returns:
        dict[str, list[tuple[int, int]]]: A dictionary keyed by
            ``"row,col"`` strings.  Each value is a list of
            ``(row, col)`` tuples identifying the cells whose
            flow direction points directly into the key cell.
    """
    flow_direction_index = self.flow_direction_index(encoding=encoding)

    rows = self.rows
    cols = self.columns

    cell_i = []
    cell_j = []
    celli_content = []
    cellj_content = []
    for i in range(rows):
        for j in range(cols):
            if not np.isnan(flow_direction_index[i, j, 0]):
                # store the indexes of not empty cells and the indexes stored inside these cells
                cell_i.append(i)
                cell_j.append(j)
                # store the index of the receiving cells
                celli_content.append(flow_direction_index[i, j, 0])
                cellj_content.append(flow_direction_index[i, j, 1])

    flow_acc_table: dict[str, list[tuple[int, int]]] = {}
    # for each cell store the directly giving cells
    for i in range(rows):
        for j in range(cols):
            if not np.isnan(flow_direction_index[i, j, 0]):
                # get the indexes of the cell and use it as a key in a dictionary
                name = str(i) + "," + str(j)
                flow_acc_table[name] = []
                for k in range(len(celli_content)):
                    # search if any cell are giving this cell
                    if i == celli_content[k] and j == cellj_content[k]:
                        flow_acc_table[name].append((cell_i[k], cell_j[k]))

    return flow_acc_table