rivabar.river

rivabar.river

River(fname=None, dirname=None, start_x=None, start_y=None, end_x=None, end_y=None, file_type=None, **kwargs)

A River class that provides a clean, object-oriented interface to rivabar functionality.

This class wraps the existing functional API without breaking any current functionality. All existing functions remain available and unchanged.

Initialize a River object.

Parameters:

Name	Type	Description	Default
fname	str	Image filename	None
dirname	(str, optional)	Directory containing the image	None
start_x	float	Starting point coordinates (UTM)	None
start_y	float	Starting point coordinates (UTM)	None
end_x	float	Ending point coordinates (UTM)	None
end_y	float	Ending point coordinates (UTM)	None
file_type	str	Type of input file ('single_tif', 'multiple_tifs', etc.)	None
**kwargs	dict	Additional parameters for map_river_banks function	{}

Source code in rivabar/river.py

def __init__(self, fname=None, dirname=None, start_x=None, start_y=None, 
             end_x=None, end_y=None, file_type=None, **kwargs):
    """
    Initialize a River object.

    Parameters
    ----------
    fname : str, optional
        Image filename
    dirname : str, optional  
        Directory containing the image
    start_x, start_y : float, optional
        Starting point coordinates (UTM)
    end_x, end_y : float, optional
        Ending point coordinates (UTM)
    file_type : str, optional
        Type of input file ('single_tif', 'multiple_tifs', etc.)
    **kwargs : dict
        Additional parameters for map_river_banks function
    """

    # Store initialization parameters
    self.fname = fname
    self.dirname = dirname
    self.start_x = start_x
    self.start_y = start_y
    self.end_x = end_x
    self.end_y = end_y
    self.file_type = file_type
    self.kwargs = kwargs

    # Graph storage - will be populated after processing
    self._D_primal = None
    self._G_rook = None  
    self._G_primal = None
    self._mndwi = None
    self._dataset = None

    # Coordinate bounds
    self._left_utm_x = None
    self._right_utm_x = None
    self._lower_utm_y = None
    self._upper_utm_y = None
    self._xs = None
    self._ys = None

    # Processing flags
    self._is_processed = False
    self._processing_successful = False

    # Add new attributes
    self._mndwi_created = False
    self._delta_x = None
    self._delta_y = None

directed_graph property

Get the directed centerline graph (D_primal).

centerline_graph property

Get the primal centerline graph (G_primal).

bankline_graph property

Get the bankline (rook) graph (G_rook).

mndwi property

Get the MNDWI water mask.

dataset property

Get the raster dataset.

bounds property

Get the UTM coordinate bounds as a dict.

centerline_coords property

Get the main centerline coordinates.

main_path property

Get the main path as a list of edge tuples.

main_channel_centerline property

Get main channel centerline coordinates as a LineString.

main_channel_banks property

Get main channel bank coordinates as LineStrings.

tributary_confluences property

Get tributary confluence locations identified before dead-end removal.

Returns a list of dicts, each with: - 'confluence_utm_coords': (x, y) UTM coordinates of the confluence point - 'confluence_pixel_coords': (row, col) pixel coordinates - 'branch_utm_coords': Nx2 array of UTM coordinates along the tributary - 'branch_pixel_coords': Nx2 array of pixel coordinates along the tributary - 'branch_length_pixels': total length of the branch in pixels (measured along edge geometry, not node count)

Sorted by branch length (largest tributaries first).

has_mndwi property

Check if MNDWI has been created.

map_river_banks(**override_kwargs)

Extract river centerlines and banklines using the existing map_river_banks function.

Parameters:

Name	Type	Description	Default
**override_kwargs	dict	Parameters to override the initialization kwargs	{}

Returns:

Type	Description
bool	True if processing was successful, False otherwise

Source code in rivabar/river.py

def map_river_banks(self, **override_kwargs):
    """
    Extract river centerlines and banklines using the existing map_river_banks function.

    Parameters
    ----------
    **override_kwargs : dict
        Parameters to override the initialization kwargs

    Returns
    -------
    bool
        True if processing was successful, False otherwise
    """

    if not self._validate_inputs():
        raise ValueError("Missing required parameters. Need fname, dirname, start/end coordinates, and file_type.")

    # Merge initialization kwargs with any overrides
    processing_kwargs = {**self.kwargs, **override_kwargs}

    print(f"Processing river: {self.fname}")
    print(f"Start point: ({self.start_x}, {self.start_y})")
    print(f"End point: ({self.end_x}, {self.end_y})")

    try:
        # Call the existing map_river_banks function
        result = map_river_banks(
            fname=self.fname,
            dirname=self.dirname,
            start_x=self.start_x,
            start_y=self.start_y, 
            end_x=self.end_x,
            end_y=self.end_y,
            file_type=self.file_type,
            **processing_kwargs
        )

        # Unpack results
        (self._D_primal, self._G_rook, self._G_primal, self._mndwi, self._dataset,
         self._left_utm_x, self._right_utm_x, self._lower_utm_y, self._upper_utm_y, 
         self._xs, self._ys) = result

        # Check if processing was successful
        if self._D_primal is not None:
            self._processing_successful = True
            self._is_processed = True
            print(f"✓ Successfully processed {self.fname}")
            print(f"  - Directed graph: {len(self._D_primal.nodes)} nodes, {len(self._D_primal.edges)} edges")
            print(f"  - Rook graph: {len(self._G_rook.nodes)} polygons")
            print(f"  - Primal graph: {len(self._G_primal.nodes)} nodes")
            return True
        else:
            self._processing_successful = False
            self._is_processed = True
            print(f"✗ Failed to process {self.fname}")
            return False

    except Exception as e:
        import traceback

        print(f"✗ Error processing {self.fname}: {str(e)}")
        print(f"Error type: {type(e).__name__}")
        print(f"Processing parameters:")
        print(f"  - File type: {self.file_type}")
        print(f"  - Start point: ({self.start_x}, {self.start_y})")
        print(f"  - End point: ({self.end_x}, {self.end_y})")
        print(f"  - Directory: {self.dirname}")
        print("Full traceback:")
        traceback.print_exc()

        self._processing_successful = False
        self._is_processed = True
        return False

split_main_path_at_points(split_points)

Split the main path into segments at the given UTM coordinates.

Each split point is projected onto the nearest vertex along the main path centerline. Edges are split within their geometry when needed, so the split is accurate even when the main path consists of very few (or a single) long edges.

The method creates lightweight synthetic edge entries in the D_primal graph for the sub-edges produced by each split so that the returned sub-paths work with get_channel_widths(path=...).

Parameters:

Name	Type	Description	Default
split_points	list of tuple	UTM coordinates (x, y) at which to split. These typically come from :func:find_common_confluences. Order does not matter; they are sorted by along-channel distance internally.	required

Returns:

Name	Type	Description
segments	list of list	Exactly len(split_points) + 1 sub-paths in along-channel order, so segments[i] always lies between the i-th and (i+1)-th split point (along-channel order). Each element is a list of (s, e, d) edge tuples that can be passed to get_channel_widths(path=...) or the pairwise analysis functions. An element can be an empty list when two split points snap to the same vertex or a split point snaps to the path's start or end.
split_info	list of dict	For each split point (in along-channel order), a dict with: - 'utm_coords': the requested split point - 'snapped_utm_coords': the nearest point on the centerline - 'along_channel_distance': cumulative distance along the path - 'snapping_distance': distance from requested to snapped point - 'input_index': index of this split point in the input split_points list (the input order can differ from the along-channel order)

Source code in rivabar/river.py

