priority_flow ============= .. py:module:: priority_flow .. autoapi-nested-parse:: PriorityFlow - Python toolkit for topographic processing for hydrologic models. This is a Python port of the original R package developed by Laura Condon and Reed Maxwell. Submodules ---------- .. toctree:: :maxdepth: 1 /autoapi/priority_flow/border_direction_fix/index /autoapi/priority_flow/cal_stream_order/index /autoapi/priority_flow/calc_flow/index /autoapi/priority_flow/d4_traverse/index /autoapi/priority_flow/data_loader/index /autoapi/priority_flow/define_subbasins/index /autoapi/priority_flow/define_watershed/index /autoapi/priority_flow/downstream_extract/index /autoapi/priority_flow/drainage_area/index /autoapi/priority_flow/find_orphan/index /autoapi/priority_flow/fix_drainage/index /autoapi/priority_flow/flat_fix/index /autoapi/priority_flow/get_border/index /autoapi/priority_flow/init_queue/index /autoapi/priority_flow/linear_distance/index /autoapi/priority_flow/peak_distance/index /autoapi/priority_flow/river_slope/index /autoapi/priority_flow/river_smoothing/index /autoapi/priority_flow/slope_calc_standard/index /autoapi/priority_flow/slope_calc_upwind/index /autoapi/priority_flow/stream_distance/index /autoapi/priority_flow/stream_traverse/index Functions --------- .. autoapisummary:: priority_flow.init_queue priority_flow.d4_traverse_b priority_flow.slope_calc_standard priority_flow.slope_calc_upwind priority_flow.load_dem priority_flow.load_watershed_mask priority_flow.load_river_mask priority_flow.load_all_test_data priority_flow.get_data_info priority_flow.fix_drainage priority_flow.path_extract priority_flow.calc_subbasins priority_flow.calc_stream_order priority_flow.stream_dist priority_flow.stream_traverse priority_flow.river_smooth priority_flow.riv_slope priority_flow.peak_dist priority_flow.lin_dist priority_flow.fix_flat priority_flow.get_border priority_flow.find_orphan priority_flow.drainage_area priority_flow.delin_watershed priority_flow.fix_border_dir priority_flow.calc_flow Package Contents ---------------- .. py:function:: init_queue(dem: numpy.ndarray, initmask: Optional[numpy.ndarray] = None, domainmask: Optional[numpy.ndarray] = None, border: Optional[numpy.ndarray] = None, d4: Tuple[int, int, int, int] = (1, 2, 3, 4)) -> Dict[str, numpy.ndarray] Initialize queue for topographic processing. Sets up a queue and initializes marked and step matrices for DEM processing This function is a port of the R function ``InitQueue`` from the PriorityFlow package. It sets up the initial queue of cells, along with ``marked``, ``step``, ``basin`` and ``direction`` matrices used during DEM processing. :param dem: 2D array of elevations for the domain in HydroFrame orientation (nx by ny). :type dem: np.ndarray :param initmask: Mask with the same shape as ``dem`` indicating the subset of cells to be considered for the initial queue (for example, only river cells). If ``None``, all border cells in the domain are eligible to be added to the queue. :type initmask: np.ndarray, optional :param domainmask: Mask defining the domain extent to be considered. If ``None``, the entire rectangular extent of ``dem`` is used as the domain. :type domainmask: np.ndarray, optional :param border: Optional pre-computed border mask. If ``None``, the border is derived from ``domainmask`` (cells on the outer edge of the domain and internal boundaries). :type border: np.ndarray, optional :param d4: Direction numbering system for the D4 neighbors, given as ``(down, left, up, right)``. Defaults to ``(1, 2, 3, 4)`` to match the original R implementation. :type d4: tuple of int, optional :returns: Dictionary containing: - ``"mask"``: the effective initialization mask used to define candidate outlet cells (``initmask`` after any defaults). - ``"queue"``: 2D array of outlet cells with columns ``(row, col, elevation)`` in column-major order, suitable for use with ``d4_traverse_b``. - ``"marked"``: 2D array indicating which cells were added to the initial queue (1 for queue cells, 0 otherwise). - ``"basins"``: 2D array assigning a unique basin number to each outlet cell (queue entry). - ``"direction"``: 2D array of flow directions at outlet points, pointing out of the domain, using the ``d4`` numbering scheme. :rtype: Dict[str, np.ndarray] .. py:function:: d4_traverse_b(dem: numpy.ndarray, queue: numpy.ndarray, marked: numpy.ndarray, mask: Optional[numpy.ndarray] = None, step: Optional[numpy.ndarray] = None, direction: Optional[numpy.ndarray] = None, basins: Optional[numpy.ndarray] = None, d4: Tuple[int, int, int, int] = (1, 2, 3, 4), printstep: bool = False, nchunk: int = 100, epsilon: float = 0.0, printflag: bool = False, *, n_chunk: Optional[int] = None) -> Dict[str, numpy.ndarray] Priority flow processing of D4 stream networks. This function is a direct port of the R function ``D4TraverseB`` from the PriorityFlow package. It processes all stream-network cells by walking upstream along D4 (cross-shaped) neighbors inside a mask. Where no D4 neighbors exist, it looks for D8 (diagonal) neighbors and effectively creates D4 ``bridges`` to those diagonal cells by filling intervening elevations. :param dem: 2D array of elevations for the domain in HydroFrame layout (nx, ny). :type dem: np.ndarray :param queue: Initial priority queue with shape ``(n, 3)``. Each row is ``(i, j, elevation)`` giving the row index, column index and elevation of a starting cell. Indices are 0-based. :type queue: np.ndarray :param marked: 2D array of the same shape as ``dem`` indicating which cells have already been marked (typically 1 for processed / queued cells and 0 otherwise). :type marked: np.ndarray :param mask: 2D mask with ones for cells to be processed and zeros for everything else. If ``None``, a mask of all ones with the same shape as ``dem`` is used. :type mask: np.ndarray, optional :param step: 2D array recording the step number at which each cell was processed. If ``None``, initialized to all zeros. :type step: np.ndarray, optional :param direction: 2D array of flow directions for processed cells. Direction codes follow the ``d4`` numbering scheme. If ``None``, initialized to ``np.nan`` everywhere. :type direction: np.ndarray, optional :param basins: 2D array of basin identifiers, typically created by an initialization script. When provided, any new cell added to the queue is assigned the same basin number as the cell that added it. If ``None``, initialized to all zeros. :type basins: np.ndarray, optional :param d4: Direction numbering system for the D4 neighbors, given as ``(down, left, up, right)``. Defaults to ``(1, 2, 3, 4)`` to match the original R implementation. :type d4: tuple of int, optional :param printstep: If ``True``, print the step number and current queue sizes while traversing. :type printstep: bool, optional :param nchunk: Queue-splitting parameter. The top ``nchunk`` values are kept in the primary queue for initial processing; additional cells are placed into a secondary queue and merged in chunks. :type nchunk: int, optional :param epsilon: Small value added to filled areas to avoid creating flat regions when raising elevations. Defaults to ``0.0``. :type epsilon: float, optional :param printflag: Optional flag preserved for API compatibility with the R function. Currently not used in the Python implementation. :type printflag: bool, optional :param n_chunk: Keyword alias for ``nchunk``. If provided, overrides ``nchunk``. :type n_chunk: int, optional :returns: Dictionary containing: - ``"dem"``: updated elevation array after filling. - ``"mask"``: mask used during processing. - ``"marked"``: updated marked array. - ``"step"``: step number at which each cell was processed. - ``"direction"``: flow direction codes for each cell. - ``"basins"``: basin identifier assigned to each cell. :rtype: Dict[str, np.ndarray] .. rubric:: Notes All arrays are 2D with shape ``(nx, ny)`` and use 0-based indexing. The traversal walks upstream, so direction codes point from each cell toward its downstream neighbor following the chosen ``d4`` numbering. .. py:function:: slope_calc_standard(dem: numpy.ndarray, direction: numpy.ndarray, dx: float, dy: float, mask: Optional[numpy.ndarray] = None, d4: Tuple[int, int, int, int] = (1, 2, 3, 4), minslope: float = 0, maxslope: float = -1, secondary_th: float = -1, printflag: bool = False) -> Dict[str, numpy.ndarray] Calculate the slopes from a DEM. This function will calculate slopes using standard or upwinding options and apply a range of smoothing options. :param dem: 2D array of elevations in **HydroFrame** layout. :type dem: np.ndarray :param direction: 2D array of D4 flow directions for each cell, using the convention encoded in ``d4``. :type direction: np.ndarray :param dx: Lateral grid cell resolutions in the x and y directions, respectively. :type dx: float :param dy: Lateral grid cell resolutions in the x and y directions, respectively. :type dy: float :param mask: 2D mask defining the domain extent to be considered. Cells with value 0 are excluded from slope calculations. If ``None``, the full rectangular domain of ``dem`` is used. :type mask: np.ndarray, optional :param d4: Direction numbering system for the D4 neighbors, given as ``(down, left, up, right)``. Defaults to ``(1, 2, 3, 4)`` to match the original R implementation. :type d4: tuple of int, optional :param minslope: Minimum absolute slope to enforce on primary-direction slopes. If ``minslope >= 0``, primary slopes with magnitude less than this value are adjusted up to ``minslope`` (with sign preserved). Defaults to ``0``. :type minslope: float, optional :param maxslope: Maximum absolute slope to enforce. If ``maxslope >= 0``, slopes with magnitude greater than this value are limited to ``±maxslope``. If set to ``-1``, no maximum is applied. :type maxslope: float, optional :param secondary_th: Secondary threshold - maximum ratio of |secondary|/|primary| to be enforced. NOTE - this scaling occurs after any max threholds are applied. Currently this is only working for two options: (1) If this is set to -1 no scaling will be applied, (2) If this is set to zero all seconeary slopes will be zero :type secondary_th: float, optional :param printflag: If True, print progress information and details about slope limiting and secondary-slope handling. :type printflag: bool, optional .. rubric:: Notes **River methods**: - ``0``: default value, no special treatment for river cells. - ``1``: scale secondary slopes along the river (requires a river mask and a river-specific secondary threshold). - ``2``: apply watershed mean slope to each river reach (requires river mask and subbasins). - ``3``: apply the stream mean slope to each reach (requires river mask and subbasins). NOTE: The river mask can differ from the rivers used to create the subbasins (for example, using a higher threshold to define subbasins and a lower threshold to define river cells). :returns: Dictionary containing: - ``"slopex"``: 2D array of adjusted x-direction face-centered slopes. - ``"slopey"``: 2D array of adjusted y-direction face-centered slopes. - ``"direction"``: the (possibly renumbered) flow-direction array used when constructing primary/secondary masks. - ``"Sinks"``: 1D array of linear indices (flattened) marking locations where slope signs were flipped to remove sinks. :rtype: Dict[str, np.ndarray] .. py:function:: slope_calc_upwind(dem: numpy.ndarray, direction: numpy.ndarray, dx: float, dy: float, mask: numpy.ndarray | None = None, borders: numpy.ndarray | None = None, borderdir: int = 1, d4: tuple[int, int, int, int] = (1, 2, 3, 4), minslope: float = 1e-05, maxslope: float = -1, secondary_th: float = -1, river_method: int = 0, river_secondary_th: float = 0, rivermask: numpy.ndarray | None = None, subbasins: numpy.ndarray | None = None, printflag: bool = False, upflag: bool = True) -> dict[str, numpy.ndarray] Calculate slopes from a DEM. This function will calculate slopes using standard or upwinding options and apply a range of smoothing options. :param dem: Digital elevation model (nx x ny) in HydroFrame layout. :type dem: np.ndarray :param direction: Flow direction (1=down, 2=left, 3=up, 4=right). May be modified in place for borders. :type direction: np.ndarray :param dx: Lateral grid cell resolution. :type dx: float :param dy: Lateral grid cell resolution. :type dy: float :param mask: Domain mask (1=active, 0=inactive). If None, full rectangular domain is used. :type mask: np.ndarray or None :param borders: Matrix with 1 for border cells default out, 2 default in, 0 for non-border. Built from mask if None. :type borders: np.ndarray or None :param borderdir: Default for border cells: 1=point out, 2=point in. :type borderdir: int :param d4: Direction codes (down, left, up, right). Default (1, 2, 3, 4). :type d4: tuple :param minslope: Minimum absolute slope for flat cells. Use -1 to disable minimum enforcement. :type minslope: float :param maxslope: Maximum absolute slope. Use -1 to disable. :type maxslope: float :param secondary_th: Max ratio |secondary|/|primary|. Use -1 to disable. :type secondary_th: float :param river_method: 0=none, 1=scale secondary along river, 2=watershed mean slope, 3=stream mean slope. :type river_method: int :param river_secondary_th: Secondary ratio for river cells when river_method 1-3. :type river_secondary_th: float :param rivermask: River mask (1=river). Required for river_method 1, 2, 3. :type rivermask: np.ndarray or None :param subbasins: Subbasin ids. Required for river_method 2, 3. :type subbasins: np.ndarray or None :param printflag: Print progress. :type printflag: bool :param upflag: If True, downwind slopes for OverlandFlow; if False, standard [i+1]-[i] (OverlandKin). :type upflag: bool .. rubric:: Notes **River methods**: - ``0``: default value, no special treatment for river cells. - ``1``: scale secondary slopes along the river (requires a river mask and a river-specific secondary threshold). - ``2``: apply watershed mean slope to each river reach (requires river mask and subbasins). - ``3``: apply the stream mean slope to each reach (requires river mask and subbasins). NOTE: The river mask can differ from the rivers used to create the subbasins (for example, using a higher threshold to define subbasins and a lower threshold to define river cells). :returns: "slopex", "slopey", "direction" (nx x ny arrays). :rtype: dict .. py:function:: load_dem() Load the Digital Elevation Model (DEM) test data. This is a small elevation dataset (215km by 172km at 1km spatial resolution) converted from the R package's TestDomain_Inputs. :returns: A 2D array of elevation values with shape (215, 172) representing the domain dimensions (nrow=215, ncol=172). :rtype: numpy.ndarray .. rubric:: Examples >>> import priority_flow.data_loader as dl >>> dem = dl.load_dem() >>> print(f"DEM shape: {dem.shape}") >>> print(f"Elevation range: {dem.min():.2f} to {dem.max():.2f}") .. py:function:: load_watershed_mask() Load the watershed mask test data. A mask showing the watershed drainage area for the test domain. :returns: A 2D array of 0's and 1's showing the watershed extent (1=inside the watershed, 0=outside the watershed) with shape (215, 172). :rtype: numpy.ndarray .. rubric:: Examples >>> import priority_flow.data_loader as dl >>> mask = dl.load_watershed_mask() >>> print(f"Watershed mask shape: {mask.shape}") >>> print(f"Watershed cells: {np.sum(mask)}") .. py:function:: load_river_mask() Load the river mask test data. A mask showing an example river network for the test domain. :returns: A 2D array of 0's and 1's showing the location of river cells (1=river, 0=non-river) with shape (215, 172). :rtype: numpy.ndarray .. rubric:: Examples >>> import priority_flow.data_loader as dl >>> river_mask = dl.load_river_mask() >>> print(f"River mask shape: {river_mask.shape}") >>> print(f"River cells: {np.sum(river_mask)}") .. py:function:: load_all_test_data() Load all test data files at once. :returns: A dictionary containing all test data arrays with keys: - 'dem': Digital Elevation Model - 'watershed_mask': Watershed drainage area mask - 'river_mask': River network mask :rtype: dict .. rubric:: Examples >>> import priority_flow.data_loader as dl >>> data = dl.load_all_test_data() >>> print(f"Available data: {list(data.keys())}") >>> print(f"DEM shape: {data['dem'].shape}") .. py:function:: get_data_info() Get information about the available test data files. :returns: A dictionary containing metadata about each data file. :rtype: dict .. py:function:: fix_drainage(dem: numpy.ndarray, direction: numpy.ndarray, mask: numpy.ndarray, bank_epsilon: float, startpoint: Union[Tuple[int, int], list, numpy.ndarray], d4: Tuple[int, int, int, int] = (1, 2, 3, 4)) -> Dict[str, numpy.ndarray] Walk upstream from a point ensuring DEM is increasing by a minum epsilon. This walks upstream from a point following a flow direction file and checks that the elevation upstream is greater than or equal to the elvation at the point + epsilon Once the function reaches a point where upstream cells pass it stops NOTE: this is only processing the immediate neigborhood and does not recurse over the entire domain, For example if you did the entire overal prioirty flow processing with an espilon of 0 then ran this functon starting at a river point with an epsilon of 0.1 this function will traverse upstream from the river bottom checking that every cell conneting to the river cell is higher by at least this ammout but once it reaches a point where every connected cell passes this test it will stop. Therefore there could still be locations higher up on the hillslope with the original epsilon of zero still. This is on purpose and this script is not intended to gobally ensure a given epsilon. :param dem: 2D array of elevations for the domain (nx, ny). :type dem: np.ndarray :param direction: 2D array of D4 flow directions for each cell. Direction codes follow the ``d4`` numbering scheme. :type direction: np.ndarray :param mask: 2D mask with ones for cells that are allowed to be processed and zeros elsewhere. Only neighbors within the mask can be adjusted. :type mask: np.ndarray :param bank_epsilon: Minimum elevation difference that upstream cells must have relative to the current cell. If an upstream neighbor is lower than ``current + bank_epsilon``, its elevation is raised to ``current + bank_epsilon``. :type bank_epsilon: float :param startpoint: Starting grid cell given as ``(row, col)`` indices (0-based) from which to begin walking upstream. :type startpoint: tuple[int, int] or list or np.ndarray :param d4: Direction numbering system for the D4 neighbors, given as ``(down, left, up, right)``. Defaults to ``(1, 2, 3, 4)`` to match the original R implementation. :type d4: tuple of int, optional :returns: Dictionary containing: - ``"dem.adj"``: 2D elevation array after local adjustments. - ``"processed"``: 2D array marking cells that were adjusted (1) or visited during processing, and 0 elsewhere. :rtype: Dict[str, np.ndarray] .. py:function:: path_extract(input: numpy.ndarray, direction: numpy.ndarray, mask: Optional[numpy.ndarray] = None, startpoint: Optional[Union[Tuple[int, int], list, numpy.ndarray]] = None, d4: Tuple[int, int, int, int] = (1, 2, 3, 4)) -> Dict[str, numpy.ndarray] Walk downstream from a starting point and extract values from a matrix. This is a port of the R function ``PathExtract`` from the PriorityFlow package. Starting from a given grid cell, it follows the D4 flow-direction grid ``direction`` downstream, collecting values from ``input`` until the path leaves the domain or exits the processing mask. :param input: 2D array of values from which to extract the stream-path sequence (for example, elevation, drainage area, or any other field defined on the same grid as ``direction``). :type input: np.ndarray :param direction: 2D array of D4 flow directions for each cell. Direction codes follow the ``d4`` numbering scheme. :type direction: np.ndarray :param mask: 2D mask with ones for cells that are allowed to be part of the path and zeros elsewhere. If ``None``, a mask of all ones with the same shape as ``direction`` is used. :type mask: np.ndarray, optional :param startpoint: Starting grid cell given as ``(row, col)`` indices (0-based). If an array-like object is provided, only the first two elements are used. :type startpoint: tuple[int, int] or list or np.ndarray, optional :param d4: Direction numbering system for the D4 neighbors, given as ``(down, left, up, right)``. Defaults to ``(1, 2, 3, 4)`` to match the original R implementation. :type d4: tuple of int, optional :returns: Dictionary containing: - ``"data"``: 1D array of values extracted from ``input`` along the downstream path, in order. - ``"path_mask"``: 2D array with non-zero entries indicating the step number at each cell visited along the path. - ``"path_list"``: 2D array of ``(row, col)`` indices for the cells visited along the path, in traversal order. :rtype: Dict[str, np.ndarray] .. py:function:: calc_subbasins(direction: numpy.ndarray, area: numpy.ndarray, mask: Optional[numpy.ndarray] = None, d4: Tuple[int, int, int, int] = (1, 2, 3, 4), riv_th: int = 50, printflag: bool = False, merge_th: int = 0) -> Dict[str, numpy.ndarray] Define subbasins and stream segments from flow direction and drainage area. This function is a port of the R function ``CalcSubbasins`` from the PriorityFlow package. It divides the domain into subbasins and individual stream segments using a flow-direction grid and a drainage-area grid. :param direction: 2D array of D4 flow directions for each cell. Direction codes follow the ``d4`` numbering scheme. :type direction: np.ndarray :param area: 2D array of drainage area values for every cell (typically in number of contributing cells or an equivalent measure). :type area: np.ndarray :param mask: 2D mask with ones for cells to be processed and zeros for everything else. If ``None``, a mask of all ones with the same shape as ``direction`` is used. :type mask: np.ndarray, optional :param d4: Direction numbering system for the D4 neighbors, given as ``(down, left, up, right)``. Defaults to ``(1, 2, 3, 4)`` to match the original R implementation. :type d4: tuple of int, optional :param riv_th: Drainage-area threshold used to classify river cells. Cells with ``area >= riv_th`` are treated as part of the river network. Defaults to ``50``. :type riv_th: int, optional :param printflag: If ``True``, print basic progress information while growing subbasins upstream from the river network. :type printflag: bool, optional :param merge_th: After all subbasins have been defined, subbasins with total area smaller than ``merge_th`` are merged into their downstream neighbors. A value of ``0`` (the default) disables merging. :type merge_th: int, optional :returns: Dictionary containing: - ``"segments"``: 2D array with a unique ID for each stream segment. - ``"subbasins"``: 2D array with a unique ID for each subbasin (segment plus its contributing area). - ``"RiverMask"``: 2D mask of river cells derived from ``area`` and ``riv_th``. - ``"summary"``: 2D array where each row summarizes a subbasin/segment (segment ID, representative coordinates, downstream segment ID and cell count). Coordinates are stored using 0-based indexing. :rtype: Dict[str, np.ndarray] .. py:function:: calc_stream_order(basin_id: numpy.ndarray, ds_id: numpy.ndarray, segments: Optional[numpy.ndarray] = None) -> Dict[str, Union[numpy.ndarray, numpy.ndarray]] Calculate Strahler stream orders using stream segments delineated using the CalcSubbasins function. :param basin_id: An array of basin IDs (can be obtained from the CalcSubbasins function Summary output column 1 "Basin_ID") :type basin_id: np.ndarray :param ds_id: Downstream basin IDs, the downstream ID of each basin, should be corresponding to the basin ID (can be obtained from the CalcSubbasins function Summary output Column 6 "Downstream_Basin_ID") :type ds_id: np.ndarray :param segments: A channel mask with segments assigned with basin IDs (can use the CalcSubbasins 'segments' output for this) :type segments: np.ndarray, optional :returns: A dictionary containing: - 'summary': A summary table with a row for every basin with three columns "Basin_ID", "Downstream_ID", "StreamOrder_number" - 'order_mask': (optional) Only available if segments is an input, this will output a spatial raster with the channel orders :rtype: Dict[str, np.ndarray] .. rubric:: Notes This function implements the Strahler stream ordering system where: - First-order streams are those without any tributaries - When two streams of the same order join, the resulting stream is of the next higher order - When streams of different orders join, the resulting stream maintains the higher order .. py:function:: stream_dist(direction: numpy.ndarray, streammask: numpy.ndarray, domainmask: Optional[numpy.ndarray] = None, d4: Tuple[int, int, int, int] = (1, 2, 3, 4)) -> Dict[str, numpy.ndarray] Find stream distance using flow-direction traversal upstream from streams. .. py:function:: stream_traverse(dem: numpy.ndarray, mask: numpy.ndarray, queue: numpy.ndarray, marked: numpy.ndarray, step: Optional[numpy.ndarray] = None, direction: Optional[numpy.ndarray] = None, basins: Optional[numpy.ndarray] = None, d4: Tuple[int, int, int, int] = (1, 2, 3, 4), printstep: bool = False, epsilon: float = 0.0) -> Dict[str, numpy.ndarray] DEM processing of stream networks. Function to process stream networks walking upstream on d4 neigbors in a river mask. Where no D4 neigbhors exist it looks for d8 neigbors and creates d4 bridges to these diagonal cells. :param dem: Digital Elevation Model matrix (nx x ny). :type dem: np.ndarray :param mask: Mask with zeros for non-river cells and 1 for river cells. :type mask: np.ndarray :param queue: Priority queue to start from; columns: x (row), y (col), elevation. :type queue: np.ndarray :param marked: Matrix of which cells have been marked already. :type marked: np.ndarray :param step: Matrix of the step number for processed cells. Defaults to zeros. :type step: np.ndarray, optional :param direction: Matrix of flow directions for processed cells. Defaults to NaN. :type direction: np.ndarray, optional :param basins: Matrix of basin numbers (e.g., from InitQueue). If provided, every newly added cell inherits the basin of the cell that adds it. :type basins: np.ndarray, optional :param d4: D4 direction numbering (down, left, up, right). Default (1,2,3,4). :type d4: tuple of int, optional :param printstep: If True, print step number and queue size. :type printstep: bool, optional :param epsilon: Amount to add to filled areas to avoid creating flats. :type epsilon: float, optional :returns: Dictionary with keys: "dem", "mask", "marked", "step", "direction", and "basins". :rtype: Dict[str, np.ndarray] .. py:function:: river_smooth(dem: numpy.ndarray, direction: numpy.ndarray, mask: Optional[numpy.ndarray] = None, river_summary: Optional[numpy.ndarray] = None, river_segments: Optional[numpy.ndarray] = None, bank_epsilon: float = 0.01, river_epsilon: float = 0.0, d4: Tuple[int, int, int, int] = (1, 2, 3, 4), printflag: bool = False) -> Dict[str, numpy.ndarray] Smooth a DEM along a predefined stream network and enforce bank slopes. This function will smooth a DEM along a stream network. It requires pre-defined stream segments and subbasins which are obtained using the CalcSubbasins function. :param dem: 2D array of elevations for the domain (nx, ny). :type dem: np.ndarray :param direction: 2D array of D4 flow directions for each cell, using the convention encoded in ``d4``. :type direction: np.ndarray :param mask: 2D mask defining the domain extent to be considered. Cells with value 0 are excluded from processing. If ``None``, the entire DEM is used. :type mask: np.ndarray, optional :param river_summary: Summary table of stream segments with one row per segment and seven columns: 1. Subbasin (segment) ID number. 2. Row index of the upstream end of the segment. 3. Column index of the upstream end of the segment. 4. Row index of the downstream end of the segment. 5. Column index of the downstream end of the segment. 6. Subbasin ID of the downstream basin (``-1`` indicates that the segment drains out of the domain). 7. Drainage area of the subbasin. Indices are 0-based (row, col) for array indexing in this Python implementation. :type river_summary: np.ndarray, optional :param river_segments: 2D array with the same shape as ``dem`` indicating the subbasin (segment) ID for each cell on the river network. Cells not on the river network should be zero. :type river_segments: np.ndarray, optional :param bank_epsilon: Minimum elevation difference between river cells and adjacent hillslope cells when walking up the banks from the river network. Passed through to ``fix_drainage``. :type bank_epsilon: float, optional :param river_epsilon: Minimum elevation difference per cell along the river segments. Used to enforce a monotonic drop in elevation from the upstream to downstream end of each segment. :type river_epsilon: float, optional :param d4: Direction numbering system for the D4 neighbors, given as ``(down, left, up, right)``. Defaults to ``(1, 2, 3, 4)`` to match the original R implementation. :type d4: tuple of int, optional :param printflag: If ``True``, print progress and diagnostic information while smoothing each river segment. :type printflag: bool, optional :returns: Dictionary containing: - ``"dem.adj"``: 2D array of DEM elevations after river smoothing and bank adjustments. - ``"processed"``: 2D mask indicating cells that were processed by this routine (1 = processed, 0 = not processed), which should match the river network mask when all segments are successfully handled. - ``"summary"``: 2D array summarizing each river segment with one row per segment and the following columns: 1. River segment ID number. 2. Row index of the top of the segment. 3. Column index of the top of the segment. 4. Row index of the bottom of the segment. 5. Column index of the bottom of the segment. 6. Length of the segment (number of cells). 7. Elevation at the top of the segment. 8. Elevation at the bottom of the segment. 9. Delta applied along the segment (i.e. (top-bottom)/length). :rtype: Dict[str, np.ndarray] .. py:function:: riv_slope(direction: numpy.ndarray, slopex: numpy.ndarray, slopey: numpy.ndarray, minslope: float, river_mask: Optional[numpy.ndarray] = None, remove_sec: bool = False) -> Dict[str, numpy.ndarray] Apply minimum slope and secondary scaling to river cells. A funtion that will apply a minimum slope threshold for the primary flow direciton along river cells. This function can also limit secondary slopes out of river cells to zero. :param direction: 2D array of D4 flow directions for each cell, using the convention ``1 = down``, ``2 = left``, ``3 = up``, ``4 = right``. :type direction: np.ndarray :param slopex: 2D array of face-centered slopes in the x-direction. These should be the same slopes produced by ``slope_calc_standard``. :type slopex: np.ndarray :param slopey: 2D array of face-centered slopes in the y-direction, compatible with ``slope_calc_standard``. :type slopey: np.ndarray :param minslope: Threshold for slope adjustment. For each river cell, the absolute value of the primary-direction slope is adjusted so that ``abs(slope) >= minslope``. :type minslope: float :param river_mask: 2D mask indicating which cells are treated as river cells (typically 1 for river, 0 for non-river). Only these cells are subject to minimum-slope and secondary-slope adjustments. :type river_mask: np.ndarray, optional :param remove_sec: If ``True``, any secondary outflows from river cells (i.e. positive slopes in non-primary directions) are set to zero. Defaults to ``False``. :type remove_sec: bool, optional :returns: Dictionary containing: - ``"slopex"``: adjusted x-direction slope field. - ``"slopey"``: adjusted y-direction slope field. - ``"adj_mask"``: mask of cells that were adjusted (non-zero where primary or secondary slopes were modified). - ``"SlopeOutlet"``: original primary outlet slope at each river cell before adjustment. - ``"SlopeOutletNew"``: primary outlet slope at each river cell after adjustment. :rtype: Dict[str, np.ndarray] .. py:function:: peak_dist(direction: numpy.ndarray, mask: Optional[numpy.ndarray] = None, d4: Tuple[int, int, int, int] = (1, 2, 3, 4), printflag: bool = False) -> Dict[str, numpy.ndarray] Calculate distance to topographic peaks. Calculates the maximum and minimum distance from a headwater cell to topographic peaks using a topological sorting approach. This function identifies headwater cells and propagates distance information downstream following flow directions. NOTE: This function has not yet been tested and is not fully implemented. :param direction: Flow direction matrix :type direction: np.ndarray :param mask: Processing mask. Defaults to processing everything :type mask: np.ndarray, optional :param d4: Directional numbering system: down, left, top, right. Defaults to (1, 2, 3, 4) :type d4: Tuple[int, int, int, int], optional :param printflag: Whether to print progress information. Defaults to False :type printflag: bool, optional :returns: A dictionary containing: - 'mindist': Matrix containing the minimum distance from headwaters to each cell - 'maxdist': Matrix containing the maximum distance from headwaters to each cell :rtype: Dict[str, np.ndarray] .. rubric:: Notes This function implements a distance calculation algorithm that: 1. Identifies headwater cells (cells with no upstream neighbors) 2. Uses topological sorting to process cells in correct order 3. Propagates distance information downstream following flow directions 4. Tracks both minimum and maximum distances from headwaters The algorithm handles: - Border detection and processing - D4 neighbor connectivity - Upstream neighbor counting - Topological sorting for correct distance propagation - Masked area handling .. py:function:: lin_dist(target_points: numpy.ndarray, mask: Optional[numpy.ndarray] = None, cell_size: float = 1.0, max_dist: Optional[float] = None, printflag: bool = False) -> numpy.ndarray Calculate the minimum linear distance from a feature set. Calculates the minimum linear distance between every point in a 2D array and a mask of target points. Note that this function assumes dx=dy, i.e., square grid cells. If this is not the case, it will not work correctly. NOTE: This function has not yet been tested and is not fully implemented. :param target_points: Matrix with a value of 1 for all cells that you would like to calculate the distance to and 0 for all other cells :type target_points: np.ndarray :param mask: Processing mask. Defaults to processing everything :type mask: np.ndarray, optional :param cell_size: The size of a grid cell. Note that this function assumes dx=dy, i.e., square grid cells. Defaults to 1.0 :type cell_size: float, optional :param max_dist: Maximum distance to check. All cells with no distance found <= max_dist will be assigned max_dist. If None, checks out to the total size of the domain :type max_dist: float, optional :param printflag: Whether to print progress information. Defaults to False :type printflag: bool, optional :returns: Matrix containing the minimum distance from each cell to the nearest target point. Cells outside the mask are assigned NaN. :rtype: np.ndarray .. rubric:: Notes This function implements a step-by-step distance checking algorithm that: 1. Creates an ordered list of distance steps and rotations to check 2. Iteratively checks each distance incrementally 3. Uses array shifting to efficiently count target points at each distance 4. Assigns distances to cells as they are found 5. Continues until all cells have been assigned distances or max_dist is reached The algorithm handles: - Square grid cells (dx = dy) - Configurable maximum distance limits - Efficient array operations for distance checking - Mask-based domain processing - Progressive distance assignment The function assumes square grid cells and uses Euclidean distance calculations. For non-square grids, results may be incorrect. .. py:function:: fix_flat(direction: numpy.ndarray, slopex: numpy.ndarray, slopey: numpy.ndarray, adj_th: float, adj_ratio: float = -1.0, mask: Optional[numpy.ndarray] = None) -> Dict[str, numpy.ndarray] Identify and fix stagnation points. A function that finds cells where the total output slope is much less than the input slopes and adjusts the outlet slope to address this. This is useful for fixing flat areas and stagnation points in flow calculations. :param direction: Nx by Ny matrix of flow directions following the convention (1=down, 2=left, 3=up, 4=right) :type direction: np.ndarray :param slopex: Nx by Ny matrix of slopes in the x direction (should be face centered slopes as calculated with SlopeCalcStan) :type slopex: np.ndarray :param slopey: Nx by Ny matrix of slopes in the y direction (should be face centered slopes as calculated with SlopeCalcStan) :type slopey: np.ndarray :param adj_th: Threshold for slope adjustment. If the total slopes out divided by the total slopes into a given cell is less than adj_th, the outlet will be scaled by adj_ratio. :type adj_th: float :param adj_ratio: Scaler value for slope adjustment. New outlet slopes will be set to initial outlet slope times adj_ratio. If no adj_ratio is provided, then the outlet slope will be set to the total inlets slopes times the adjustment threshold. Defaults to -1.0. :type adj_ratio: float, optional :param mask: Nx by Ny matrix indicating the active domain, 1=active, 0=inactive. If no mask is provided, the function will be applied over the entire input matrices. :type mask: np.ndarray, optional :returns: A dictionary containing: - 'outslope': Nx by Ny matrix with the sum of the slopes pointed out of every cell (sum of slope magnitudes, all positive) - 'inslope': Nx by Ny matrix with the sum of the slopes pointed into every cell (sum of slope magnitudes, all positive) - 'OutIn_ratio': Nx by Ny matrix with the ratio of outslope to inslope (cells outside mask or with 0 inslope are assigned -1) - 'slopex': The adjusted Nx by Ny matrix of slopex values - 'slopey': The adjusted Nx by Ny matrix of slopey values - 'adj_mask': Nx by Ny matrix indicating cells that were adjusted (1=adjusted, 0=not adjusted) - 'SlopeOutlet': Nx by Ny matrix of the outlet slope for every grid cell - 'SlopeOutletNew': Nx by Ny matrix of the outlet slope for every grid cell after processing :rtype: Dict[str, np.ndarray] .. rubric:: Notes This function implements a stagnation point detection and correction algorithm that: 1. Calculates total input and output slopes for each cell 2. Identifies cells with low output/input slope ratios 3. Adjusts outlet slopes to meet minimum threshold requirements 4. Provides comprehensive slope analysis and adjustment tracking The algorithm handles: - Face-centered slope calculations - Direction-based slope adjustments - Threshold-based correction criteria - Comprehensive output reporting - Mask-based domain processing .. py:function:: get_border(mask_mat: numpy.ndarray) -> numpy.ndarray Identify border cells. Function that reads in a mask and returns a mask of the border cells for the irregular boundary. Border cells are identified as cells that have fewer than 4 valid neighbors within the domain. :param mask_mat: Matrix mask with values of 0 for cells outside the domain and 1 for cells inside the domain :type mask_mat: np.ndarray :returns: Matrix mask where 1 indicates border cells and 0 indicates non-border cells :rtype: np.ndarray .. rubric:: Notes This function implements a border detection algorithm that: 1. Identifies cells at the domain boundary 2. Counts valid neighbors for each cell 3. Classifies cells with fewer than 4 neighbors as border cells 4. Handles irregular domain boundaries The algorithm uses D4 connectivity (4-directional) to identify border cells: - Cells with 4 valid neighbors are interior cells (not borders) - Cells with fewer than 4 valid neighbors are border cells - Only cells within the domain mask are considered Border cells are typically: - Domain edge cells - Cells adjacent to masked areas - Cells in irregular boundary regions - Cells with incomplete neighbor connectivity .. py:function:: find_orphan(dem: numpy.ndarray, mask: numpy.ndarray, marked: numpy.ndarray) -> Dict[str, Union[int, numpy.ndarray, None]] Find orphan branches in river network processing. Function to look for unprocessed river cells that have D8 neighbors on the river network or on the boundary. :param dem: Digital elevation model matrix :type dem: np.ndarray :param mask: River network mask (1 for river cells, 0 for non-river cells) :type mask: np.ndarray :param marked: Matrix of grid cells that have been processed (1 for processed, 0 for unprocessed) :type marked: np.ndarray :returns: A dictionary containing: - 'norphan': Number of orphaned branches found - 'queue': Queue of marked neighbors to be processed (numpy array with columns: x, y, elevation) :rtype: Dict[str, Union[int, np.ndarray, None]] .. rubric:: Notes This function implements an orphan detection algorithm that: 1. Identifies unprocessed river cells (masked but not marked) 2. Checks D8 connectivity to find marked neighbors 3. Counts marked neighbors for each orphan cell 4. Creates a processing queue from marked neighbors of orphan cells The algorithm uses D8 connectivity (8-directional) to ensure complete neighbor checking and proper orphan identification. .. py:function:: drainage_area(direction: numpy.ndarray, mask: Optional[numpy.ndarray] = None, d4: Tuple[int, int, int, int] = (1, 2, 3, 4), printflag: bool = False) -> numpy.ndarray Calculate drainage area (number of upstream cells) for each grid cell given a flow direction file. :param direction: 2D array of D4 flow directions for each cell. Direction codes follow the ``d4`` numbering scheme. :type direction: np.ndarray :param mask: 2D mask with ones for cells to be processed and zeros for everything else. Cells outside the mask are excluded from the drainage-area computation. If ``None``, a mask of all ones with the same shape as ``direction`` is used. :type mask: np.ndarray, optional :param d4: Direction numbering system for the D4 neighbors, given as ``(down, left, up, right)``. Defaults to ``(1, 2, 3, 4)`` to match the original R implementation. :type d4: tuple of int, optional :param printflag: If ``True``, print basic progress information as the algorithm walks downstream through the queue of cells. :type printflag: bool, optional :returns: 2D array of drainage areas with the same shape as ``direction``. Each entry gives the number of cells (including itself) that drain to that cell, multiplied by ``mask`` so cells outside the mask have value 0. :rtype: np.ndarray .. py:function:: delin_watershed(outlets: numpy.ndarray, direction: numpy.ndarray, d4: Tuple[int, int, int, int] = (1, 2, 3, 4), printflag: bool = False) -> Dict[str, numpy.ndarray] Function to define the watershed for a point or set of outlet points based on the flow direction file. NOTE: This function is not yet tested. :param outlets: x,y coordinates of the outlet points or points to mask upstream areas for. This should be a 2D array with a separate row for each point. :type outlets: np.ndarray :param direction: Flow direction matrix :type direction: np.ndarray :param d4: Directional numbering system: down, left, top, right. Defaults to (1, 2, 3, 4) :type d4: Tuple[int, int, int, int], optional :param printflag: Optional flag to print out the number of cells in the queue during iterations. Defaults to False :type printflag: bool, optional :returns: A dictionary containing: - 'watershed': Binary mask of the watershed area (1 for watershed cells, 0 otherwise) - 'xrange': Range of x coordinates covered by the watershed as numpy array [min_x, max_x] - 'yrange': Range of y coordinates covered by the watershed as numpy array [min_y, max_y] :rtype: Dict[str, np.ndarray] .. rubric:: Notes This function implements an upstream traversal algorithm that: 1. Starts from specified outlet points 2. Traverses upstream following flow directions 3. Marks all cells that contribute flow to the outlet points 4. Returns the complete watershed boundary and extent The algorithm uses a queue-based approach to efficiently process large watersheds. .. py:function:: fix_border_dir(direction: numpy.ndarray, dem: numpy.ndarray, d4: Tuple[int, int, int, int] = (1, 2, 3, 4)) -> numpy.ndarray Fix the directions of border cells. Originally all border cells are initialized to point out of the domain. This function checks the slopes and areas to determine if they need to swap direction to point into the domain instead. NOTE: This function is not yet tested. :param direction: Flow direction matrix [nx, ny] - will be modified in place :type direction: np.ndarray :param dem: Digital Elevation Model matrix [nx, ny] :type dem: np.ndarray :param d4: Direction numbering system for down, left, top, right. Defaults to (1, 2, 3, 4). :type d4: tuple, optional :returns: The modified direction array :rtype: np.ndarray .. rubric:: Notes This is a Python port of the R function FixBorderDir from PriorityFlow. The function analyzes slopes at domain boundaries to determine whether border cells should drain inward or outward. WARNING: This function modifies the input direction array in place for backwards compatibility with the R version. Direction mapping: - d4[0] = 1: Down - d4[1] = 2: Left - d4[2] = 3: Up - d4[3] = 4: Right .. py:function:: calc_flow(file_path: str, run_name: str, file_nums: List[int], slopex_file: str, slopey_file: str, overland_type: str, mannings: Union[float, numpy.ndarray], epsilon: float = 1e-05, dx: float = 1.0, dy: float = 1.0, mask: Optional[numpy.ndarray] = None) -> None Calculate overland flow from ParFlow pressure file outputs. This function will write three pfb outputs for every timestep: (1) outflow - volumetric outflow from each grid cell [l3/t] (2) q_east - volumetric flow across the east face of each cell (Note: this will have dimensions nx+1, ny to account for the western edge of the domain) (3) q_north - volumetric flow across the north face of each cell. (Note this will have dimensions nx, ny+1 to account for the southern edge of the domain) :param file_path: Directory path where pressure and slope files are located, also where flow files will be written :type file_path: str :param run_name: ParFlow run name used to read pressure files (standard naming: runname.out.press.00000.pfb) :type run_name: str :param file_nums: List of file numbers to be processed :type file_nums: List[int] :param slopex_file: Name of the x slopes file (must be a pfb file) :type slopex_file: str :param slopey_file: Name of the y slopes file (must be a pfb file) :type slopey_file: str :param overland_type: Type of overland flow used for ParFlow simulation. Choices: "OverlandFlow" or "OverlandKinematic" :type overland_type: str :param mannings: Manning roughness coefficient, can be either a value or a matrix :type mannings: Union[float, np.ndarray] :param epsilon: Epsilon used only in OverlandKinematic formulation. Set with Solver.OverlandKinematic.Epsilon key in ParFlow. Defaults to 1e-5. :type epsilon: float, optional :param dx: Grid size in x direction [l]. Defaults to 1.0. :type dx: float, optional :param dy: Grid size in y direction [l]. Defaults to 1.0. :type dy: float, optional :param mask: Mask with ones for cells to be processed and zeros for everything else. Defaults to mask of all 1's. Should be [nx, ny] matrix where [0,0] is lower left corner. :type mask: np.ndarray, optional .. rubric:: Notes This is a Python port of the R function CalcFlow from PriorityFlow. The function calculates overland flow using either OverlandFlow or OverlandKinematic formulations. File naming convention expected: - Pressure files: {run_name}.out.press.{file_num:05d}.pfb - Output files: {run_name}.out.{outflow/q_east/q_north}.{file_num:05d}.pfb