Please add D-Inf flow direction

[debpal]

Kind of request

Adding new functionality

Enhancement Description

pyflwdir is the best module I have worked with for watershed delineation. It is very fast and easier to handle for large raster. However, I have not found any function regarding D-inf flow direction which is also used widely. If possible, please add this functionality

Use case

No response

Additional Context

No response

[DirkEilander]

Pyflwdir only supports D8 flow at this time, and there are currently no plans to add D-inf flow. I would however be open to suggestions how to implement this from the user community.

[debpal]

@DirkEilander Thanks for the reply. D-inf flow direction is available in TauDEM (https://hydrology.usu.edu/taudem/taudem5/documentation.html). If we can go through with their documentation, it is written as D-infinity flow direction grid (flow direction grid measured in radians, counter clockwise from east). When I asked ChatGPT for help, it has suggested the following code. However, to make the code faster or make it more robust, help from an expert software developed like you is essential. I hope you could find a solution.

thanks for your response!

D-inf flow direction is well-documented in TauDEM (documentation link), which states D-infinity flow direction grid (flow direction grid measured in radians, counter clockwise from east).

Based on this, I asked ChatGPT for guidance, and it provided the following Python code as a starting point for calculating the D-inf flow direction. While it’s functional, I believe there’s room for improvement, especially in terms of performance and robustness when working with large datasets.

Since you are an expert software developer, I was wondering if you might have suggestions on how to:

a) Make the implementation faster and more efficient.
b) Integrate this functionality into your Python package, pyflwdir, for broader usability.

import numpy as np
import rasterio
from math import atan2, pi

# Load DEM using rasterio
with rasterio.open("path_to_dem.tif") as src:
    dem = src.read(1)  # Read DEM as a NumPy array
    transform = src.transform  # To handle spatial reference

# compute Slope and Direction for Each Neighbor
def compute_slope_and_direction(dem, cell_size):
    rows, cols = dem.shape
    flow_dir = np.zeros((rows, cols))  # To store flow direction
    slope = np.zeros((rows, cols))    # To store slopes

    # Offsets for 8 neighboring cells (D-infinity considers all neighbors)
    offsets = [
        (-1, 0),  # N
        (-1, 1),  # NE
        (0, 1),   # E
        (1, 1),   # SE
        (1, 0),   # S
        (1, -1),  # SW
        (0, -1),  # W
        (-1, -1), # NW
    ]

    for row in range(1, rows - 1):  # Avoid edges
        for col in range(1, cols - 1):
            max_slope = -np.inf
            flow_angle = 0

            # Elevation at the current cell
            center_elevation = dem[row, col]

            # Check all neighbors
            for dr, dc in offsets:
                neighbor_row = row + dr
                neighbor_col = col + dc
                neighbor_elevation = dem[neighbor_row, neighbor_col]

                # Distance to neighbor
                delta_x = cell_size * np.sqrt(dc**2 + dr**2)
                delta_z = center_elevation - neighbor_elevation

                # Compute slope
                current_slope = delta_z / delta_x

                # Check if this is the steepest slope
                if current_slope > max_slope:
                    max_slope = current_slope

                    # Calculate direction in radians
                    flow_angle = atan2(dr, dc)  # atan2(y, x), counterclockwise from east

            # Store results
            slope[row, col] = max_slope
            flow_dir[row, col] = flow_angle if flow_angle >= 0 else flow_angle + 2 * pi

    return slope, flow_dir

# Calculate slope and D-inf flow direction
slope, flow_dir = compute_slope_and_direction(dem, transform.a)

[DirkEilander]

This seems to return the angle to a single neighboring cell instead of the d-inf angle which is "the steepest downward slope on planar triangular facets on a block centered grid". And even if we implement a d-inf slope method in pyflwdir, this would not be compatible with the Flwdir class, which currently assumes each cell has one downstream neighbor.