def split_main_path_at_points(self, split_points):
    """
    Split the main path into segments at the given UTM coordinates.

    Each split point is projected onto the nearest vertex along the main
    path centerline. Edges are split within their geometry when needed,
    so the split is accurate even when the main path consists of very
    few (or a single) long edges.

    The method creates lightweight synthetic edge entries in the
    ``D_primal`` graph for the sub-edges produced by each split so that
    the returned sub-paths work with ``get_channel_widths(path=...)``.

    Parameters
    ----------
    split_points : list of tuple
        UTM coordinates ``(x, y)`` at which to split. These typically come
        from :func:`find_common_confluences`. Order does not matter; they
        are sorted by along-channel distance internally.

    Returns
    -------
    segments : list of list
        Exactly ``len(split_points) + 1`` sub-paths in along-channel
        order, so ``segments[i]`` always lies between the i-th and
        (i+1)-th split point (along-channel order). Each element is a
        list of ``(s, e, d)`` edge tuples that can be passed to
        ``get_channel_widths(path=...)`` or the pairwise analysis
        functions. An element can be an empty list when two split points
        snap to the same vertex or a split point snaps to the path's
        start or end.
    split_info : list of dict
        For each split point (in along-channel order), a dict with:
        - 'utm_coords': the requested split point
        - 'snapped_utm_coords': the nearest point on the centerline
        - 'along_channel_distance': cumulative distance along the path
        - 'snapping_distance': distance from requested to snapped point
        - 'input_index': index of this split point in the input
          ``split_points`` list (the input order can differ from the
          along-channel order)
    """
    self._check_processed()
    path = self.main_path
    if path is None:
        raise ValueError("No main path available.")

    # Build a continuous coordinate array and track which edge / local
    # index each vertex belongs to
    all_x, all_y = [], []
    edge_indices = []   # index into path
    local_indices = []  # index within that edge's geometry
    for idx, (s, e, d) in enumerate(path):
        geom = self._D_primal[s][e][d]['geometry']
        xs, ys = geom.xy
        n = len(xs)
        all_x.extend(xs)
        all_y.extend(ys)
        edge_indices.extend([idx] * n)
        local_indices.extend(range(n))

    all_x = np.array(all_x)
    all_y = np.array(all_y)
    edge_indices = np.array(edge_indices)
    local_indices = np.array(local_indices)

    # Cumulative along-channel distance
    dx = np.diff(all_x)
    dy = np.diff(all_y)
    ds = np.sqrt(dx**2 + dy**2)
    s_cum = np.concatenate([[0], np.cumsum(ds)])

    # For each split point, find the nearest centerline vertex
    centerline_pts = np.column_stack([all_x, all_y])
    split_data = []
    for input_idx, sp in enumerate(split_points):
        sp = np.array(sp)
        dists = np.sqrt(np.sum((centerline_pts - sp)**2, axis=1))
        nearest_idx = int(np.argmin(dists))
        split_data.append({
            'utm_coords': (float(sp[0]), float(sp[1])),
            'snapped_utm_coords': (float(all_x[nearest_idx]),
                                   float(all_y[nearest_idx])),
            'along_channel_distance': float(s_cum[nearest_idx]),
            'snapping_distance': float(dists[nearest_idx]),
            'input_index': input_idx,
            'global_idx': nearest_idx,
            'edge_index': int(edge_indices[nearest_idx]),
            'local_index': int(local_indices[nearest_idx]),
        })

    # Sort by along-channel distance (input order may differ)
    split_data.sort(key=lambda d: d['along_channel_distance'])

    # ---- Build segments by slicing edges at the split vertices ----------
    # We process the path edge by edge. For each edge we track which split
    # points fall inside it (by edge_index) and slice accordingly.
    splits_by_edge = {}  # edge_index -> list of local_index values
    for sd in split_data:
        ei = sd['edge_index']
        li = sd['local_index']
        splits_by_edge.setdefault(ei, []).append(li)

    # Ensure unique ids for synthetic edges we add to the graph
    max_key = 0
    for s_node, e_node, k in self._D_primal.edges(keys=True):
        if isinstance(k, int) and k > max_key:
            max_key = k
    next_key = max_key + 1000  # offset to avoid collisions

    segments = [[]]  # start with one empty segment; new ones added at splits

    for path_idx, (s_node, e_node, d_key) in enumerate(path):
        edge_data = self._D_primal[s_node][e_node][d_key]
        geom_coords = list(edge_data['geometry'].coords)
        hw_keys = list(edge_data['half_widths'].keys())
        hw0 = list(edge_data['half_widths'][hw_keys[0]])
        hw1 = list(edge_data['half_widths'][hw_keys[1]])
        n_pts = len(geom_coords)

        if path_idx not in splits_by_edge:
            # No splits in this edge — add it whole to current segment
            segments[-1].append((s_node, e_node, d_key))
            continue

        # Split this edge at the local indices. Duplicates are kept so
        # that every split point produces exactly one segment boundary
        # (coincident split points yield empty segments) and segments
        # stay positionally aligned with the split points.
        cut_indices = sorted(splits_by_edge[path_idx])

        # Slice the edge into pieces. Boundaries are inclusive vertex
        # indices; the boundary vertex is shared between consecutive
        # pieces so there are no gaps. A cut at the edge's first/last
        # vertex yields a single-point piece, which is skipped, placing
        # the segment boundary at the edge's end without losing it.
        boundaries = [0] + cut_indices + [n_pts - 1]
        for i in range(len(boundaries) - 1):
            start = boundaries[i]
            end = boundaries[i + 1]
            sub_coords = geom_coords[start:end + 1]
            sub_hw0 = hw0[start:end + 1]
            sub_hw1 = hw1[start:end + 1]

            if len(sub_coords) >= 2:
                # Create a synthetic edge in D_primal for this sub-segment
                syn_key = next_key
                next_key += 1
                self._D_primal.add_edge(s_node, e_node, key=syn_key,
                    geometry=LineString(sub_coords),
                    half_widths={hw_keys[0]: sub_hw0, hw_keys[1]: sub_hw1},
                    mm_len=edge_data.get('mm_len', 0),
                    width=edge_data.get('width', 0),
                    _synthetic=True,
                )
                segments[-1].append((s_node, e_node, syn_key))

            # Start a new segment after each cut (except the last piece)
            if i < len(boundaries) - 2:
                segments.append([])

    # Empty segments are kept: segments[i] must always lie between split
    # points i-1 and i (along-channel order) for callers that pair
    # segments with confluences by position

    # Build split_info
    split_info = []
    for sd in split_data:
        split_info.append({
            'utm_coords': sd['utm_coords'],
            'snapped_utm_coords': sd['snapped_utm_coords'],
            'along_channel_distance': sd['along_channel_distance'],
            'snapping_distance': sd['snapping_distance'],
            'input_index': sd['input_index'],
        })

    return segments, split_info

get_channel_widths(path=None, pixel_size=None)

Get channel widths along the main path.

Parameters:

Name	Type	Description	Default
path	list	Custom path as list of edge tuples. If None, uses main path.	None
pixel_size	float	Pixel size in meters, used to convert widths to meters. If None, taken from the dataset transform; required for River objects whose raster data has been cleared (e.g., loaded from older pickles that did not save the transform).	None

Returns:

Type	Description
ndarray	Channel widths along the path

Source code in rivabar/river.py

def get_channel_widths(self, path=None, pixel_size=None):
    """
    Get channel widths along the main path.

    Parameters
    ----------
    path : list, optional
        Custom path as list of edge tuples. If None, uses main path.
    pixel_size : float, optional
        Pixel size in meters, used to convert widths to meters. If None,
        taken from the dataset transform; required for River objects
        whose raster data has been cleared (e.g., loaded from older
        pickles that did not save the transform).

    Returns
    -------
    np.ndarray
        Channel widths along the path
    """
    self._check_processed()

    if path is None:
        path = self.main_path
    if path is None:
        raise ValueError("No main path available and no custom path provided.")

    if pixel_size is None:
        if self._dataset is not None and getattr(self._dataset, 'transform', None) is not None:
            pixel_size = self._dataset.transform[0]
        else:
            raise ValueError("No dataset transform available (raster data was cleared, "
                             "e.g. on a River loaded from an older pickle); "
                             "pass pixel_size explicitly.")

    xl, yl, w1l, w2l, w, s = get_channel_widths_along_path(self._D_primal, path)

    return s, np.array(w)*pixel_size # convert to meters

analyze_wavelength_and_width(path=None, ax=None, **kwargs)

Analyze channel wavelength and width statistics.

Parameters:

Name	Type	Description	Default
path	list	Custom path as list of edge tuples. If None, uses main path.	None
ax	Axes	Axes to plot on. If None, creates new figure.	None

Returns:

Type	Description
dict	Dictionary with analysis results

Source code in rivabar/river.py

def analyze_wavelength_and_width(self, path=None, ax=None, **kwargs):
    """
    Analyze channel wavelength and width statistics.

    Parameters
    ----------
    path : list, optional
        Custom path as list of edge tuples. If None, uses main path.
    ax : matplotlib.Axes, optional
        Axes to plot on. If None, creates new figure.

    Returns
    -------
    dict
        Dictionary with analysis results
    """
    self._check_processed()

    if path is None:
        path = self.main_path
    if path is None:
        raise ValueError("No main path available and no custom path provided.")

    return analyze_width_and_wavelength(self._D_primal, path, ax, **kwargs)

plot_overview(figsize=(12, 8), **kwargs)

Plot an overview of the river with centerlines and banks.

Parameters:

Name	Type	Description	Default
figsize	tuple	Figure size (width, height)	(12, 8)
**kwargs	dict	Additional arguments passed to plot_im_and_lines	{}

Returns:

Type	Description
(fig, ax)	Matplotlib figure and axes objects

Source code in rivabar/river.py

def plot_overview(self, figsize=(12, 8), **kwargs):
    """
    Plot an overview of the river with centerlines and banks.

    Parameters
    ----------
    figsize : tuple, optional
        Figure size (width, height)
    **kwargs : dict
        Additional arguments passed to plot_im_and_lines

    Returns
    -------
    fig, ax
        Matplotlib figure and axes objects
    """
    self._check_processed()

    fig, ax = plot_im_and_lines(
        self._mndwi, 
        self._left_utm_x, self._right_utm_x, 
        self._lower_utm_y, self._upper_utm_y,
        self._G_rook, self._G_primal, self._D_primal,
        start_x=self.start_x, start_y=self.start_y,
        end_x=self.end_x, end_y=self.end_y,
        **kwargs
    )

    ax.set_title(f'River: {self.fname}')
    fig.set_size_inches(figsize)

    return fig, ax

plot_tributaries(ax=None, color='r', linewidth=2, markersize=8, min_length=None, label=True)

Plot tributary branches and their confluence points.

Parameters:

Name	Type	Description	Default
ax	Axes	Axes to plot on. If None, creates an overview plot first.	None
color	str	Color for tributary lines and markers (default 'r').	'r'
linewidth	float	Line width for tributary branches (default 2).	2
markersize	float	Marker size for confluence points (default 8).	8
min_length	float	If set, only plot tributaries with branch_length_pixels >= this value.	None
label	bool	If True, annotate each confluence with the tributary's pixel length (default True).	True

Returns:

Name	Type	Description
ax	Axes

Source code in rivabar/river.py

def plot_tributaries(self, ax=None, color='r', linewidth=2, markersize=8,
                     min_length=None, label=True):
    """
    Plot tributary branches and their confluence points.

    Parameters
    ----------
    ax : matplotlib.Axes, optional
        Axes to plot on. If None, creates an overview plot first.
    color : str, optional
        Color for tributary lines and markers (default 'r').
    linewidth : float, optional
        Line width for tributary branches (default 2).
    markersize : float, optional
        Marker size for confluence points (default 8).
    min_length : float, optional
        If set, only plot tributaries with branch_length_pixels >= this value.
    label : bool, optional
        If True, annotate each confluence with the tributary's pixel length
        (default True).

    Returns
    -------
    ax : matplotlib.Axes
    """
    self._check_processed()
    tribs = self.tributary_confluences
    if not tribs:
        print("No tributaries found.")
        return ax

    if ax is None:
        fig, ax = self.plot_overview()
        plot_graph_w_colors(self._D_primal, ax)

    for trib in tribs:
        length = trib['branch_length_pixels']
        if min_length is not None and length < min_length:
            continue
        if 'branch_utm_coords' in trib:
            coords = trib['branch_utm_coords']
            ax.plot(coords[:, 0], coords[:, 1], '-', color=color,
                    linewidth=linewidth)
        if 'confluence_utm_coords' in trib:
            cx, cy = trib['confluence_utm_coords']
            ax.plot(cx, cy, 'o', color=color, markersize=markersize)
            if label:
                ax.annotate(f'{length:.0f} px', (cx, cy),
                            textcoords='offset points', xytext=(5, 5),
                            fontsize=8, color=color)
    return ax

plot_directed_graph(ax=None, **kwargs)

Plot the directed centerline graph with colors.

Parameters:

Name	Type	Description	Default
ax	Axes	Axes to plot on. If None, uses current axes.	None
**kwargs	dict	Additional arguments passed to plot_graph_w_colors	{}

Returns:

Type	Description
ax	Matplotlib axes object

Source code in rivabar/river.py

def plot_directed_graph(self, ax=None, **kwargs):
    """
    Plot the directed centerline graph with colors.

    Parameters
    ----------
    ax : matplotlib.Axes, optional
        Axes to plot on. If None, uses current axes.
    **kwargs : dict
        Additional arguments passed to plot_graph_w_colors

    Returns
    -------
    ax
        Matplotlib axes object
    """
    self._check_processed()

    if ax is None:
        ax = plt.gca()

    return plot_graph_w_colors(self._D_primal, ax, **kwargs)

create_mndwi(mndwi_threshold=0.01, delete_pixels_polys=False, small_hole_threshold=64, remove_smaller_components=True, solidity_filter=False, **kwargs)

Create MNDWI water mask from the river's input data.

This allows you to create and visualize the water mask before running the full centerline extraction process.

Parameters:

Name	Type	Description	Default
mndwi_threshold	float	Threshold for water detection. Default 0.01.	0.01
delete_pixels_polys	bool or list	Polygons to mask out (e.g., bridges). Default False.	False
small_hole_threshold	int	Minimum hole size to keep. Default 64.	64
remove_smaller_components	bool	Remove small water bodies. Default True.	True
solidity_filter	bool	Filter by solidity. Default False.	False
**kwargs	dict	Additional parameters for create_mndwi	{}

Returns:

Type	Description
ndarray	Binary water mask

Source code in rivabar/river.py

def create_mndwi(self, mndwi_threshold=0.01, delete_pixels_polys=False, 
                 small_hole_threshold=64, remove_smaller_components=True, 
                 solidity_filter=False, **kwargs):
    """
    Create MNDWI water mask from the river's input data.

    This allows you to create and visualize the water mask before running
    the full centerline extraction process.

    Parameters
    ----------
    mndwi_threshold : float, optional
        Threshold for water detection. Default 0.01.
    delete_pixels_polys : bool or list, optional
        Polygons to mask out (e.g., bridges). Default False.
    small_hole_threshold : int, optional
        Minimum hole size to keep. Default 64.
    remove_smaller_components : bool, optional
        Remove small water bodies. Default True.
    solidity_filter : bool, optional
        Filter by solidity. Default False.
    **kwargs : dict
        Additional parameters for create_mndwi

    Returns
    -------
    numpy.ndarray
        Binary water mask
    """

    if not self._validate_inputs():
        raise ValueError("Missing required parameters. Need fname, dirname, and file_type.")

    print(f"Creating MNDWI for: {self.fname}")

    # Call the create_mndwi function
    result = create_mndwi(
        dirname=self.dirname,
        fname=self.fname,
        file_type=self.file_type,
        mndwi_threshold=mndwi_threshold,
        delete_pixels_polys=delete_pixels_polys,
        small_hole_threshold=small_hole_threshold,
        remove_smaller_components=remove_smaller_components,
        solidity_filter=solidity_filter,
        **kwargs
    )

    # Unpack and store results
    (self._mndwi, self._left_utm_x, self._upper_utm_y, self._right_utm_x, 
     self._lower_utm_y, self._delta_x, self._delta_y, self._dataset) = result

    # Update processing flags
    self._mndwi_created = True

    print(f"✓ MNDWI created: {self._mndwi.shape}")
    print(f"  - Water pixels: {np.sum(self._mndwi):,}")
    print(f"  - Total pixels: {self._mndwi.size:,}")
    print(f"  - Water percentage: {100*np.sum(self._mndwi)/self._mndwi.size:.1f}%")

    return self._mndwi

plot_mndwi(figsize=(12, 8), add_start_end_points=True, **kwargs)

Plot the MNDWI water mask.

Parameters:

Name	Type	Description	Default
figsize	tuple	Figure size. Default (12, 8).	(12, 8)
add_start_end_points	bool	Whether to plot start/end points. Default True.	True
**kwargs	dict	Additional arguments for plt.imshow	{}

Returns:

Type	Description
(fig, ax)	Matplotlib figure and axes

Source code in rivabar/river.py

def plot_mndwi(self, figsize=(12, 8), add_start_end_points=True, **kwargs):
    """
    Plot the MNDWI water mask.

    Parameters
    ----------
    figsize : tuple, optional
        Figure size. Default (12, 8).
    add_start_end_points : bool, optional
        Whether to plot start/end points. Default True.
    **kwargs : dict
        Additional arguments for plt.imshow

    Returns
    -------
    fig, ax
        Matplotlib figure and axes
    """

    if self._mndwi is None:
        raise RuntimeError("MNDWI not created yet. Call create_mndwi() first.")

    fig, ax = plt.subplots(figsize=figsize)

    # Plot MNDWI
    ax.imshow(
        self._mndwi, 
        extent=[self._left_utm_x, self._right_utm_x, 
               self._lower_utm_y, self._upper_utm_y],
        cmap='gray_r',
        **kwargs
    )

    # Add start/end points if provided
    if add_start_end_points and all([self.start_x, self.start_y, self.end_x, self.end_y]):
        ax.plot(self.start_x, self.start_y, 'go', markersize=8, label='Start')
        ax.plot(self.end_x, self.end_y, 'ro', markersize=8, label='End')
        ax.legend()

    ax.set_xlabel('UTM X (m)')
    ax.set_ylabel('UTM Y (m)')
    ax.set_title(f'MNDWI Water Mask: {self.fname}')

    return fig, ax

get_start_end_points_interactive()

Interactively select start and end points on the MNDWI plot.

Returns:

Type	Description
tuple	(start_x, start_y, end_x, end_y)

Source code in rivabar/river.py

def get_start_end_points_interactive(self):
    """
    Interactively select start and end points on the MNDWI plot.

    Returns
    -------
    tuple
        (start_x, start_y, end_x, end_y)
    """

    if self._mndwi is None:
        raise RuntimeError("MNDWI not created yet. Call create_mndwi() first.")

    fig, ax = self.plot_mndwi(add_start_end_points=False)

    print("Click two points: first for START, second for END")
    points = plt.ginput(n=2, timeout=-1)  # Set timeout to -1 to prevent timing out

    if len(points) == 2:
        self.start_x, self.start_y = points[0]
        self.end_x, self.end_y = points[1]

        # Add points to plot
        ax.plot(self.start_x, self.start_y, 'go', markersize=8, label='Start')
        ax.plot(self.end_x, self.end_y, 'ro', markersize=8, label='End')
        ax.legend()
        plt.draw()

        print(f"Start point: ({self.start_x:.1f}, {self.start_y:.1f})")
        print(f"End point: ({self.end_x:.1f}, {self.end_y:.1f})")

        return self.start_x, self.start_y, self.end_x, self.end_y
    else:
        print("Need exactly 2 points!")
        return None

summary()

Print a summary of the river analysis.

Source code in rivabar/river.py

def summary(self):
    """Print a summary of the river analysis."""
    print(f"\n=== River Summary: {self.fname} ===")
    print(f"Processed: {self._is_processed}")
    print(f"Successful: {self._processing_successful}")

    if self._processing_successful:
        print(f"Start point: ({self.start_x:.1f}, {self.start_y:.1f})")
        print(f"End point: ({self.end_x:.1f}, {self.end_y:.1f})")
        print(f"Directed graph: {len(self._D_primal.nodes)} nodes, {len(self._D_primal.edges)} edges")
        print(f"Rook graph: {len(self._G_rook.nodes)} polygons") 
        print(f"Primal graph: {len(self._G_primal.nodes)} nodes")
        if self._mndwi is not None:
            print(f"MNDWI shape: {self._mndwi.shape}")

        # Main path info
        if self.main_path:
            print(f"Main path: {len(self.main_path)} edges")
            try:
                s, widths = self.get_channel_widths()
                print(f"Channel width: {np.mean(widths):.1f} ± {np.std(widths):.1f} m")
                print(f"Channel length: {s[-1]:.1f} m")
            except ValueError as e:
                print(f"Channel width: unavailable ({e})")

    print("=" * (20 + len(self.fname)))

collect_stats(pixel_size=None)

Collect summary statistics for this river.

Thin wrapper around :func:collect_river_stats; see its docstring for the returned fields.

Source code in rivabar/river.py

def collect_stats(self, pixel_size=None):
    """
    Collect summary statistics for this river.

    Thin wrapper around :func:`collect_river_stats`; see its docstring
    for the returned fields.
    """
    return collect_river_stats(self, pixel_size=pixel_size)

clear_raster_data()

Clear heavy raster data to save memory while keeping graph results.

This removes the MNDWI array but preserves all graph data, coordinates, and analysis results.

Source code in rivabar/river.py

def clear_raster_data(self):
    """
    Clear heavy raster data to save memory while keeping graph results.

    This removes the MNDWI array but preserves all graph data,
    coordinates, and analysis results.
    """
    # Clear heavy raster data
    self._mndwi = None

    # Keep everything else: graphs, coordinates, bounds, etc.
    print(f"Cleared raster data for {self.fname}")

to_geopandas()

Export river data to GeoDataFrames.

Returns:

Type	Description
dict	Dictionary containing GeoDataFrames for different components

Source code in rivabar/river.py

def to_geopandas(self):
    """
    Export river data to GeoDataFrames.

    Returns
    -------
    dict
        Dictionary containing GeoDataFrames for different components
    """
    self._check_processed()

    crs = getattr(self._dataset, 'crs', None) if self._dataset is not None else None
    if crs is None:
        print("Warning: no CRS available (raster data was cleared); "
              "GeoDataFrames will have no CRS")

    result = {}

    # Main centerline
    if self.main_channel_centerline:
        result['centerline'] = gpd.GeoDataFrame(
            [{'geometry': self.main_channel_centerline, 'name': self.fname}],
            crs=crs
        )

    # Banks
    banks = self.main_channel_banks
    if banks:
        bank_data = [
            {'geometry': banks['left_bank'], 'side': 'left', 'name': self.fname},
            {'geometry': banks['right_bank'], 'side': 'right', 'name': self.fname}
        ]
        result['banks'] = gpd.GeoDataFrame(bank_data, crs=crs)

    # Channel polygons (from rook graph)
    polygons = []
    for node_id, data in self._G_rook.nodes(data=True):
        if 'bank_polygon' in data and data['bank_polygon'] is not None:
            polygons.append({
                'geometry': data['bank_polygon'],
                'node_id': node_id,
                'name': self.fname
            })

    if polygons:
        result['polygons'] = gpd.GeoDataFrame(polygons, crs=crs)

    return result

save_results(filepath)

Save river analysis results to file.

Source code in rivabar/river.py

def save_results(self, filepath):
    """Save river analysis results to file."""
    import pickle

    # Create a results dictionary without heavy raster data
    results = {
        'fname': self.fname,
        'dirname': self.dirname,
        'start_x': self.start_x, 'start_y': self.start_y,
        'end_x': self.end_x, 'end_y': self.end_y,
        'file_type': self.file_type,
        'kwargs': self.kwargs,

        # Graph data (these are the important results)
        'D_primal': self._D_primal,
        'G_rook': self._G_rook,
        'G_primal': self._G_primal,

        # Coordinate data
        'left_utm_x': self._left_utm_x,
        'right_utm_x': self._right_utm_x,
        'lower_utm_y': self._lower_utm_y,
        'upper_utm_y': self._upper_utm_y,
        'xs': self._xs,
        'ys': self._ys,

        # Processing flags
        'processing_successful': self._processing_successful,
        'is_processed': self._is_processed,

        # Dataset metadata (so loaded rivers can convert widths to
        # meters and export with a CRS without the raster data)
        'dataset_crs': (str(self._dataset.crs)
                        if self._dataset is not None and getattr(self._dataset, 'crs', None) is not None
                        else None),
        'dataset_transform': (self._dataset.transform
                              if self._dataset is not None and getattr(self._dataset, 'transform', None) is not None
                              else None),
        'dataset_shape': (self._dataset.shape
                          if self._dataset is not None and getattr(self._dataset, 'shape', None) is not None
                          else None),
    }

    with open(filepath, 'wb') as f:
        pickle.dump(results, f)

    print(f"Saved results to {filepath}")

load_results(filepath) classmethod

Load river results from file and create a River instance.

Source code in rivabar/river.py

@classmethod
def load_results(cls, filepath):
    """Load river results from file and create a River instance."""
    import pickle

    with open(filepath, 'rb') as f:
        results = pickle.load(f)

    # Create river instance
    river = cls(
        fname=results['fname'],
        dirname=results['dirname'],
        start_x=results['start_x'],
        start_y=results['start_y'],
        end_x=results['end_x'],
        end_y=results['end_y'],
        file_type=results['file_type'],
        **results['kwargs']
    )

    # Restore processed data
    river._D_primal = results['D_primal']
    river._G_rook = results['G_rook']
    river._G_primal = results['G_primal']
    river._left_utm_x = results['left_utm_x']
    river._right_utm_x = results['right_utm_x']
    river._lower_utm_y = results['lower_utm_y']
    river._upper_utm_y = results['upper_utm_y']
    river._xs = results['xs']
    river._ys = results['ys']
    river._processing_successful = results['processing_successful']
    river._is_processed = results['is_processed']

    # MNDWI is not loaded (saving memory); rebuild a minimal dataset
    # from the saved metadata so width conversion and CRS export work
    river._mndwi = None
    if results.get('dataset_crs') or results.get('dataset_transform'):
        from .data_io import MinimalDataset
        river._dataset = MinimalDataset(
            results.get('dataset_crs'),
            results.get('dataset_transform'),
            results.get('dataset_shape'),
        )
    else:
        river._dataset = None

    print(f"Loaded results from {filepath}")
    return river

get_memory_usage()

Get approximate memory usage of stored data.

Source code in rivabar/river.py

def get_memory_usage(self):
    """Get approximate memory usage of stored data."""
    memory_mb = 0
    if self._mndwi is not None:
        memory_mb += self._mndwi.nbytes / 1024**2
    if self._dataset is not None:
        # Rough estimate for dataset
        memory_mb += 50  # Approximate

    return f"{memory_mb:.1f} MB"

get_memory_breakdown()

Get detailed memory usage breakdown.

Source code in rivabar/river.py

def get_memory_breakdown(self):
    """Get detailed memory usage breakdown."""
    breakdown = {}
    total_mb = 0

    if self._mndwi is not None:
        mndwi_mb = self._mndwi.nbytes / 1024**2
        breakdown['MNDWI'] = f"{mndwi_mb:.1f} MB"
        total_mb += mndwi_mb
    else:
        breakdown['MNDWI'] = "Not loaded"

    if self._dataset is not None:
        dataset_mb = 50  # Rough estimate
        breakdown['Dataset'] = f"{dataset_mb:.1f} MB"
        total_mb += dataset_mb
    else:
        breakdown['Dataset'] = "Not loaded"

    # Graph data is typically small
    if self._D_primal is not None:
        breakdown['Graphs'] = "~1-5 MB"
        total_mb += 3  # Rough estimate
    else:
        breakdown['Graphs'] = "Not processed"

    breakdown['Total'] = f"{total_mb:.1f} MB"
    return breakdown

batch_process_landsat_scenes(path_number, row_number, start_x, start_y, end_x, end_y, years=None, max_cloud_cover=10, n_scenes_per_year=3, clear_rasters=True, save_individual=True, save_dir='river_results', skip_scenes=None, download_mndwi=False, mndwi_dir='mndwi_rasters', download_false_color=False, false_color_dir='false_color_images', **processing_kwargs) classmethod

Batch process multiple Landsat scenes from Google Earth Engine.

Parameters:

Name	Type	Description	Default
path_number	int	Landsat WRS path number	required
row_number	int	Landsat WRS row number	required
start_x	float	Starting point coordinates (UTM)	required
start_y	float	Starting point coordinates (UTM)	required
end_x	float	Ending point coordinates (UTM)	required
end_y	float	Ending point coordinates (UTM)	required
years	list or range	Years to process. If None, uses 1984-2023	None
max_cloud_cover	float	Maximum cloud cover percentage (default 10)	10
n_scenes_per_year	int	Maximum scenes per year (default 3)	3
clear_rasters	bool	Whether to clear raster data after processing (default True)	True
save_individual	bool	Whether to save each river individually to prevent data loss (default True)	True
save_dir	str	Directory to save individual river files (default 'river_results')	'river_results'
skip_scenes	list	List of scene IDs to skip (useful for avoiding problematic scenes that crash the kernel)	None
download_mndwi	bool	Whether to download MNDWI rasters to local folder (default False)	False
mndwi_dir	str	Directory to save MNDWI rasters (default 'mndwi_rasters')	'mndwi_rasters'
download_false_color	bool	Whether to download false color images (SWIR2-NIR-Green) to local folder (default False)	False
false_color_dir	str	Directory to save false color images (default 'false_color_images')	'false_color_images'
**processing_kwargs	dict	Additional arguments for map_river_banks	{}

Returns:

Name	Type	Description
successful_rivers	list	List of successfully processed River instances (or file paths if save_individual=True)
failed_scenes	list	List of scene IDs that failed processing

Source code in rivabar/river.py

@classmethod
def batch_process_landsat_scenes(cls, path_number, row_number, start_x, start_y, end_x, end_y,
                               years=None, max_cloud_cover=10, n_scenes_per_year=3, 
                               clear_rasters=True, save_individual=True, save_dir='river_results',
                               skip_scenes=None, download_mndwi=False, mndwi_dir='mndwi_rasters',
                               download_false_color=False, false_color_dir='false_color_images',
                               **processing_kwargs):
    """
    Batch process multiple Landsat scenes from Google Earth Engine.

    Parameters
    ----------
    path_number : int
        Landsat WRS path number
    row_number : int
        Landsat WRS row number  
    start_x, start_y : float
        Starting point coordinates (UTM)
    end_x, end_y : float
        Ending point coordinates (UTM)
    years : list or range, optional
        Years to process. If None, uses 1984-2023
    max_cloud_cover : float, optional
        Maximum cloud cover percentage (default 10)
    n_scenes_per_year : int, optional
        Maximum scenes per year (default 3)
    clear_rasters : bool, optional
        Whether to clear raster data after processing (default True)
    save_individual : bool, optional
        Whether to save each river individually to prevent data loss (default True)
    save_dir : str, optional
        Directory to save individual river files (default 'river_results')
    skip_scenes : list, optional
        List of scene IDs to skip (useful for avoiding problematic scenes that crash the kernel)
    download_mndwi : bool, optional
        Whether to download MNDWI rasters to local folder (default False)
    mndwi_dir : str, optional
        Directory to save MNDWI rasters (default 'mndwi_rasters')
    download_false_color : bool, optional
        Whether to download false color images (SWIR2-NIR-Green) to local folder (default False)
    false_color_dir : str, optional
        Directory to save false color images (default 'false_color_images')
    **processing_kwargs : dict
        Additional arguments for map_river_banks

    Returns
    -------
    successful_rivers : list
        List of successfully processed River instances (or file paths if save_individual=True)
    failed_scenes : list
        List of scene IDs that failed processing
    """
    import ee
    import tempfile
    import os
    import gc
    import geemap
    import pickle
    from datetime import datetime

    if years is None:
        years = range(1984, 2024)

    if skip_scenes is None:
        skip_scenes = []

    # Create save directory if saving individually
    if save_individual:
        os.makedirs(save_dir, exist_ok=True)
        print(f"💾 Individual rivers will be saved to: {save_dir}")

    # Create MNDWI directory if downloading rasters
    if download_mndwi:
        os.makedirs(mndwi_dir, exist_ok=True)
        print(f"🗺️  MNDWI rasters will be saved to: {mndwi_dir}")

    # Create false color directory if downloading images
    if download_false_color:
        os.makedirs(false_color_dir, exist_ok=True)
        print(f"🌈 False color images will be saved to: {false_color_dir}")

    if skip_scenes:
        print(f"⏭️  Will skip {len(skip_scenes)} problematic scenes: {skip_scenes}")

    successful_rivers = []
    failed_scenes = []
    skipped_scenes = []
    scene_count = 0

    print(f"🛰️  Batch processing Landsat scenes for Path {path_number}, Row {row_number}")
    print(f"Years: {min(years)}-{max(years)}, Max cloud cover: {max_cloud_cover}%")
    print(f"Processing parameters: {processing_kwargs}")

    # Process each year
    for year in years:
        start_date = f'{year}-01-01'
        end_date = f'{year}-12-31'

        print(f"\n--- Processing year {year} ---")

        # Determine which Landsat mission to use
        if year < 2013:
            # Landsat 5
            collection_id = "LANDSAT/LT05/C02/T1_L2"
            green_band = 'SR_B2'
            nir_band = 'SR_B4'
            swir_band = 'SR_B5'
            swir2_band = 'SR_B7'
        elif year < 2021:
            # Landsat 8
            collection_id = "LANDSAT/LC08/C02/T1_L2" 
            green_band = 'SR_B3'
            nir_band = 'SR_B5'
            swir_band = 'SR_B6'
            swir2_band = 'SR_B7'
        else:
            # Landsat 9
            collection_id = "LANDSAT/LC09/C02/T1_L2"
            green_band = 'SR_B3'
            nir_band = 'SR_B5'
            swir_band = 'SR_B6'
            swir2_band = 'SR_B7'

        # Get image collection
        try:
            collection = (ee.ImageCollection(collection_id)
                        .filterDate(start_date, end_date)
                        .filter(ee.Filter.eq('WRS_PATH', path_number))
                        .filter(ee.Filter.eq('WRS_ROW', row_number))
                        .filter(ee.Filter.lt("CLOUD_COVER", max_cloud_cover))
                        .sort('CLOUD_COVER', True)
                        .limit(n_scenes_per_year))

            imgs = collection.toList(collection.size())
            n_images = collection.size().getInfo()

            print(f"Found {n_images} scenes for {year}")

            # Process each image
            for i in range(n_images):
                scene_count += 1

                try:
                    im = ee.Image(imgs.get(i))
                    scene_id = im.getString('system:index').getInfo()
                    acquisition_date = im.getString('DATE_ACQUIRED').getInfo()
                    cloud_cover = im.getNumber('CLOUD_COVER').getInfo()

                    print(f"  Processing scene {scene_count}: {scene_id}")
                    print(f"    Date: {acquisition_date}, Cloud cover: {cloud_cover:.1f}%")

                    # Check if this scene should be skipped
                    if scene_id in skip_scenes:
                        print(f"    ⏭️  Skipping {scene_id} (in skip list)")
                        skipped_scenes.append(scene_id)
                        continue

                    # Check if this scene was already processed (if saving individually)
                    if save_individual:
                        river_file = os.path.join(save_dir, f"river_{scene_id}.pkl")
                        if os.path.exists(river_file):
                            print(f"    ⏭️  Skipping {scene_id} (already processed)")
                            successful_rivers.append(river_file)
                            continue

                    # Mask clouds and cloud shadows using QA_PIXEL band
                    # Bit 3 = cloud, bit 4 = cloud shadow
                    qa = im.select('QA_PIXEL')
                    cloud_mask = (qa.bitwiseAnd(1 << 3).eq(0)    # not cloud
                                   .And(qa.bitwiseAnd(1 << 4).eq(0)))  # not cloud shadow
                    im_masked = im.updateMask(cloud_mask)

                    # Calculate MNDWI
                    mndwi = im_masked.normalizedDifference([green_band, swir_band])

                    # Download to temporary file
                    with tempfile.NamedTemporaryFile(suffix='.tif', delete=False) as tmp_file:
                        tmp_path = tmp_file.name

                    try:
                        # Download the MNDWI image
                        print(f"    Downloading MNDWI...")
                        geemap.download_ee_image(mndwi, tmp_path)

                        # Save MNDWI raster to permanent location if requested
                        if download_mndwi:
                            mndwi_filename = f"mndwi_{scene_id}.tif"
                            mndwi_path = os.path.join(mndwi_dir, mndwi_filename)

                            # Check if MNDWI file already exists
                            if not os.path.exists(mndwi_path):
                                import shutil
                                shutil.copy2(tmp_path, mndwi_path)
                                print(f"    🗺️  Saved MNDWI: {mndwi_path}")
                            else:
                                print(f"    🗺️  MNDWI already exists: {mndwi_path}")

                        # Create and save false color image if requested
                        if download_false_color:
                            false_color_filename = f"false_color_{scene_id}.tif"
                            false_color_path = os.path.join(false_color_dir, false_color_filename)

                            # Check if false color image already exists
                            if not os.path.exists(false_color_path):
                                print(f"    Creating false color image...")

                                # Create false color composite: SWIR2, NIR, Green
                                false_color = im.select([swir2_band, nir_band, green_band])

                                # Create temporary file for false color image
                                with tempfile.NamedTemporaryFile(suffix='.tif', delete=False) as tmp_fc_file:
                                    tmp_fc_path = tmp_fc_file.name

                                try:
                                    print(f"    Downloading false color image...")
                                    geemap.download_ee_image(false_color, tmp_fc_path)

                                    import shutil
                                    shutil.copy2(tmp_fc_path, false_color_path)
                                    print(f"    🌈 Saved false color: {false_color_path}")

                                    # Clean up temporary false color file
                                    os.unlink(tmp_fc_path)

                                except Exception as e:
                                    print(f"    ❌ Failed to download false color image: {e}")
                                    if os.path.exists(tmp_fc_path):
                                        os.unlink(tmp_fc_path)
                            else:
                                print(f"    🌈 False color already exists: {false_color_path}")

                        # Create River instance and process
                        print(f"    Processing with rivabar...")

                        river = cls(
                            fname=os.path.basename(tmp_path),
                            dirname=os.path.dirname(tmp_path) + '/',
                            start_x=start_x,
                            start_y=start_y,
                            end_x=end_x,
                            end_y=end_y,
                            file_type='water_index',
                            **processing_kwargs
                        )

                        # Add metadata
                        river.scene_id = scene_id
                        river.acquisition_date = acquisition_date
                        river.cloud_cover = cloud_cover
                        river.landsat_mission = collection_id.split('/')[1]
                        river.year = year

                        # Process the river
                        memory_before = river.get_memory_usage()
                        success = river.map_river_banks()

                        if success:
                            print(f"    ✅ Success! Memory: {memory_before} → {river.get_memory_usage()}")

                            # Clear raster data to save memory
                            if clear_rasters:
                                river.clear_raster_data()
                                print(f"    🧹 Cleared rasters, memory: {river.get_memory_usage()}")

                            # Save individual river if requested
                            if save_individual:
                                river_file = os.path.join(save_dir, f"river_{scene_id}.pkl")

                                # Create a safe copy of the river for pickling
                                # Store the dataset temporarily and remove it
                                temp_dataset = river._dataset
                                river._dataset = None

                                # Create a results dictionary for saving
                                river_data = {
                                    'river': river,
                                    'scene_id': scene_id,
                                    'acquisition_date': acquisition_date,
                                    'cloud_cover': cloud_cover,
                                    'landsat_mission': collection_id.split('/')[1],
                                    'year': year,
                                    'processing_timestamp': datetime.now(),
                                    'path_number': path_number,
                                    'row_number': row_number,
                                    # Store dataset info separately
                                    'dataset_crs': str(temp_dataset.crs) if temp_dataset else None,
                                    'dataset_transform': temp_dataset.transform if temp_dataset else None,
                                    'dataset_shape': temp_dataset.shape if temp_dataset else None
                                }

                                try:
                                    with open(river_file, 'wb') as f:
                                        pickle.dump(river_data, f)
                                    print(f"    💾 Saved: {river_file}")
                                    successful_rivers.append(river_file)  # Store file path
                                except Exception as pickle_error:
                                    print(f"    ⚠️ Error saving {river_file}: {pickle_error}")
                                    # If saving fails, still add to successful_rivers list
                                    successful_rivers.append(river)
                                finally:
                                    # Restore the dataset reference
                                    river._dataset = temp_dataset
                            else:
                                successful_rivers.append(river)  # Store river object

                            print(f"    📊 Collected river (total: {len(successful_rivers)})")
                        else:
                            print(f"    ❌ Processing failed")
                            failed_scenes.append(scene_id)

                    finally:
                        # Clean up temporary file
                        if os.path.exists(tmp_path):
                            os.unlink(tmp_path)

                except Exception as e:
                    print(f"    ❌ Error processing scene: {str(e)}")
                    failed_scenes.append(f"{year}_{i}")

                # Force garbage collection
                gc.collect()

        except Exception as e:
            print(f"❌ Error accessing collection for {year}: {str(e)}")
            continue

    print(f"\n🎉 Batch processing complete!")
    print(f"✅ Successful: {len(successful_rivers)} rivers")
    print(f"❌ Failed: {len(failed_scenes)} scenes")
    print(f"⏭️  Skipped: {len(skipped_scenes)} scenes")

    if failed_scenes:
        print(f"Failed scenes: {failed_scenes}")

    if skipped_scenes:
        print(f"Skipped scenes: {skipped_scenes}")

    if save_individual:
        print(f"💾 Individual river files saved in: {save_dir}")
        print(f"📁 To load all rivers: rivers = River.load_batch_results('{save_dir}')")

    if download_mndwi:
        print(f"🗺️  MNDWI rasters saved in: {mndwi_dir}")

    return successful_rivers, failed_scenes

load_batch_results(save_dir, min_file_size_kb=0) classmethod

Load all saved rivers from a batch processing directory.

Parameters:

Name	Type	Description	Default
save_dir	str	Directory containing saved river files	required
min_file_size_kb	float	Minimum file size in kilobytes to load (default 0, loads all files)	0

Returns:

Type	Description
list	List of River instances loaded from files

Source code in rivabar/river.py

@classmethod
def load_batch_results(cls, save_dir, min_file_size_kb=0):
    """
    Load all saved rivers from a batch processing directory.

    Parameters
    ----------
    save_dir : str
        Directory containing saved river files
    min_file_size_kb : float, optional
        Minimum file size in kilobytes to load (default 0, loads all files)

    Returns
    -------
    list
        List of River instances loaded from files
    """
    import pickle
    import os
    import glob

    river_files = glob.glob(os.path.join(save_dir, "river_*.pkl"))
    rivers = []
    skipped_files = 0

    print(f"📂 Found {len(river_files)} river files in {save_dir}")
    if min_file_size_kb > 0:
        print(f"🔍 Only loading files larger than {min_file_size_kb} KB")

    for river_file in sorted(river_files):
        # Check file size if minimum size is specified
        if min_file_size_kb > 0:
            file_size_kb = os.path.getsize(river_file) / 1024  # Convert bytes to KB
            if file_size_kb < min_file_size_kb:
                skipped_files += 1
                continue

        try:
            with open(river_file, 'rb') as f:
                river_data = pickle.load(f)

            # Extract river object and restore metadata
            if isinstance(river_data, dict) and 'river' in river_data:
                river = river_data['river']
                # Restore metadata attributes
                river.scene_id = river_data.get('scene_id', 'unknown')
                river.acquisition_date = river_data.get('acquisition_date', 'unknown')
                river.cloud_cover = river_data.get('cloud_cover', 0)
                river.landsat_mission = river_data.get('landsat_mission', 'unknown')
                river.year = river_data.get('year', 0)

                # Create a minimal dataset-like object for CRS info if available
                if river_data.get('dataset_crs') and river._dataset is None:
                    from .data_io import MinimalDataset
                    river._dataset = MinimalDataset(
                        river_data.get('dataset_crs'),
                        river_data.get('dataset_transform'),
                        river_data.get('dataset_shape')
                    )
            else:
                # Backward compatibility - if it's just a river object
                river = river_data

            rivers.append(river)

        except Exception as e:
            print(f"❌ Error loading {river_file}: {e}")

    print(f"✅ Successfully loaded {len(rivers)} rivers")
    if min_file_size_kb > 0:
        print(f"⏭️ Skipped {skipped_files} files smaller than {min_file_size_kb} KB")
    return rivers

batch_process_from_file_list(dirname, fnames, start_x, start_y, end_x, end_y, clear_rasters=True, **processing_kwargs) classmethod

Batch process rivers from a list of existing files.

Parameters:

Name	Type	Description	Default
dirname	str	Directory containing the files	required
fnames	list	List of filenames to process	required
start_x	float	Starting point coordinates (UTM)	required
start_y	float	Starting point coordinates (UTM)	required
end_x	float	Ending point coordinates (UTM)	required
end_y	float	Ending point coordinates (UTM)	required
clear_rasters	bool	Whether to clear raster data after processing (default True)	True
**processing_kwargs	dict	Additional arguments for map_river_banks	{}

Returns:

Name	Type	Description
successful_rivers	list	List of successfully processed River instances
failed_files	list	List of filenames that failed processing

Source code in rivabar/river.py

@classmethod  
def batch_process_from_file_list(cls, dirname, fnames, start_x, start_y, end_x, end_y,
                               clear_rasters=True, **processing_kwargs):
    """
    Batch process rivers from a list of existing files.

    Parameters
    ----------
    dirname : str
        Directory containing the files
    fnames : list
        List of filenames to process
    start_x, start_y : float
        Starting point coordinates (UTM)
    end_x, end_y : float
        Ending point coordinates (UTM)
    clear_rasters : bool, optional
        Whether to clear raster data after processing (default True)
    **processing_kwargs : dict
        Additional arguments for map_river_banks

    Returns
    -------
    successful_rivers : list
        List of successfully processed River instances
    failed_files : list
        List of filenames that failed processing
    """
    import gc

    successful_rivers = []
    failed_files = []

    print(f"🗂️  Batch processing {len(fnames)} files from {dirname}")

    for i, fname in enumerate(fnames):
        print(f"\n--- Processing {i+1}/{len(fnames)}: {fname} ---")

        try:
            # Create River instance
            river = cls(
                fname=fname,
                dirname=dirname,
                start_x=start_x,
                start_y=start_y,
                end_x=end_x,
                end_y=end_y,
                file_type='water_index',
                **processing_kwargs
            )

            # Process the river
            memory_before = river.get_memory_usage()
            success = river.map_river_banks()

            if success:
                print(f"✅ Success! Memory: {memory_before} → {river.get_memory_usage()}")

                # Clear raster data to save memory
                if clear_rasters:
                    river.clear_raster_data()
                    print(f"🧹 Cleared rasters, memory: {river.get_memory_usage()}")

                successful_rivers.append(river)
                print(f"📊 Collected river (total: {len(successful_rivers)})")
            else:
                print(f"❌ Processing failed")
                failed_files.append(fname)

        except Exception as e:
            print(f"❌ Error processing {fname}: {str(e)}")
            failed_files.append(fname)

        # Force garbage collection
        gc.collect()

    print(f"\n🎉 Batch processing complete!")
    print(f"✅ Successful: {len(successful_rivers)} rivers")
    print(f"❌ Failed: {len(failed_files)} files")

    return successful_rivers, failed_files 

save_planetscope_result(save_dir='planetscope_results', scene_id=None, acquisition_date=None, cloud_cover=None, source_files=None, processing_metadata=None)

Save river mapping result for PlanetScope data as a pickle file.

Parameters:

Name	Type	Description	Default
save_dir	str	Directory to save the result file (default: 'planetscope_results')	'planetscope_results'
scene_id	str	Scene identifier. If None, will be extracted from filename or generated	None
acquisition_date	str	Acquisition date in 'YYYY-MM-DD' format. If None, will be extracted from filename	None
cloud_cover	float	Cloud cover percentage	None
source_files	list	List of source files used to create the mosaic	None
processing_metadata	dict	Additional processing metadata	None

Returns:

Type	Description
str	Path to the saved pickle file

Examples:

>>> river.save_planetscope_result(
...     scene_id='rio_beni_20250612',
...     acquisition_date='2025-06-12',
...     cloud_cover=5.2,
...     source_files=['20250612_151136_17_24f4_3B_AnalyticMS_SR_8b_clip.tif']
... )

Source code in rivabar/river.py

def save_planetscope_result(self, save_dir='planetscope_results', scene_id=None, 
                          acquisition_date=None, cloud_cover=None, 
                          source_files=None, processing_metadata=None):
    """
    Save river mapping result for PlanetScope data as a pickle file.

    Parameters
    ----------
    save_dir : str, optional
        Directory to save the result file (default: 'planetscope_results')
    scene_id : str, optional
        Scene identifier. If None, will be extracted from filename or generated
    acquisition_date : str, optional
        Acquisition date in 'YYYY-MM-DD' format. If None, will be extracted from filename
    cloud_cover : float, optional
        Cloud cover percentage
    source_files : list, optional
        List of source files used to create the mosaic
    processing_metadata : dict, optional
        Additional processing metadata

    Returns
    -------
    str
        Path to the saved pickle file

    Examples
    --------
    >>> river.save_planetscope_result(
    ...     scene_id='rio_beni_20250612',
    ...     acquisition_date='2025-06-12',
    ...     cloud_cover=5.2,
    ...     source_files=['20250612_151136_17_24f4_3B_AnalyticMS_SR_8b_clip.tif']
    ... )
    """
    import os
    import pickle
    from datetime import datetime

    # Create save directory if it doesn't exist
    os.makedirs(save_dir, exist_ok=True)

    # Extract metadata from filename if not provided
    if scene_id is None:
        if hasattr(self, 'fname') and self.fname:
            # Extract date from filename like 'ddwi_mosaic_small.tif' or PlanetScope pattern
            fname_base = os.path.splitext(self.fname)[0]
            if 'mosaic' in fname_base.lower():
                scene_id = f"planetscope_mosaic_{datetime.now().strftime('%Y%m%d')}"
            else:
                scene_id = fname_base
        else:
            scene_id = f"planetscope_river_{datetime.now().strftime('%Y%m%d_%H%M%S')}"

    if acquisition_date is None and source_files:
        # Try to extract date from first source file
        for fname in source_files:
            if fname.startswith('20'):  # PlanetScope format: 20250612_...
                try:
                    date_str = fname[:8]  # YYYYMMDD
                    year = int(date_str[:4])
                    month = int(date_str[4:6])
                    day = int(date_str[6:8])
                    acquisition_date = f"{year:04d}-{month:02d}-{day:02d}"
                    break
                except (ValueError, IndexError):
                    continue

    if acquisition_date is None:
        acquisition_date = datetime.now().strftime('%Y-%m-%d')

    # Create filename for saving
    safe_scene_id = scene_id.replace('/', '_').replace('\\', '_')
    river_file = os.path.join(save_dir, f"river_{safe_scene_id}.pkl")

    # Store dataset and MNDWI temporarily and remove them for pickling to save space
    temp_dataset = self._dataset
    temp_mndwi = self._mndwi
    self._dataset = None
    self._mndwi = None

    try:
        # Create results dictionary for saving
        river_data = {
            'river': self,
            'scene_id': scene_id,
            'acquisition_date': acquisition_date,
            'cloud_cover': cloud_cover or 0.0,
            'satellite_mission': 'PlanetScope',
            'processing_timestamp': datetime.now(),
            'source_files': source_files or [],
            'processing_metadata': processing_metadata or {},
            # Store dataset info separately
            'dataset_crs': str(temp_dataset.crs) if temp_dataset else None,
            'dataset_transform': temp_dataset.transform if temp_dataset else None,
            'dataset_shape': temp_dataset.shape if temp_dataset else None,
            # MNDWI info without storing the actual array
            'mndwi_shape': temp_mndwi.shape if temp_mndwi is not None else None,
            'mndwi_dtype': str(temp_mndwi.dtype) if temp_mndwi is not None else None,
            'water_pixel_count': int(np.sum(temp_mndwi)) if temp_mndwi is not None else None,
            # PlanetScope specific metadata
            'data_type': 'planetscope',
            'mosaic_info': {
                'n_tiles': len(source_files) if source_files else 1,
                'water_index_type': getattr(self, 'water_index_type', 'ddwi'),
                'file_type': getattr(self, 'file_type', 'water_index')
            }
        }

        # Add river metadata attributes for consistency with Landsat processing
        self.scene_id = scene_id
        self.acquisition_date = acquisition_date
        self.cloud_cover = cloud_cover or 0.0
        self.satellite_mission = 'PlanetScope'

        # Save the pickle file
        with open(river_file, 'wb') as f:
            pickle.dump(river_data, f)

        # Get file size for reporting
        import os
        file_size_mb = os.path.getsize(river_file) / (1024**2)

        print(f"💾 Saved PlanetScope river result: {river_file}")
        print(f"   Scene ID: {scene_id}")
        print(f"   Acquisition date: {acquisition_date}")
        if source_files:
            print(f"   Source tiles: {len(source_files)}")
        print(f"   File size: {file_size_mb:.1f} MB (MNDWI excluded)")
        if temp_mndwi is not None:
            mndwi_size_mb = temp_mndwi.nbytes / (1024**2)
            print(f"   MNDWI would add: {mndwi_size_mb:.1f} MB")

        return river_file

    except Exception as e:
        print(f"⚠️ Error saving PlanetScope result: {e}")
        raise
    finally:
        # Restore the dataset and MNDWI references
        self._dataset = temp_dataset
        self._mndwi = temp_mndwi

extract_planetscope_metadata_from_files(file_list) staticmethod

Extract metadata from PlanetScope file names.

Parameters:

Name	Type	Description	Default
file_list	list	List of PlanetScope file paths	required

Returns:

Type	Description
dict	Dictionary containing extracted metadata

Source code in rivabar/river.py

@staticmethod
def extract_planetscope_metadata_from_files(file_list):
    """
    Extract metadata from PlanetScope file names.

    Parameters
    ----------
    file_list : list
        List of PlanetScope file paths

    Returns
    -------
    dict
        Dictionary containing extracted metadata
    """
    metadata = {
        'acquisition_dates': [],
        'item_ids': [],
        'unique_dates': set()
    }

    for filepath in file_list:
        filename = os.path.basename(filepath)

        # PlanetScope pattern: YYYYMMDD_HHMMSS_XX_XXXX_3B_AnalyticMS_SR_8b_clip.tif
        if filename.startswith('20') and len(filename) > 8:
            try:
                # Extract date
                date_str = filename[:8]  # YYYYMMDD
                year = int(date_str[:4])
                month = int(date_str[4:6])
                day = int(date_str[6:8])
                acquisition_date = f"{year:04d}-{month:02d}-{day:02d}"

                # Extract time if available
                if len(filename) > 15 and filename[8] == '_':
                    time_str = filename[9:15]  # HHMMSS
                    if time_str.isdigit() and len(time_str) == 6:
                        hour = int(time_str[:2])
                        minute = int(time_str[2:4])
                        second = int(time_str[4:6])
                        full_datetime = f"{acquisition_date}T{hour:02d}:{minute:02d}:{second:02d}"
                    else:
                        full_datetime = acquisition_date
                else:
                    full_datetime = acquisition_date

                metadata['acquisition_dates'].append(full_datetime)
                metadata['unique_dates'].add(acquisition_date)

                # Extract item ID (everything before the first underscore after date_time)
                parts = filename.split('_')
                if len(parts) >= 4:
                    item_id = '_'.join(parts[:4])  # YYYYMMDD_HHMMSS_XX_XXXX
                    metadata['item_ids'].append(item_id)

            except (ValueError, IndexError):
                continue

    # Set primary acquisition date (most common or first)
    if metadata['unique_dates']:
        metadata['primary_acquisition_date'] = sorted(metadata['unique_dates'])[0]
    else:
        from datetime import datetime
        metadata['primary_acquisition_date'] = datetime.now().strftime('%Y-%m-%d')

    return metadata 

collect_river_stats(river, pixel_size=None)

Collect summary statistics from a processed River object.

Parameters:

Name	Type	Description	Default
river	River	A processed River object.	required
pixel_size	float	Pixel size in meters, forwarded to get_channel_widths. Required for rivers whose raster data has been cleared (e.g., loaded from pickles without a saved transform).	None

Returns:

Name	Type	Description
stats	dict	Dictionary with graph sizes (rook and primal), island degree statistics, bank-node degrees, island area/length statistics, channel width statistics, centerline length, and sinuosity. Width/centerline entries are NaN when they cannot be computed.

Source code in rivabar/river.py

def collect_river_stats(river, pixel_size=None):
    """
    Collect summary statistics from a processed River object.

    Parameters
    ----------
    river : River
        A processed River object.
    pixel_size : float, optional
        Pixel size in meters, forwarded to ``get_channel_widths``. Required
        for rivers whose raster data has been cleared (e.g., loaded from
        pickles without a saved transform).

    Returns
    -------
    stats : dict
        Dictionary with graph sizes (rook and primal), island degree
        statistics, bank-node degrees, island area/length statistics,
        channel width statistics, centerline length, and sinuosity.
        Width/centerline entries are NaN when they cannot be computed.
    """
    river._check_processed()

    stats = {}

    G_rook = river._G_rook
    G_primal = river._G_primal

    # Graph sizes
    stats['n_nodes_rook'] = G_rook.number_of_nodes()
    stats['n_edges_rook'] = G_rook.number_of_edges()
    stats['n_nodes_primal'] = G_primal.number_of_nodes()
    stats['n_edges_primal'] = G_primal.number_of_edges()

    # Degree distribution (islands only, skip bank nodes 0 and 1)
    degrees = [G_rook.degree(n) for n in G_rook if n >= 2]
    stats['degree_mean'] = np.mean(degrees) if degrees else np.nan
    stats['degree_median'] = np.median(degrees) if degrees else np.nan
    stats['degree_std'] = np.std(degrees) if degrees else np.nan

    # First two nodes (banks)
    stats['degree_node0'] = G_rook.degree(0) if 0 in G_rook else np.nan
    stats['degree_node1'] = G_rook.degree(1) if 1 in G_rook else np.nan

    # Island areas and lengths
    island_areas = []
    island_lengths = []
    for node in G_rook:
        if node < 2:  # skip the two bank nodes
            continue
        poly = G_rook.nodes[node].get('bank_polygon')
        if poly is None or not hasattr(poly, 'area'):
            # bank_polygon can be missing or an empty placeholder for
            # islands whose polygon could not be built
            continue
        if poly.area > 0:
            island_areas.append(poly.area)
            island_lengths.append(poly.length)

    stats['island_area_mean'] = np.mean(island_areas) if island_areas else np.nan
    stats['island_area_median'] = np.median(island_areas) if island_areas else np.nan
    stats['island_area_std'] = np.std(island_areas) if island_areas else np.nan
    stats['island_length_mean'] = np.mean(island_lengths) if island_lengths else np.nan
    stats['island_length_median'] = np.median(island_lengths) if island_lengths else np.nan
    stats['island_length_std'] = np.std(island_lengths) if island_lengths else np.nan
    stats['n_islands'] = len(island_areas)

    # Channel widths and centerline
    try:
        s, widths = river.get_channel_widths(pixel_size=pixel_size)
        stats['width_mean'] = np.nanmean(widths)
        stats['width_median'] = np.nanmedian(widths)
        stats['width_std'] = np.nanstd(widths)
        stats['centerline_length'] = s[-1]

        cl = river.main_channel_centerline
        if cl is not None:
            straight_line = np.sqrt((cl.coords[-1][0] - cl.coords[0][0])**2 +
                                    (cl.coords[-1][1] - cl.coords[0][1])**2)
            stats['sinuosity'] = cl.length / straight_line if straight_line > 0 else np.nan
        else:
            stats['sinuosity'] = np.nan
    except Exception as e:
        print(f"Width/centerline failed: {e}")
        stats['width_mean'] = np.nan
        stats['width_median'] = np.nan
        stats['width_std'] = np.nan
        stats['centerline_length'] = np.nan
        stats['sinuosity'] = np.nan

    return stats