rivabar.temporal_analysis

rivabar.temporal_analysis

calculate_iou(poly1, poly2)

Calculate the Intersection over Union (IoU) metric between two polygons.

The IoU is defined as the area of the intersection of the two polygons divided by the area of their union. If the area of the union is 0, the function returns 0.

Parameters:

Name	Type	Description	Default
poly1	Polygon	The first polygon.	required
poly2	Polygon	The second polygon.	required

Returns:

Type	Description
float	The IoU of the two polygons. If the area of the union is 0, it returns 0.

Source code in rivabar/temporal_analysis.py

def calculate_iou(poly1, poly2):
    """
    Calculate the Intersection over Union (IoU) metric between two polygons.

    The IoU is defined as the area of the intersection of the two polygons 
    divided by the area of their union. If the area of the union is 0, the 
    function returns 0.

    Parameters
    ----------
    poly1 : shapely.geometry.Polygon
        The first polygon.
    poly2 : shapely.geometry.Polygon
        The second polygon.

    Returns
    -------
    float
        The IoU of the two polygons. If the area of the union is 0, it returns 0.
    """
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", RuntimeWarning)
        intersection = poly1.intersection(poly2).area
    union = poly1.union(poly2).area
    return intersection / union if union != 0 else 0

modified_iou(poly1, poly2)

Calculates a modified version of the Intersection over Union (IoU) metric between two polygons.

Instead of dividing the area of the intersection by the area of the union of the two polygons, it divides the area of the intersection by the area of the smaller polygon. This gives a measure of the fraction of the smaller polygon that is inside the larger polygon.

Parameters:

Name	Type	Description	Default
poly1	Polygon	The first polygon.	required
poly2	Polygon	The second polygon.	required

Returns:

Type	Description
float	The fraction of the smaller polygon that is inside the larger polygon. If the area of the smaller polygon is 0, it returns 0.

Source code in rivabar/temporal_analysis.py

def modified_iou(poly1, poly2):
    """
    Calculates a modified version of the Intersection over Union (IoU) metric between two polygons.

    Instead of dividing the area of the intersection by the area of the union of the two polygons, 
    it divides the area of the intersection by the area of the smaller polygon. This gives a measure 
    of the fraction of the smaller polygon that is inside the larger polygon.

    Parameters
    ----------
    poly1 : shapely.geometry.Polygon
        The first polygon.
    poly2 : shapely.geometry.Polygon
        The second polygon.

    Returns
    -------
    float
        The fraction of the smaller polygon that is inside the larger polygon. If the area of the 
        smaller polygon is 0, it returns 0.
    """
    # Calculate the intersection
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", RuntimeWarning)
        intersection = poly1.intersection(poly2).area
    # Determine the smaller polygon
    smaller_poly_area = min(poly1.area, poly2.area)
    # Calculate the fraction of the smaller polygon that is inside the larger polygon
    return intersection / smaller_poly_area if smaller_poly_area != 0 else 0

cluster_polygons(gdf, iou_threshold, max_days=2 * 365)

Cluster polygons based on Intersection over Union (IoU) and time difference.

Parameters:

Name	Type	Description	Default
gdf	GeoDataFrame	A GeoDataFrame containing polygon geometries and associated attributes.	required
iou_threshold	float	The IoU threshold above which polygons are considered adjacent.	required
max_days	int	The maximum number of days difference allowed between polygons to be considered for clustering (default is 2*365).	2 * 365

Returns:

Name	Type	Description
G	Graph	A graph where nodes represent polygons and edges represent adjacency based on IoU and time difference.
clusters	list of sets	A list of sets, where each set contains the indices of polygons that form a cluster.

Source code in rivabar/temporal_analysis.py

def cluster_polygons(gdf, iou_threshold, max_days=2*365):
    """
    Cluster polygons based on Intersection over Union (IoU) and time difference.

    Parameters
    ----------
    gdf : GeoDataFrame
        A GeoDataFrame containing polygon geometries and associated attributes.
    iou_threshold : float
        The IoU threshold above which polygons are considered adjacent.
    max_days : int, optional
        The maximum number of days difference allowed between polygons to be considered for clustering (default is 2*365).

    Returns
    -------
    G : networkx.Graph
        A graph where nodes represent polygons and edges represent adjacency based on IoU and time difference.
    clusters : list of sets
        A list of sets, where each set contains the indices of polygons that form a cluster.
    """
    # Create a spatial index for the polygons
    sindex = gdf.sindex
    # Create a graph to represent adjacency (IoU above threshold)
    G = nx.Graph()
    for i, row1 in tqdm(gdf.iterrows()):
        if row1['type'] != 1 and row1['type'] != 0: # ignore the main banks
            poly1 = row1.geometry
            n_days1 = row1.n_days
            # Possible matches index
            with warnings.catch_warnings():
                warnings.simplefilter("ignore", RuntimeWarning)
                possible_matches_index = list(sindex.intersection(poly1.bounds))
                possible_matches = gdf.iloc[possible_matches_index]
                precise_matches = possible_matches[possible_matches.geometry.intersects(poly1)]        
            for j, row2 in precise_matches.iterrows():
                if i < j:  # Avoid duplicate pairs
                    poly2 = row2.geometry
                    n_days2 = row2.n_days
                    if np.abs(n_days2-n_days1) <= max_days: # time difference less than 'max_days'
                        iou = modified_iou(poly1, poly2)
                        if iou > iou_threshold:
                            G.add_edge(i, j)
    # Find clusters (connected components) in the graph
    clusters = list(nx.connected_components(G))
    return G, clusters

get_ch_and_bar_areas(gdf, xmin, xmax, ymin, ymax)

Calculate channel and bar areas within a specified area of interest (AOI) over time.

Parameters:

Name	Type	Description	Default
gdf	GeoDataFrame	A GeoDataFrame containing geometries and attributes of river banks and bars.	required
xmin	float	Minimum x-coordinate of the AOI.	required
xmax	float	Maximum x-coordinate of the AOI.	required
ymin	float	Minimum y-coordinate of the AOI.	required
ymax	float	Maximum y-coordinate of the AOI.	required

Returns:

Name	Type	Description
dates	list	List of unique dates corresponding to the n_days in the GeoDataFrame.
all_bars	list	List of geometries representing all bars within the AOI for each date.
chs	list	List of geometries representing channels within the AOI for each date.
ch_belts	list	List of geometries representing channel belts within the AOI for each date.
bar_areas	list	List of areas of bars within the AOI for each date.
ch_areas	list	List of areas of channels within the AOI for each date.

Source code in rivabar/temporal_analysis.py

def get_ch_and_bar_areas(gdf, xmin, xmax, ymin, ymax):
    """
    Calculate channel and bar areas within a specified area of interest (AOI) over time.

    Parameters
    ----------
    gdf : GeoDataFrame
        A GeoDataFrame containing geometries and attributes of river banks and bars.
    xmin : float
        Minimum x-coordinate of the AOI.
    xmax : float
        Maximum x-coordinate of the AOI.
    ymin : float
        Minimum y-coordinate of the AOI.
    ymax : float
        Maximum y-coordinate of the AOI.

    Returns
    -------
    dates : list
        List of unique dates corresponding to the n_days in the GeoDataFrame.
    all_bars : list
        List of geometries representing all bars within the AOI for each date.
    chs : list
        List of geometries representing channels within the AOI for each date.
    ch_belts : list
        List of geometries representing channel belts within the AOI for each date.
    bar_areas : list
        List of areas of bars within the AOI for each date.
    ch_areas : list
        List of areas of channels within the AOI for each date.
    """
    ch_areas = []
    bar_areas = []
    chs = []
    ch_belts = []
    all_bars = []
    dates = []
    aoi_poly = Polygon([(xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax)])
    for n_days in tqdm(gdf.n_days.unique()):
        date = gdf[(gdf['n_days']== n_days) & (gdf['type']==0)].date.dt.year
        bank0 = gdf[(gdf['n_days']== n_days) & (gdf['type']==0)].geometry.values[0]
        bank0 = bank0.intersection(aoi_poly)
        bank1 = gdf[(gdf['n_days']== n_days) & (gdf['type']==1)].geometry.values[0]
        bank1 = bank1.intersection(aoi_poly)
        ch_belt = aoi_poly.difference(unary_union([bank0, bank1]))
        bars = []
        for i in gdf.index:
            if gdf.n_days[i] == n_days and gdf['type'][i] != 0 and gdf['type'][i] != 1:
                bars.append(gdf.geometry[i])
        bars = unary_union(bars)
        bars = bars.intersection(aoi_poly)
        bar_areas.append(bars.area)
        all_bars.append(bars)
        ch = ch_belt.difference(bars)
        ch_areas.append(ch.area)
        chs.append(ch)
        ch_belts.append(ch_belt)
        dates.append(date)
    return dates, all_bars, chs, ch_belts, bar_areas, ch_areas

create_and_plot_bars(rivers, ts1, ts2, ax1=None, ax2=None, depo_cmap='Blues', erosion_cmap='Reds', alpha=0.5, aoi=None, colorbar=True, color_scale_timestep=None)

Create preserved scroll bar polygons and erosion polygons from a list of rook neighborhood graphs and plot them. It also handles the creation of the color map for the plot.

Parameters:

Name	Type	Description	Default
rivers	list	A list of river objects with rook neighborhood graphs.	required
ts1	int	The first time step to consider (inclusive).	required
ts2	int	The last time step to consider (inclusive).	required
ax1	Axes	The axes on which to plot the deposition polygons. Defaults to None.	None
ax2	Axes	The axes on which to plot the erosion polygons. Defaults to None.	None
depo_cmap	str	The name of the color map to use for the deposition polygons. Defaults to "Blues".	'Blues'
erosion_cmap	str	The name of the color map to use for the erosion polygons. Defaults to "Reds".	'Reds'
aoi	list or tuple	Area of interest defined as [xmin, xmax, ymin, ymax]. If provided, all channel polygons will be cropped to this area before processing. Defaults to None.	None
colorbar	bool	Whether to show colorbar. Defaults to True.	True
color_scale_timestep	int	Alternative final timestep (absolute index into rivers) to use for color scaling instead of ts2, so that colors stay consistent across plots/animation frames with different ts2 values. The actual analysis still uses the original ts1 and ts2. Defaults to None.	None

Returns:

Name	Type	Description
chs	list	A list of channel polygons (cropped to AOI if provided).
bars	list	A list of scroll bar polygons.
erosions_final	list	A list of final erosion polygons.
aoi_dates	list	A list of datetime objects corresponding to the channels that remain after AOI cropping. If no AOI is provided, this equals the original dates list.
aoi_centerlines	list	A list of main channel centerline segments (LineStrings) corresponding to the AOI. If no AOI is provided, returns all available centerlines.

Example

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20,15), sharex=True, sharey=True) chs, bars, erosions, dates, centerlines = rb.create_and_plot_bars(rivers, 0, len(rivers)-1, ax1=ax1, ax2=ax2)

With area of interest:

aoi = [xmin, xmax, ymin, ymax] chs, bars, erosions, aoi_dates, aoi_centerlines = rb.create_and_plot_bars(rivers, 5, 15, ax1=ax1, ax2=ax2, aoi=aoi)

Source code in rivabar/temporal_analysis.py

def create_and_plot_bars(rivers, ts1, ts2, ax1=None, ax2=None, depo_cmap="Blues", erosion_cmap="Reds", alpha=0.5, aoi=None, colorbar=True, color_scale_timestep=None):
    """
    Create preserved scroll bar polygons and erosion polygons from a list of rook neighborhood graphs 
    and plot them. It also handles the creation of the color map for the plot.

    Parameters
    ----------
    rivers : list
        A list of river objects with rook neighborhood graphs.
    ts1 : int
        The first time step to consider (inclusive).
    ts2 : int
        The last time step to consider (inclusive).
    ax1 : matplotlib.axes.Axes, optional
        The axes on which to plot the deposition polygons. Defaults to None.
    ax2 : matplotlib.axes.Axes, optional
        The axes on which to plot the erosion polygons. Defaults to None.
    depo_cmap : str, optional
        The name of the color map to use for the deposition polygons. Defaults to "Blues".
    erosion_cmap : str, optional
        The name of the color map to use for the erosion polygons. Defaults to "Reds".
    aoi : list or tuple, optional
        Area of interest defined as [xmin, xmax, ymin, ymax]. If provided, all channel polygons
        will be cropped to this area before processing. Defaults to None.
    colorbar : bool, optional
        Whether to show colorbar. Defaults to True.
    color_scale_timestep : int, optional
        Alternative final timestep (absolute index into ``rivers``) to use for
        color scaling instead of ts2, so that colors stay consistent across
        plots/animation frames with different ts2 values. The actual analysis
        still uses the original ts1 and ts2. Defaults to None.

    Returns
    -------
    chs : list
        A list of channel polygons (cropped to AOI if provided).
    bars : list
        A list of scroll bar polygons.
    erosions_final : list
        A list of final erosion polygons.
    aoi_dates : list
        A list of datetime objects corresponding to the channels that remain after AOI cropping.
        If no AOI is provided, this equals the original dates list.
    aoi_centerlines : list
        A list of main channel centerline segments (LineStrings) corresponding to the AOI.
        If no AOI is provided, returns all available centerlines.

    Example
    -------
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20,15), sharex=True, sharey=True)
    chs, bars, erosions, dates, centerlines = rb.create_and_plot_bars(rivers, 0, len(rivers)-1, ax1=ax1, ax2=ax2)

    # With area of interest:
    aoi = [xmin, xmax, ymin, ymax]
    chs, bars, erosions, aoi_dates, aoi_centerlines = rb.create_and_plot_bars(rivers, 5, 15, ax1=ax1, ax2=ax2, aoi=aoi)
    """
    # Validate input parameters
    if ts1 < 0:
        print('ts1 must be >= 0!')
        return None, None, None, None, None
    if ts2 <= ts1:
        print('ts2 must be greater than ts1!')
        return None, None, None, None, None
    if ts2 > len(rivers) - 1:
        print('ts2 must be <= len(rivers) - 1!')
        return None, None, None, None, None
    if ts2 - ts1 < 2:
        print('Need at least 2 time steps (ts2 - ts1 >= 2)!')
        return None, None, None, None, None

    # Calculate the number of time steps to process
    ts = ts2 - ts1 + 1

    # Calculate alternative timestep for color scaling if provided
    if color_scale_timestep is not None:
        color_ts = color_scale_timestep - ts1 + 1
    else:
        color_ts = ts

    # Create AOI polygon if provided
    aoi_polygon = None
    if aoi is not None:
        from shapely.geometry import box
        xmin, xmax, ymin, ymax = aoi
        aoi_polygon = box(xmin, ymin, xmax, ymax)
        print(f"Using AOI: xmin={xmin}, xmax={xmax}, ymin={ymin}, ymax={ymax}")

    dates = []
    for river in rivers[ts1:ts2+1]:
        dates.append(datetime.strptime(river.acquisition_date, "%Y-%m-%d"))
    df = pd.DataFrame(dates, columns=['date'])
    # Determine the range of dates
    min_date = df['date'].min()
    max_date = df['date'].max()
    # Generate a list of January 1st dates within the range
    jan_first_dates = pd.date_range(start=min_date, end=max_date, freq='AS')
    # Create a dataframe for the January 1st dates
    jan_first_df = pd.DataFrame(jan_first_dates, columns=['date'])
    # Concatenate the original dataframe with the January 1st dataframe
    combined_df = pd.concat([df, jan_first_df], ignore_index=True)
    combined_df['date_type'] = 0
    combined_df.iloc[len(df):, 1] = 1
    # Sort by datetime column
    df = combined_df.sort_values(by='date').reset_index(drop=True)
    df['timedelta'] = df['date'] - df['date'].min()
    df['n_days'] = df['timedelta'].dt.days
    n_days = df.n_days.values

    bars = [] # these are 'scroll' bars - shapely MultiPolygon objects that correspond to one time step
    chs = [] # list of channels - shapely Polygon objects
    all_chs = [] # list of merged channels (to be used for erosion)
    erosions = []
    erosion_all_chs = []
    cmap = plt.get_cmap(depo_cmap)

    # create list of channels:
    for i in trange(ts-1):
        ch1 = create_channel_nw_polygon(rivers[ts1+i]._G_rook, buffer=10, dataset=rivers[ts1+i]._dataset)
        ch1 = ch1.buffer(0)
        if type(ch1) == MultiPolygon:
            areas = []
            for geom in ch1.geoms:
                areas.append(geom.area)
            ch1 = ch1.geoms[np.argmax(areas)]
        # print(ts1+i+1)
        ch2 = create_channel_nw_polygon(rivers[ts1+i+1]._G_rook, buffer=10, dataset=rivers[ts1+i+1]._dataset)
        ch2 = ch2.buffer(0)
        if type(ch2) == MultiPolygon:
            areas = []
            for geom in ch2.geoms:
                areas.append(geom.area)
            ch2 = ch2.geoms[np.argmax(areas)]
        chs.append(ch1)
    chs.append(ch2) # append last channel

    # Extract centerlines from rivers
    centerlines = []
    for i, river in enumerate(rivers[ts1:ts2+1]):
        try:
            centerline = river.main_channel_centerline
            if centerline is not None:
                centerlines.append(centerline)
            else:
                print(f"Warning: No centerline available for river {ts1+i}")
                centerlines.append(None)
        except AttributeError:
            print(f"Warning: River {ts1+i} does not have main_channel_centerline attribute")
            centerlines.append(None)

    # Crop channels and centerlines to AOI if provided
    if aoi_polygon is not None:
        print(f"Cropping {len(chs)} channel polygons and {len(centerlines)} centerlines to AOI...")
        cropped_chs = []
        cropped_centerlines = []

        for i, (ch, centerline) in enumerate(zip(chs, centerlines)):
            try:
                # Crop channel polygon
                cropped_ch = ch.intersection(aoi_polygon)

                # Handle different geometry types returned by intersection
                if cropped_ch.is_empty:
                    print(f"Warning: Channel {ts1+i} has no intersection with AOI")
                    cropped_chs.append(None)
                elif type(cropped_ch) == MultiPolygon:
                    # Keep the largest polygon if multiple are returned
                    areas = [geom.area for geom in cropped_ch.geoms if type(geom) == Polygon]
                    if areas:
                        largest_idx = np.argmax(areas)
                        largest_poly = [geom for geom in cropped_ch.geoms if type(geom) == Polygon][largest_idx]
                        cropped_chs.append(largest_poly)
                    else:
                        print(f"Warning: Channel {ts1+i} intersection contains no valid polygons")
                        cropped_chs.append(None)
                elif type(cropped_ch) == Polygon:
                    cropped_chs.append(cropped_ch)
                else:
                    print(f"Warning: Channel {ts1+i} intersection returned unexpected geometry type: {type(cropped_ch)}")
                    cropped_chs.append(None)

                # Crop centerline
                if centerline is not None:
                    try:
                        cropped_centerline = centerline.intersection(aoi_polygon)
                        if cropped_centerline.is_empty:
                            cropped_centerlines.append(None)
                        elif hasattr(cropped_centerline, 'geoms'):
                            # MultiLineString - take the longest segment
                            segments = [geom for geom in cropped_centerline.geoms if hasattr(geom, 'length')]
                            if segments:
                                longest_segment = max(segments, key=lambda x: x.length)
                                cropped_centerlines.append(longest_segment)
                            else:
                                cropped_centerlines.append(None)
                        else:
                            # Single LineString
                            cropped_centerlines.append(cropped_centerline)
                    except Exception as e:
                        print(f"Warning: Error cropping centerline {ts1+i}: {e}")
                        cropped_centerlines.append(None)
                else:
                    cropped_centerlines.append(None)

            except Exception as e:
                print(f"Error cropping channel {ts1+i}: {e}")
                cropped_chs.append(None)
                cropped_centerlines.append(None)

        # Filter out None values and update lists
        valid_indices = [i for i, ch in enumerate(cropped_chs) if ch is not None]
        chs = [cropped_chs[i] for i in valid_indices]
        aoi_centerlines = [cropped_centerlines[i] for i in valid_indices]

        # Update dates to match the valid channels
        aoi_dates = [dates[i] for i in valid_indices]

        # Update ts to reflect the number of valid cropped channels
        original_ts = ts
        ts = len(chs)

        if ts < 2:
            print(f"Error: Only {ts} valid channels after cropping to AOI. Need at least 2 channels for analysis.")
            return None, None, None, None, None

        if ts != original_ts:
            print(f"Note: Using {ts} channels after AOI cropping (originally {original_ts})")

        # Print centerline summary
        valid_centerlines = [cl for cl in aoi_centerlines if cl is not None]
        print(f"Note: {len(valid_centerlines)} valid centerlines in AOI (out of {len(aoi_centerlines)})")
    else:
        # If no AOI provided, use original dates and centerlines
        aoi_dates = dates
        aoi_centerlines = centerlines

    # create list of merged channels:
    all_ch = chs[ts-1]
    all_chs.append(all_ch)
    for i in trange(2, ts): 
        all_ch = all_ch.union(chs[ts-i])
        all_chs.append(all_ch)
    erosion_all_ch = chs[0]
    erosion_all_chs.append(erosion_all_ch)
    for i in trange(1, ts-1): 
        erosion_all_ch = erosion_all_ch.union(chs[i])
        erosion_all_chs.append(erosion_all_ch)

    # create scroll bars and plot them:
    for i in trange(ts-1): # create scroll bars
        bar = chs[i].difference(all_chs[ts-i-2]) # scroll bar defined by difference
        bars.append(bar)
        erosion = chs[i+1].difference(erosion_all_chs[i])
        erosions.append(erosion)
        if ax1: # plotting
            color = cmap(i/float(color_ts))
            if type(bar) == MultiPolygon or type(bar) == GeometryCollection:
                for b in bar.geoms:
                    if type(b) == Polygon:
                        if len(b.interiors) == 0:
                            ax1.fill(b.exterior.xy[0], b.exterior.xy[1], facecolor=color, edgecolor='k', linewidth=0.2, alpha=alpha)
                        else:
                            plot_polygon(ax1, b, facecolor=color, edgecolor='k', linewidth=0.2, alpha=alpha)
            if type(bar) == Polygon:
                if len(bar.interiors) == 0:
                    ax1.fill(bar.exterior.xy[0], bar.exterior.xy[1], facecolor=color, edgecolor='k', linewidth=0.2, alpha=alpha)
                else:
                    plot_polygon(ax1, bar, facecolor=color, edgecolor='k', linewidth=0.2, alpha=alpha)
    if ax1:
        # Use the last channel (accounting for potential AOI cropping)
        final_channel = chs[-1] if chs else ch2
        plot_polygon(ax1, final_channel, facecolor='lightblue', edgecolor='k', alpha=alpha)
        ax1.set_aspect('equal')
        # date = datetime.strptime(rivers[min(ts2, len(rivers)-1)].acquisition_date, "%Y-%m-%d")
        # ax1.set_title('Deposition, ' + date.strftime('%m')+'/'+date.strftime('%d')+'/'+date.strftime('%Y'))
        if colorbar:
            xlim = ax1.get_xlim()
            ylim = ax1.get_ylim()
            dummy_data = np.vstack((n_days, n_days))
            dummy_im = ax1.imshow(dummy_data, cmap=cmap)
            dummy_im.remove() # remove the image
            ax1.set_xlim(xlim)
            ax1.set_ylim(ylim)
            cbar = plt.colorbar(dummy_im, ax=ax1, shrink=0.5)
            cbar.set_label('age', fontsize=12)
            # Define the tick values and labels you want
            tick_values = df.n_days[df.date_type == 1].values
            tick_labels = df.date[df.date_type == 1].dt.year.values
            # Set the ticks and tick labels on the colorbar
            cbar.set_ticks(tick_values)
            cbar.set_ticklabels(tick_labels)

    # compute final erosional polygons:
    erosions_final = []
    cmap = plt.get_cmap(erosion_cmap)
    for i in trange(ts-1):
        color = cmap(i/float(color_ts))
        if i==0:
            all_bars = chs[i]
        else:
            all_bars = all_bars.union(bars[i]) # all depositional area up until this point in time
        # erosion = erosions[i].difference(all_bars) # remove net depositional areas from final erosion
        erosion = erosions[i].buffer(0).difference(all_bars.buffer(0))
        erosions_final.append(erosion)
        if ax2: # plotting
            if type(erosion) == MultiPolygon or type(erosion) == GeometryCollection:
                for b in erosion.geoms:
                    if type(b) == Polygon:
                        if len(b.interiors) == 0:
                            ax2.fill(b.exterior.xy[0], b.exterior.xy[1], facecolor=color, edgecolor='k', linewidth=0.2, alpha=alpha)
                            # ax2.fill(b.exterior.xy[0], b.exterior.xy[1], facecolor='xkcd:faded red', edgecolor='k', linewidth=0.2, alpha=alpha)
                        else:
                            plot_polygon(ax2, b, facecolor=color, edgecolor='k', linewidth=0.2, alpha=alpha)
                            # plot_polygon(ax2, b, facecolor='xkcd:faded red', edgecolor='k', linewidth=0.2, alpha=alpha)
            if type(erosion) == Polygon:
                if len(erosion.interiors) == 0:
                    ax2.fill(erosion.exterior.xy[0], erosion.exterior.xy[1], facecolor=color, edgecolor='k', linewidth=0.2, alpha=alpha)
                    # ax2.fill(erosion.exterior.xy[0], erosion.exterior.xy[1], facecolor='xkcd:faded red', edgecolor='k', linewidth=0.2, alpha=alpha)
                else:
                    plot_polygon(ax2, erosion, facecolor=color, edgecolor='k', linewidth=0.2, alpha=alpha)
                    # plot_polygon(ax2, erosion, facecolor='xkcd:faded red', edgecolor='k', linewidth=0.2, alpha=alpha)
    if ax2:
        plot_polygon(ax2, chs[0], facecolor='lightblue', edgecolor='k', alpha=alpha)
        ax2.set_aspect('equal')
        # date = datetime.strptime(rivers[min(ts2, len(rivers)-1)].acquisition_date, "%Y-%m-%d")
        # ax2.set_title('Erosion, ' + date.strftime('%m')+'/'+date.strftime('%d')+'/'+date.strftime('%Y'))
        if colorbar:
            xlim = ax2.get_xlim()
            ylim = ax2.get_ylim()
            dummy_data = np.vstack((n_days, n_days))
            dummy_im = ax2.imshow(dummy_data, cmap=cmap)
            dummy_im.remove() # remove the image
            ax2.set_xlim(xlim)
            ax2.set_ylim(ylim)
            cbar = plt.colorbar(dummy_im, ax=ax2, shrink=0.5)
            cbar.set_label('age', fontsize=12)
            # Define the tick values and labels you want
            tick_values = df.n_days[df.date_type == 1].values
            tick_labels = df.date[df.date_type == 1].dt.year.values
            # Set the ticks and tick labels on the colorbar
            cbar.set_ticks(tick_values)
            cbar.set_ticklabels(tick_labels)
        # plt.tight_layout()
    return chs, bars, erosions_final, aoi_dates, aoi_centerlines

create_geodataframe_from_bank_polygons(G_rooks, crs)

Creates a GeoDataFrame from bank polygons.

Parameters:

Name	Type	Description	Default
G_rooks	list	A list of graph objects, each containing nodes with 'bank_polygon' attributes.	required
crs	str	Coordinate reference system in EPSG code format (e.g., '4326').	required

Returns:

Name	Type	Description
gdf	GeoDataFrame	A GeoDataFrame containing the bank polygons with additional attributes: - 'geometry': The bank polygons. - 'date': The date extracted from the graph object names. - 'timedelta': The time difference from the earliest date. - 'n_days': The number of days since the earliest date. - 'length': The length of each polygon. - 'type': The type of bank (0, 1, or 2).

Source code in rivabar/temporal_analysis.py

def create_geodataframe_from_bank_polygons(G_rooks, crs):
    """
    Creates a GeoDataFrame from bank polygons.

    Parameters
    ----------
    G_rooks : list
        A list of graph objects, each containing nodes with 'bank_polygon' attributes.
    crs : str
        Coordinate reference system in EPSG code format (e.g., '4326').

    Returns
    -------
    gdf : geopandas.GeoDataFrame
        A GeoDataFrame containing the bank polygons with additional attributes:
        - 'geometry': The bank polygons.
        - 'date': The date extracted from the graph object names.
        - 'timedelta': The time difference from the earliest date.
        - 'n_days': The number of days since the earliest date.
        - 'length': The length of each polygon.
        - 'type': The type of bank (0, 1, or 2).
    """
    bank_polys = []
    bank_type = []
    dates = []
    for G_rook in G_rooks:
        dates.append(datetime.strptime(G_rook.name[12:20], "%Y%m%d"))
    date = []
    for j in range(len(G_rooks)):
        if len(G_rook) > 0:
            for i in range(0, len(G_rooks[j])):
                bank_polys.append(G_rooks[j].nodes()[i]['bank_polygon'])
                date.append(dates[j])           
                if i == 0:
                    bank_type.append(0)
                elif i == 1:
                    bank_type.append(1)
                else:
                    bank_type.append(2)
    gdf = gpd.GeoDataFrame(bank_polys, columns = ['geometry'])
    gdf['date'] = date
    gdf['timedelta'] = gdf['date'] - gdf['date'].min()
    gdf['n_days'] = gdf['timedelta'].dt.days
    gdf['length'] = gdf.length
    gdf['type'] = bank_type
    gdf.crs = 'epsg:'+crs
    gdf = gdf.sort_values(by='date')
    gdf = gdf.reset_index(drop=True)
    return gdf

create_dataframe_from_bank_polygons(rivers)

Creates a GeoDataFrame from bank polygons extracted from river objects.

Parameters:

Name	Type	Description	Default
rivers	list	A list of river objects, each containing a graph (_G_rook) with nodes that have 'bank_polygon' attributes and an acquisition_date property.	required

Returns:

Name	Type	Description
gdf	GeoDataFrame	A GeoDataFrame containing the bank polygons with additional attributes: - 'geometry': The bank polygons. - 'year', 'month', 'day': Individual date components. - 'date': The full date as datetime object. - 'timedelta': The time difference from the earliest date. - 'n_days': The number of days since the earliest date. - 'length': The length of each polygon. - 'type': The type of bank (0, 1, or 2).

Notes

Invalid geometries are fixed using the buffer(0) technique.

Source code in rivabar/temporal_analysis.py

def create_dataframe_from_bank_polygons(rivers):
    """
    Creates a GeoDataFrame from bank polygons extracted from river objects.

    Parameters
    ----------
    rivers : list
        A list of river objects, each containing a graph (_G_rook) with nodes 
        that have 'bank_polygon' attributes and an acquisition_date property.

    Returns
    -------
    gdf : geopandas.GeoDataFrame
        A GeoDataFrame containing the bank polygons with additional attributes:
        - 'geometry': The bank polygons.
        - 'year', 'month', 'day': Individual date components.
        - 'date': The full date as datetime object.
        - 'timedelta': The time difference from the earliest date.
        - 'n_days': The number of days since the earliest date.
        - 'length': The length of each polygon.
        - 'type': The type of bank (0, 1, or 2).

    Notes
    -----
    Invalid geometries are fixed using the buffer(0) technique.
    """
    bank_polys = []
    year = []
    month = []
    day = []
    bank_type = []
    for j in range(len(rivers)):
        for i in range(0, len(rivers[j]._G_rook)):
            bank_polys.append(rivers[j]._G_rook.nodes()[i]['bank_polygon'])
            year.append(int(rivers[j].acquisition_date[:4]))
            month.append(int(rivers[j].acquisition_date[5:7]))
            day.append(int(int(rivers[j].acquisition_date[8:])))             
            if i == 0:
                bank_type.append(0)
            elif i == 1:
                bank_type.append(1)
            else:
                bank_type.append(2)
    gdf = geopandas.GeoDataFrame(bank_polys, columns = ['geometry'])
    gdf['year'] = year
    gdf['month'] = month
    gdf['day'] = day
    gdf['date'] = pd.to_datetime(gdf[['year', 'month', 'day']])
    gdf['timedelta'] = gdf['date'] - gdf['date'].min()
    gdf['n_days'] = gdf['timedelta'].dt.days
    gdf['length'] = gdf.length
    gdf['type'] = bank_type
    gdf.crs = 'epsg:'+str(rivers[0]._dataset.crs.to_epsg())
    gdf = gdf.sort_values(by='date')
    gdf = gdf.reset_index(drop=True)
    for i in range(0, len(gdf.geometry)):
        if not gdf.geometry[i].is_valid:
            gdf.geometry[i] = gdf.geometry[i].buffer(0)
    return gdf

get_landsat_scene_crs(path_number, row_number, year=2020)

Get the CRS used by Landsat scenes for a specific path/row.

Requires the optional earthengine-api package and an authenticated Earth Engine session.

Parameters:

Name	Type	Description	Default
path_number	int	WRS path and row	required
row_number	int	WRS path and row	required
year	int	Year to sample (default 2020)	2020

Returns:

Name	Type	Description
crs_string	str	EPSG code or CRS string used by Landsat scenes

Source code in rivabar/temporal_analysis.py

def get_landsat_scene_crs(path_number, row_number, year=2020):
    """
    Get the CRS used by Landsat scenes for a specific path/row.

    Requires the optional earthengine-api package and an authenticated
    Earth Engine session.

    Parameters
    ----------
    path_number, row_number : int
        WRS path and row
    year : int, optional
        Year to sample (default 2020)

    Returns
    -------
    crs_string : str
        EPSG code or CRS string used by Landsat scenes
    """
    import ee

    # Get a sample scene to extract CRS
    if year < 2013:
        collection_id = "LANDSAT/LT05/C02/T1_L2"
    elif year < 2021:
        collection_id = "LANDSAT/LC08/C02/T1_L2"
    else:
        collection_id = "LANDSAT/LC09/C02/T1_L2"

    collection = (ee.ImageCollection(collection_id)
                  .filterDate(f'{year}-01-01', f'{year}-12-31')
                  .filter(ee.Filter.eq('WRS_PATH', path_number))
                  .filter(ee.Filter.eq('WRS_ROW', row_number))
                  .first())

    # Get the CRS
    crs = collection.projection().crs().getInfo()
    return crs

convert_to_landsat_crs(point0_lonlat, point1_lonlat, path_number, row_number)

Convert coordinates to match the Landsat scene's CRS.

Source code in rivabar/temporal_analysis.py

def convert_to_landsat_crs(point0_lonlat, point1_lonlat, path_number, row_number):
    """
    Convert coordinates to match the Landsat scene's CRS.
    """

    # Get the Landsat scene CRS
    landsat_crs = get_landsat_scene_crs(path_number, row_number)
    print(f"Landsat scene CRS: {landsat_crs}")

    # Create transformer
    transformer = Transformer.from_crs("EPSG:4326", landsat_crs, always_xy=True)

    # Convert coordinates
    start_x, start_y = transformer.transform(point0_lonlat[0], point0_lonlat[1])
    end_x, end_y = transformer.transform(point1_lonlat[0], point1_lonlat[1])

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

    return start_x, start_y, end_x, end_y, landsat_crs

collect_river_endpoints(m)

Interactive map tool to collect start and end points for river analysis.

Source code in rivabar/temporal_analysis.py

def collect_river_endpoints(m):
    """Interactive map tool to collect start and end points for river analysis."""
    import ipywidgets as widgets
    from IPython.display import display

    # Data storage
    endpoints = {'start': None, 'end': None}

    # Status displays
    status = widgets.HTML(value="<b>Click to set river START point</b>")
    start_point = widgets.Text(description='Start point:')
    end_point = widgets.Text(description='End point:')

    # Create button to confirm points
    confirm_button = widgets.Button(
        description='Confirm Points',
        disabled=True,
        button_style='success'
    )

    # Create button to reset points
    reset_button = widgets.Button(
        description='Reset',
        button_style='warning'
    )

    def handle_click(**kwargs):
        if kwargs.get('type') == 'click':
            coords = kwargs.get('coordinates')
            lon, lat = round(coords[0], 6), round(coords[1], 6)

            # Add marker based on which point we're capturing
            if endpoints['start'] is None:
                # Mark start point
                endpoints['start'] = [lon, lat]
                m.add_marker(location=[lat, lon]) #, popup=widgets.HTML("START"))
                start_point.value = f"[{lon}, {lat}]"
                status.value = "<b>Click to set river END point</b>"

            elif endpoints['end'] is None:
                # Mark end point
                endpoints['end'] = [lon, lat]
                m.add_marker(location=[lat, lon]) #, popup=widgets.HTML("END"))
                end_point.value = f"[{lon}, {lat}]"
                status.value = "<b>Points captured! Click 'Confirm' to use these points.</b>"
                confirm_button.disabled = False

    def handle_reset(b):
        """Reset the map and points."""
        endpoints['start'] = None
        endpoints['end'] = None
        start_point.value = ""
        end_point.value = ""
        status.value = "<b>Click to set river START point</b>"
        confirm_button.disabled = True
        # Reset map
        # m.clear_markers()

    def handle_confirm(b):
        """Return the confirmed endpoints."""
        status.value = "<b>✅ Points confirmed! Ready to use in your analysis.</b>"

        # Format points as code snippet for easy copy-paste
        point_code = f"""
        # River endpoints (longitude, latitude)
        point0 = {endpoints['start']}  # Start point
        point1 = {endpoints['end']}  # End point
        """

        # Display copy-pasteable code
        display(widgets.HTML(f"<pre>{point_code}</pre>"))

    # Attach event handlers
    m.on_interaction(handle_click)
    reset_button.on_click(handle_reset)
    confirm_button.on_click(handle_confirm)

    # Display map and controls
    display(m)
    display(widgets.VBox([
        status,
        widgets.HBox([start_point, end_point]),
        widgets.HBox([confirm_button, reset_button])
    ]))

    return endpoints

plot_deposition_erosion_with_dates(bars, erosions, dates, centerlines, figsize=(15, 8))

Plot deposition and erosion with bar widths proportional to time intervals. Rates are normalized by centerline length.

Parameters:

Name	Type	Description	Default
bars	list	List of deposition polygons for each time step (length n-1)	required
erosions	list	List of erosion polygons for each time step (length n-1)	required
dates	list	List of acquisition dates (length n)	required
centerlines	list	List of centerline geometries (LineString or MultiLineString) for each time step (length n)	required
figsize	tuple	Figure size (width, height)	(15, 8)

Source code in rivabar/temporal_analysis.py

def plot_deposition_erosion_with_dates(bars, erosions, dates, centerlines, figsize=(15, 8)):
    """
    Plot deposition and erosion with bar widths proportional to time intervals.
    Rates are normalized by centerline length.

    Parameters
    ----------
    bars : list
        List of deposition polygons for each time step (length n-1)
    erosions : list  
        List of erosion polygons for each time step (length n-1)
    dates : list
        List of acquisition dates (length n)
    centerlines : list
        List of centerline geometries (LineString or MultiLineString) for each time step (length n)
    figsize : tuple
        Figure size (width, height)
    """
    import pandas as pd
    import matplotlib.dates as mdates
    from shapely.geometry import LineString, MultiLineString

    # Convert dates to pandas datetime if they aren't already
    dates_pd = pd.to_datetime(dates)

    # Calculate time intervals between consecutive dates
    time_intervals = []
    interval_centers = []

    for i in range(len(dates_pd) - 1):
        # Time interval between consecutive acquisitions
        interval_days = (dates_pd[i+1] - dates_pd[i]).days
        time_intervals.append(interval_days)

        # Center point of the interval (where we'll place the bar)
        center_date = dates_pd[i] + (dates_pd[i+1] - dates_pd[i]) / 2
        interval_centers.append(center_date)

    # Calculate centerline lengths for normalization
    centerline_lengths = []
    for centerline in centerlines:
        if centerline is not None:
            if isinstance(centerline, MultiLineString):
                # Sum lengths of all line segments
                total_length = sum(line.length for line in centerline.geoms)
            elif isinstance(centerline, LineString):
                total_length = centerline.length
            else:
                total_length = 0
            centerline_lengths.append(total_length)
        else:
            centerline_lengths.append(0)

    # For each time interval, use the average centerline length
    interval_centerline_lengths = []
    for i in range(len(centerline_lengths) - 1):
        # Average length between consecutive time steps
        avg_length = (centerline_lengths[i] + centerline_lengths[i+1]) / 2
        interval_centerline_lengths.append(avg_length)

    # Calculate areas for each time step
    deposition_areas = []
    erosion_areas = []

    for i, (bar_polys, erosion_polys) in enumerate(zip(bars, erosions)):
        # Calculate total deposition area
        depo_area = 0
        if bar_polys is not None:
            if hasattr(bar_polys, 'geoms'):  # MultiPolygon
                for poly in bar_polys.geoms:
                    if hasattr(poly, 'area'):
                        depo_area += poly.area
            elif hasattr(bar_polys, 'area'):  # Single Polygon
                depo_area = bar_polys.area

        # Calculate total erosion area  
        ero_area = 0
        if erosion_polys is not None:
            if hasattr(erosion_polys, 'geoms'):  # MultiPolygon
                for poly in erosion_polys.geoms:
                    if hasattr(poly, 'area'):
                        ero_area += poly.area
            elif hasattr(erosion_polys, 'area'):  # Single Polygon
                ero_area = erosion_polys.area

        deposition_areas.append(depo_area)
        erosion_areas.append(ero_area)

    # Create the plot
    fig, ax = plt.subplots(figsize=figsize)

    # Plot bars with widths proportional to time intervals
    for i, (center_date, interval_days, depo_area, ero_area, centerline_length) in enumerate(
        zip(interval_centers, time_intervals, deposition_areas, erosion_areas, interval_centerline_lengths)):

        # Bar width proportional to time interval
        width = interval_days * 1.0

        # Normalize rates by centerline length (avoid division by zero)
        if centerline_length > 0:
            depo_rate_normalized = (depo_area/interval_days) / centerline_length
            ero_rate_normalized = (ero_area/interval_days) / centerline_length
        else:
            depo_rate_normalized = 0
            ero_rate_normalized = 0

        # Plot deposition (positive)
        ax.bar(center_date, depo_rate_normalized, width=width, 
               color='xkcd:blue', edgecolor='k', alpha=0.8, label='Deposition' if i == 0 else "")

        # Plot erosion (negative)
        ax.bar(center_date, -ero_rate_normalized, width=width, 
               color='xkcd:faded red', edgecolor='k', alpha=0.8, label='Erosion' if i == 0 else "")

    # Customize the plot
    ax.set_xlabel('Date', fontsize=12)
    ax.set_ylabel('Rate of Change (m/day)', fontsize=12)
    ax.set_title('Deposition and Erosion Over Time', fontsize=14)

    # Add horizontal line at zero
    ax.axhline(y=0, color='black', linestyle='-', linewidth=0.8)

    # # Add vertical lines at acquisition times
    # for date in dates_pd:
    #     ax.axvline(x=date, color='gray', linestyle='-', alpha=0.5, linewidth=1)

    # Format x-axis dates
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
    ax.xaxis.set_major_locator(mdates.YearLocator())

    # Add minor ticks for better resolution
    ax.xaxis.set_minor_locator(mdates.MonthLocator())

    # Rotate labels for better readability
    plt.xticks(rotation=45)

    # Add legend
    ax.legend()

    # Set axis limits with some padding
    ax.set_xlim(dates_pd[0] - pd.Timedelta(days=30), 
                dates_pd[-1] + pd.Timedelta(days=30))

    # Adjust layout
    plt.tight_layout()

    # Print summary statistics
    print(f"Number of time intervals: {len(bars)}")
    print(f"Date range: {dates_pd[0].strftime('%Y-%m-%d')} to {dates_pd[-1].strftime('%Y-%m-%d')}")
    print(f"Average time interval: {np.mean(time_intervals):.1f} days")
    print(f"Average centerline length: {np.mean(interval_centerline_lengths):.2f} m")
    print(f"Total deposition: {np.sum(deposition_areas):.2f} m²")
    print(f"Total erosion: {np.sum(erosion_areas):.2f} m²")
    print(f"Net change: {np.sum(deposition_areas) - np.sum(erosion_areas):.2f} m²")

    # Create summary DataFrame with normalized rates
    normalized_depo_rates = []
    normalized_ero_rates = []

    for i, (depo_area, ero_area, interval_days, centerline_length) in enumerate(
        zip(deposition_areas, erosion_areas, time_intervals, interval_centerline_lengths)):

        if centerline_length > 0:
            normalized_depo_rates.append((depo_area/interval_days) / centerline_length)
            normalized_ero_rates.append((ero_area/interval_days) / centerline_length)
        else:
            normalized_depo_rates.append(0)
            normalized_ero_rates.append(0)

    summary_df = pd.DataFrame({
        'start_date': dates_pd[:-1],
        'end_date': dates_pd[1:],
        'center_date': interval_centers,
        'interval_days': time_intervals,
        'centerline_length_m': interval_centerline_lengths,
        'deposition_m2': deposition_areas,
        'erosion_m2': erosion_areas,
        'deposition_rate_m_per_day': normalized_depo_rates,
        'erosion_rate_m_per_day': normalized_ero_rates,
        'net_change_m2': np.array(deposition_areas) - np.array(erosion_areas),
        'net_rate_m_per_day': np.array(normalized_depo_rates) - np.array(normalized_ero_rates)
    })

    return fig, ax, summary_df

map_graphs_over_time(D_primal_t1, G_rook_t1, D_primal_t2, G_rook_t2, rook_iou_threshold=0.5, primal_dist_threshold=100.0, primal_edge_sim_threshold=0.5, primal_edge_buffer=50.0)

Maps graph components (nodes and edges) between two time steps (t1 and t2). Mapping is based on the similarity of the locations and shapes of centerlines, centerline nodes, and banklines/islands. Args: D_primal_t1 (nx.MultiDiGraph): Directed centerline graph at time t1. G_rook_t1 (nx.Graph): Bankline rook graph at time t1. D_primal_t2 (nx.MultiDiGraph): Directed centerline graph at time t2. G_rook_t2 (nx.Graph): Bankline rook graph at time t2. rook_iou_threshold (float): IoU threshold for mapping G_rook nodes. primal_dist_threshold (float): Distance threshold for mapping D_primal nodes. primal_edge_sim_threshold (float): Similarity threshold for mapping D_primal edges. primal_edge_buffer (float): Buffer distance for D_primal edge similarity calculation. Returns: dict: A dictionary containing mappings for G_rook nodes, D_primal nodes, and D_primal edges.

Source code in rivabar/temporal_analysis.py

def map_graphs_over_time(D_primal_t1, G_rook_t1, D_primal_t2, G_rook_t2,
                         rook_iou_threshold=0.5, primal_dist_threshold=100.0,
                         primal_edge_sim_threshold=0.5, primal_edge_buffer=50.0):
    """
    Maps graph components (nodes and edges) between two time steps (t1 and t2).
    Mapping is based on the similarity of the locations and shapes of centerlines, 
    centerline nodes, and banklines/islands.
    Args:
        D_primal_t1 (nx.MultiDiGraph): Directed centerline graph at time t1.
        G_rook_t1 (nx.Graph): Bankline rook graph at time t1.
        D_primal_t2 (nx.MultiDiGraph): Directed centerline graph at time t2.
        G_rook_t2 (nx.Graph): Bankline rook graph at time t2.
        rook_iou_threshold (float): IoU threshold for mapping G_rook nodes.
        primal_dist_threshold (float): Distance threshold for mapping D_primal nodes.
        primal_edge_sim_threshold (float): Similarity threshold for mapping D_primal edges.
        primal_edge_buffer (float): Buffer distance for D_primal edge similarity calculation.
    Returns:
        dict: A dictionary containing mappings for G_rook nodes, D_primal nodes, and D_primal edges.
    """
    print("Mapping G_rook nodes...")
    rook_node_mapping = _map_rook_nodes(G_rook_t1, G_rook_t2, iou_threshold=rook_iou_threshold)

    print("Mapping D_primal nodes...")
    primal_node_mapping = _map_primal_nodes(D_primal_t1, D_primal_t2, distance_threshold=primal_dist_threshold)

    print("Mapping D_primal edges...")
    primal_edge_mapping = _map_primal_edges(D_primal_t1, D_primal_t2, 
                                            similarity_threshold=primal_edge_sim_threshold, 
                                            buffer_dist=primal_edge_buffer)

    mappings = {
        'rook_nodes': rook_node_mapping,
        'primal_nodes': primal_node_mapping,
        'primal_edges': primal_edge_mapping,
    }

    print("Graph mapping complete.")
    return mappings

calculate_node_displacement_deviation(D_primal_t1, D_primal_t2, primal_node_mapping)

Calculates how much the node displacement vector deviates from the average orientation of the connected centerline edges, and the displacement distance.

For each mapped node, this function computes two vectors: 1. The displacement vector, from the node's position at t1 to its new position at t2. 2. The average orientation vector of all connected centerline edges at t1.

It then calculates the angle between these two vectors and the length of the displacement.

Args: D_primal_t1 (nx.MultiGraph): Centerline graph at time t1. D_primal_t2 (nx.MultiGraph): Centerline graph at time t2. primal_node_mapping (dict): A dictionary mapping node IDs from t1 to t2.

Returns: tuple: A tuple containing two dictionaries: - deviations (dict): Node IDs from t1 to deviation angles in degrees. - distances (dict): Node IDs from t1 to displacement distances.

Source code in rivabar/temporal_analysis.py

def calculate_node_displacement_deviation(D_primal_t1, D_primal_t2, primal_node_mapping):
    """
    Calculates how much the node displacement vector deviates from the average 
    orientation of the connected centerline edges, and the displacement distance.

    For each mapped node, this function computes two vectors:
    1. The displacement vector, from the node's position at t1 to its new position at t2.
    2. The average orientation vector of all connected centerline edges at t1.

    It then calculates the angle between these two vectors and the length of the displacement.

    Args:
        D_primal_t1 (nx.MultiGraph): Centerline graph at time t1.
        D_primal_t2 (nx.MultiGraph): Centerline graph at time t2.
        primal_node_mapping (dict): A dictionary mapping node IDs from t1 to t2.

    Returns:
        tuple: A tuple containing two dictionaries:
               - deviations (dict): Node IDs from t1 to deviation angles in degrees.
               - distances (dict): Node IDs from t1 to displacement distances.
    """
    deviations = {}
    distances = {}

    for n1, n2 in tqdm(primal_node_mapping.items(), desc="Calculating node deviations"):
        if not D_primal_t1.has_node(n1) or not D_primal_t2.has_node(n2):
            continue

        p1 = np.array(D_primal_t1.nodes[n1]['geometry'].coords[0])
        p2 = np.array(D_primal_t2.nodes[n2]['geometry'].coords[0])
        displacement_vector = p2 - p1
        displacement_norm = np.linalg.norm(displacement_vector)
        distances[n1] = displacement_norm

        if displacement_norm == 0:
            deviations[n1] = 0.0
            continue

        normalized_displacement = displacement_vector / displacement_norm

        orientation_vectors = []
        # The graph may be a MultiGraph, which is undirected. Use .edges()
        for u, v, _, _ in D_primal_t1.edges(n1, data=True, keys=True):
            other_node = v if u == n1 else u

            p_other = np.array(D_primal_t1.nodes[other_node]['geometry'].coords[0])
            vector = p_other - p1

            vector_norm = np.linalg.norm(vector)
            if vector_norm > 0:
                orientation_vectors.append(vector / vector_norm)

        if not orientation_vectors:
            deviations[n1] = None
            continue

        avg_orientation_vector = np.mean(orientation_vectors, axis=0)
        avg_orientation_norm = np.linalg.norm(avg_orientation_vector)

        if avg_orientation_norm == 0:
            deviations[n1] = None
            continue

        normalized_avg_orientation = avg_orientation_vector / avg_orientation_norm

        dot_product = np.clip(np.dot(normalized_displacement, normalized_avg_orientation), -1.0, 1.0)
        angle_rad = np.arccos(dot_product)

        deviations[n1] = np.degrees(angle_rad)

    return deviations, distances

filter_rivers_by_length(rivers, std_threshold=0.3, pixel_size=30.0)

Filter rivers by removing anomalously short channels.

Computes each river's main-channel centerline length and mean channel width, then drops rivers whose centerline length falls more than std_threshold standard deviations below the mean length across all rivers. Useful for cleaning batch-processed scenes where the extraction only captured part of the reach.

Parameters:

Name	Type	Description	Default
rivers	list of River	Processed River objects.	required
std_threshold	float	Number of standard deviations below the mean length to use as the cutoff (default 0.3; e.g. 1.0 is more permissive).	0.3
pixel_size	float	Pixel size in meters used to convert widths for rivers whose raster data has been cleared (default 30). When a river still has its dataset (or a saved transform), the pixel size is taken from there.	30.0

Returns:

Name	Type	Description
filtered_rivers	list of River	Rivers whose centerline length exceeds the threshold.
filtered_lengths	list of float	Centerline lengths of the filtered rivers (m).
valid_indices	list of int	Indices of the filtered rivers in the input list.
ch_lengths	list of float	Centerline lengths of all input rivers (0 where unavailable).
ch_widths	list of float	Mean channel widths of all input rivers (m; 0 where unavailable).
threshold	float	The length threshold that was applied (m).

Source code in rivabar/temporal_analysis.py

def filter_rivers_by_length(rivers, std_threshold=0.3, pixel_size=30.0):
    """
    Filter rivers by removing anomalously short channels.

    Computes each river's main-channel centerline length and mean channel
    width, then drops rivers whose centerline length falls more than
    ``std_threshold`` standard deviations below the mean length across all
    rivers. Useful for cleaning batch-processed scenes where the extraction
    only captured part of the reach.

    Parameters
    ----------
    rivers : list of River
        Processed River objects.
    std_threshold : float, optional
        Number of standard deviations below the mean length to use as the
        cutoff (default 0.3; e.g. 1.0 is more permissive).
    pixel_size : float, optional
        Pixel size in meters used to convert widths for rivers whose raster
        data has been cleared (default 30). When a river still has its
        dataset (or a saved transform), the pixel size is taken from there.

    Returns
    -------
    filtered_rivers : list of River
        Rivers whose centerline length exceeds the threshold.
    filtered_lengths : list of float
        Centerline lengths of the filtered rivers (m).
    valid_indices : list of int
        Indices of the filtered rivers in the input list.
    ch_lengths : list of float
        Centerline lengths of all input rivers (0 where unavailable).
    ch_widths : list of float
        Mean channel widths of all input rivers (m; 0 where unavailable).
    threshold : float
        The length threshold that was applied (m).
    """
    ch_lengths = []
    ch_widths = []
    for river in rivers:
        try:
            cl = river.main_channel_centerline
        except Exception:
            cl = None  # unprocessed/failed rivers count as zero-length
        if cl is not None:
            ch_lengths.append(cl.length)
            ps = pixel_size
            if river._dataset is not None and getattr(river._dataset, 'transform', None) is not None:
                ps = river._dataset.transform[0]
            try:
                s, widths = river.get_channel_widths(pixel_size=ps)
                ch_widths.append(float(np.nanmean(widths)))
            except Exception as e:
                print(f"Could not compute widths for {getattr(river, 'scene_id', None) or river.fname}: {e}")
                ch_widths.append(0)
        else:
            print(f"No main channel centerline for {getattr(river, 'scene_id', None) or river.fname}")
            ch_lengths.append(0)
            ch_widths.append(0)

    # Filter out anomalously low channel lengths
    mean_length = np.mean(ch_lengths)
    std_length = np.std(ch_lengths)
    threshold = mean_length - std_threshold * std_length

    valid_indices = [i for i, length in enumerate(ch_lengths) if length > threshold]
    filtered_rivers = [rivers[i] for i in valid_indices]
    filtered_lengths = [ch_lengths[i] for i in valid_indices]

    print(f"Original number of rivers: {len(rivers)}")
    print(f"Filtered number of rivers: {len(filtered_rivers)}")
    print(f"Removed {len(rivers) - len(filtered_rivers)} anomalously short channels")

    return filtered_rivers, filtered_lengths, valid_indices, ch_lengths, ch_widths, threshold

find_common_confluences(rivers, min_scene_fraction=0.5, width_scale_factor=3.0, min_branch_length=None)

Identify persistent tributary confluence locations across multiple river scenes.

Collects tributary confluences from all rivers, clusters them spatially (using a distance threshold scaled to the mean channel width), and retains only clusters that appear in at least min_scene_fraction of the scenes.

Parameters:

Name	Type	Description	Default
rivers	list of River	Processed River objects (typically from the same Landsat path/row). Each must have been processed with map_river_banks().	required
min_scene_fraction	float	Minimum fraction of scenes in which a confluence must appear to be considered persistent (default 0.5).	0.5
width_scale_factor	float	The clustering distance threshold is set to width_scale_factor * mean_channel_width (default 3.0).	3.0
min_branch_length	float	If set, only consider tributary branches with branch_length_pixels >= min_branch_length.	None

Returns:

Name	Type	Description
common_confluences	list of dict	Each dict contains: - 'utm_coords': (x, y) centroid of the cluster in UTM - 'n_scenes': number of scenes the confluence was found in - 'scene_fraction': fraction of scenes - 'member_coords': list of (x, y) coordinates from individual scenes

Source code in rivabar/temporal_analysis.py

def find_common_confluences(rivers, min_scene_fraction=0.5, width_scale_factor=3.0,
                            min_branch_length=None):
    """
    Identify persistent tributary confluence locations across multiple river scenes.

    Collects tributary confluences from all rivers, clusters them spatially
    (using a distance threshold scaled to the mean channel width), and retains
    only clusters that appear in at least ``min_scene_fraction`` of the scenes.

    Parameters
    ----------
    rivers : list of River
        Processed River objects (typically from the same Landsat path/row).
        Each must have been processed with ``map_river_banks()``.
    min_scene_fraction : float, optional
        Minimum fraction of scenes in which a confluence must appear to be
        considered persistent (default 0.5).
    width_scale_factor : float, optional
        The clustering distance threshold is set to
        ``width_scale_factor * mean_channel_width`` (default 3.0).
    min_branch_length : float, optional
        If set, only consider tributary branches with
        ``branch_length_pixels >= min_branch_length``.

    Returns
    -------
    common_confluences : list of dict
        Each dict contains:
        - 'utm_coords': (x, y) centroid of the cluster in UTM
        - 'n_scenes': number of scenes the confluence was found in
        - 'scene_fraction': fraction of scenes
        - 'member_coords': list of (x, y) coordinates from individual scenes
    """
    # Collect all confluence points with scene indices
    all_points = []  # list of (x, y, river_index)
    valid_rivers = []
    for i, river in enumerate(rivers):
        if not river._is_processed or not river._processing_successful:
            continue
        valid_rivers.append(i)
        for trib in river.tributary_confluences:
            if 'confluence_utm_coords' not in trib:
                continue
            if min_branch_length is not None and trib['branch_length_pixels'] < min_branch_length:
                continue
            x, y = trib['confluence_utm_coords']
            all_points.append((x, y, i))

    if not all_points:
        print("No tributary confluences found across the provided rivers.")
        return []

    n_valid = len(valid_rivers)
    min_scenes = max(1, int(np.ceil(min_scene_fraction * n_valid)))

    # Estimate mean channel width across all valid rivers (in meters)
    mean_width, measured = _estimate_mean_channel_width([rivers[i] for i in valid_rivers])
    if not measured:
        print(f"Could not estimate channel width; using default clustering distance of {mean_width * width_scale_factor:.0f} m")

    cluster_dist = width_scale_factor * mean_width
    print(f"Mean channel width: {mean_width:.0f} m, clustering distance: {cluster_dist:.0f} m")

    # Agglomerative clustering of confluence points
    coords = np.array([(p[0], p[1]) for p in all_points])
    scene_ids = np.array([p[2] for p in all_points])

    # Use scipy hierarchical clustering with the distance threshold
    from scipy.cluster.hierarchy import fcluster, linkage
    if len(coords) == 1:
        labels = np.array([1])
    else:
        Z = linkage(coords, method='average', metric='euclidean')
        labels = fcluster(Z, t=cluster_dist, criterion='distance')

    # Analyze clusters
    common_confluences = []
    for label in np.unique(labels):
        mask = labels == label
        cluster_coords = coords[mask]
        cluster_scenes = scene_ids[mask]

        # Count unique scenes in this cluster
        unique_scenes = set(cluster_scenes)
        n_scenes = len(unique_scenes)
        scene_fraction = n_scenes / n_valid

        if n_scenes >= min_scenes:
            centroid = cluster_coords.mean(axis=0)
            common_confluences.append({
                'utm_coords': (float(centroid[0]), float(centroid[1])),
                'n_scenes': n_scenes,
                'scene_fraction': scene_fraction,
                'member_coords': [(float(c[0]), float(c[1])) for c in cluster_coords],
            })

    # Sort by along-channel position (approximate: use distance from first river's start point)
    if common_confluences and valid_rivers:
        ref_river = rivers[valid_rivers[0]]
        sx, sy = ref_river.start_x, ref_river.start_y
        common_confluences.sort(
            key=lambda c: (c['utm_coords'][0] - sx)**2 + (c['utm_coords'][1] - sy)**2
        )

    print(f"Found {len(common_confluences)} persistent confluences "
          f"(out of {len(np.unique(labels))} total clusters, "
          f"threshold: present in >= {min_scenes}/{n_valid} scenes)")

    return common_confluences

match_river_segments(rivers, common_confluences, max_snapping_distance=None, width_scale_factor=2.0, min_rivers_per_segment=4)

Split all rivers at common confluences and group corresponding segments.

For each river, the main path is split at the common confluence locations. A segment is kept only if its bounding confluence points snap within max_snapping_distance of the river's centerline. Segments are then grouped by index (i.e., segment 0 from all qualifying rivers form one group, segment 1 from all qualifying rivers form another, etc.). Groups with fewer than min_rivers_per_segment rivers are dropped.

Parameters:

Name	Type	Description	Default
rivers	list of River	Processed River objects.	required
common_confluences	list of dict	Output of :func:find_common_confluences.	required
max_snapping_distance	float	Maximum allowed distance (in metres) between a confluence point and the nearest point on a river's centerline. Segments bounded by a confluence that exceeds this threshold are excluded for that river. If None (default), the threshold is set to width_scale_factor * mean_channel_width.	None
width_scale_factor	float	Multiplier for mean channel width to set the default snapping threshold (default 2.0). Ignored if max_snapping_distance is given.	2.0
min_rivers_per_segment	int	Minimum number of rivers that must contribute a segment for it to be included in the output (default 4).	4

Returns:

Name	Type	Description
segment_groups	list of dict	One entry per valid segment index. Each dict contains: 'segment_index': int — position along the channel (0 = most upstream). 'rivers': list of River — rivers that have a valid segment here. 'paths': list of list — the sub-path edge lists, one per river (same order as 'rivers'). 'n_rivers': int — number of rivers in this group. 'upstream_confluence': (x, y) or None for the first segment. 'downstream_confluence': (x, y) or None for the last segment.
rejected	dict	Keyed by (river_index, segment_index) with the reason string.

Source code in rivabar/temporal_analysis.py

def match_river_segments(rivers, common_confluences, max_snapping_distance=None,
                         width_scale_factor=2.0, min_rivers_per_segment=4):
    """
    Split all rivers at common confluences and group corresponding segments.

    For each river, the main path is split at the common confluence locations.
    A segment is kept only if its bounding confluence points snap within
    ``max_snapping_distance`` of the river's centerline. Segments are then
    grouped by index (i.e., segment 0 from all qualifying rivers form one
    group, segment 1 from all qualifying rivers form another, etc.). Groups
    with fewer than ``min_rivers_per_segment`` rivers are dropped.

    Parameters
    ----------
    rivers : list of River
        Processed River objects.
    common_confluences : list of dict
        Output of :func:`find_common_confluences`.
    max_snapping_distance : float, optional
        Maximum allowed distance (in metres) between a confluence point and
        the nearest point on a river's centerline. Segments bounded by a
        confluence that exceeds this threshold are excluded for that river.
        If *None* (default), the threshold is set to
        ``width_scale_factor * mean_channel_width``.
    width_scale_factor : float, optional
        Multiplier for mean channel width to set the default snapping
        threshold (default 2.0). Ignored if *max_snapping_distance* is given.
    min_rivers_per_segment : int, optional
        Minimum number of rivers that must contribute a segment for it to be
        included in the output (default 4).

    Returns
    -------
    segment_groups : list of dict
        One entry per valid segment index. Each dict contains:

        - ``'segment_index'``: int — position along the channel (0 = most
          upstream).
        - ``'rivers'``: list of River — rivers that have a valid segment here.
        - ``'paths'``: list of list — the sub-path edge lists, one per river
          (same order as ``'rivers'``).
        - ``'n_rivers'``: int — number of rivers in this group.
        - ``'upstream_confluence'``: (x, y) or *None* for the first segment.
        - ``'downstream_confluence'``: (x, y) or *None* for the last segment.
    rejected : dict
        Keyed by ``(river_index, segment_index)`` with the reason string.
    """
    split_points = [c['utm_coords'] for c in common_confluences]
    n_confluences = len(split_points)
    n_segments = n_confluences + 1  # segments between / outside confluences

    # Estimate mean channel width for default threshold
    if max_snapping_distance is None:
        mean_width, _ = _estimate_mean_channel_width(rivers)
        max_snapping_distance = width_scale_factor * mean_width
        print(f"Max snapping distance: {max_snapping_distance:.0f} m "
              f"({width_scale_factor}x mean width of {mean_width:.0f} m)")

    # Split each river at all confluence points, then merge segments
    # across unreachable confluences so that the split boundaries align
    # with the confluence markers in the visualization.
    #
    # For each river we:
    #   1. Split at ALL confluence points (to obtain snapping distances).
    #   2. Identify which confluences are reachable (snap within threshold).
    #   3. Collapse segments across unreachable confluences by merging
    #      their edge lists.
    #   4. Assign each merged segment to the standard group keyed by its
    #      most-downstream reachable confluence boundary.
    #
    # Group key: (up_conf_idx or None, down_conf_idx or None).
    from collections import defaultdict
    group_collector = defaultdict(lambda: {'rivers': [], 'paths': []})
    rejected = {}

    for ri, river in enumerate(rivers):
        if not river._is_processed or not river._processing_successful:
            rejected[(ri, 'all')] = 'river not processed'
            continue
        if river.main_path is None:
            rejected[(ri, 'all')] = 'no main path'
            continue

        try:
            segments, split_info = river.split_main_path_at_points(split_points)
        except Exception as e:
            rejected[(ri, 'all')] = f'split failed: {e}'
            continue

        # split_info and segments are in along-channel order, which can
        # differ from the input order of common_confluences on sinuous
        # rivers; conf_input_idx maps each along-channel boundary back to
        # its index in common_confluences so that group keys are consistent
        # across rivers.
        snap_dists = [info['snapping_distance'] for info in split_info]
        conf_reachable = [d <= max_snapping_distance for d in snap_dists]
        conf_input_idx = [info.get('input_index', j)
                          for j, info in enumerate(split_info)]

        # Defensive: segments should always have n_confluences + 1 entries
        while len(segments) < n_segments:
            segments.append([])

        # If no confluences are reachable, reject the entire river
        if not any(conf_reachable):
            rejected[(ri, 'all')] = (
                f'no reachable confluences '
                f'(snap distances: {[f"{d:.0f}" for d in snap_dists]}, '
                f'threshold: {max_snapping_distance:.0f} m)')
            continue

        # Merge segments across unreachable confluences.
        # Walk through segments 0..n_segments-1.  At each reachable
        # confluence boundary, finalise the current merged segment;
        # at unreachable boundaries, keep accumulating edges.
        merged = []          # list of (edge_list, downstream_conf_index_or_None)
        current_path = list(segments[0]) if segments[0] else []

        for ci in range(n_confluences):
            next_seg = segments[ci + 1] if (ci + 1) < len(segments) else []
            if conf_reachable[ci]:
                # Valid boundary — save current merged segment, keyed by
                # the confluence's index in common_confluences
                merged.append((current_path, conf_input_idx[ci]))
                current_path = list(next_seg)
            else:
                # Unreachable — merge next segment into current
                current_path.extend(next_seg)

        # Final segment (extends to river end)
        merged.append((current_path, None))

        # Assign each merged segment to a group keyed by its actual
        # upstream and downstream confluence boundaries.  When a
        # confluence is unreachable the segment spans across it, so
        # the key reflects the true (possibly wider) extent — e.g.
        # (None, 1) instead of the standard (0, 1).
        upstream_boundary = None  # None = river start
        for path, down_ci in merged:
            if not path:
                # Advance upstream boundary even for empty segments
                upstream_boundary = down_ci
                continue
            group_key = (upstream_boundary, down_ci)
            group_collector[group_key]['rivers'].append(river)
            group_collector[group_key]['paths'].append(path)
            upstream_boundary = down_ci

    # Apply minimum count and build output
    segment_groups = []
    for group_key in sorted(group_collector.keys(),
                            key=lambda k: (k[0] if k[0] is not None else -1)):
        g = group_collector[group_key]
        up_idx, down_idx = group_key

        if len(g['rivers']) < min_rivers_per_segment:
            for river in g['rivers']:
                ri = rivers.index(river)
                rejected[(ri, group_key)] = (
                    f'segment group too small '
                    f'({len(g["rivers"])}/{min_rivers_per_segment})')
            continue

        upstream_conf = split_points[up_idx] if up_idx is not None else None
        downstream_conf = split_points[down_idx] if down_idx is not None else None

        segment_groups.append({
            'segment_index': group_key,
            'rivers': g['rivers'],
            'paths': g['paths'],
            'n_rivers': len(g['rivers']),
            'upstream_confluence': upstream_conf,
            'downstream_confluence': downstream_conf,
        })

    # Summary
    total_valid = sum(g['n_rivers'] for g in segment_groups)
    print(f"Matched {len(segment_groups)} segment groups, "
          f"{total_valid} river-segment pairs "
          f"(rejected {len(rejected)})")
    for g in segment_groups:
        up = g['segment_index'][0]
        down = g['segment_index'][1]
        up_label = f'confluence {up}' if up is not None else 'start'
        down_label = f'confluence {down}' if down is not None else 'end'
        print(f"  {up_label} -> {down_label}: {g['n_rivers']} rivers")

    return segment_groups, rejected

match_rivers_to_images(rivers, image_directory, tolerance_days=1)

Match processed rivers to georeferenced images by acquisition date.

Scans image_directory for GeoTIFFs (e.g., false-color images downloaded with River.batch_process_landsat_scenes(download_false_color=True)), parses acquisition dates from the filenames, and pairs each image with the river whose acquisition date is closest (within tolerance_days).

Parameters:

Name	Type	Description	Default
rivers	list of River	Processed River objects with an acquisition_date attribute ('YYYY-MM-DD').	required
image_directory	str or Path	Directory containing .tif images with dates in their filenames. Supported filename patterns include Landsat scene/product IDs (e.g. LC08_231064_20200515, LC08_L2SP_232060_20140219_...) and generic YYYYMMDD or YYYY-MM-DD dates.	required
tolerance_days	int	Maximum allowed difference between image and river acquisition dates (default 1).	1

Returns:

Name	Type	Description
matched_rivers	list of River	Rivers with a matching image, in image-date order.
image_files	list of pathlib.Path	The matched image files (same order and length as matched_rivers).
dates	list of datetime.datetime	Image acquisition dates (same order and length as matched_rivers).

Source code in rivabar/temporal_analysis.py

def match_rivers_to_images(rivers, image_directory, tolerance_days=1):
    """
    Match processed rivers to georeferenced images by acquisition date.

    Scans *image_directory* for GeoTIFFs (e.g., false-color images downloaded
    with ``River.batch_process_landsat_scenes(download_false_color=True)``),
    parses acquisition dates from the filenames, and pairs each image with
    the river whose acquisition date is closest (within *tolerance_days*).

    Parameters
    ----------
    rivers : list of River
        Processed River objects with an ``acquisition_date`` attribute
        ('YYYY-MM-DD').
    image_directory : str or pathlib.Path
        Directory containing .tif images with dates in their filenames.
        Supported filename patterns include Landsat scene/product IDs
        (e.g. ``LC08_231064_20200515``, ``LC08_L2SP_232060_20140219_...``)
        and generic ``YYYYMMDD`` or ``YYYY-MM-DD`` dates.
    tolerance_days : int, optional
        Maximum allowed difference between image and river acquisition
        dates (default 1).

    Returns
    -------
    matched_rivers : list of River
        Rivers with a matching image, in image-date order.
    image_files : list of pathlib.Path
        The matched image files (same order and length as matched_rivers).
    dates : list of datetime.datetime
        Image acquisition dates (same order and length as matched_rivers).
    """
    from pathlib import Path
    import re

    image_directory = Path(image_directory)
    tif_files = sorted(set(list(image_directory.glob('*.tif')) +
                           list(image_directory.glob('*.TIF'))))
    if not tif_files:
        raise ValueError(f"No TIF files found in {image_directory}")

    date_patterns = [
        # Landsat scene IDs: LC08_231064_20200515, false_color_LC08_231064_20200515
        r'(?:false_color_)?L[CETM]\d{2}_\d{6}_(\d{8})',
        # Landsat product IDs: LC08_L2SP_232060_20140219_20200911_02_T1
        r'L[CETM]\d{2}_L\d\w{2}_\d{6}_(\d{8})',
        # Generic date patterns: YYYYMMDD
        r'(\d{8})',
        # Date with separators: YYYY-MM-DD, YYYY_MM_DD
        r'(\d{4}[-_]\d{2}[-_]\d{2})',
    ]

    dated_files = []
    for file_path in tif_files:
        for pattern in date_patterns:
            match = re.search(pattern, file_path.name)
            if match:
                date_str = match.group(1).replace('-', '').replace('_', '')
                try:
                    date_obj = datetime.strptime(date_str, '%Y%m%d')
                except ValueError:
                    continue
                dated_files.append((file_path, date_obj))
                break
        else:
            print(f"Could not parse a date from {file_path.name}; skipping")
    dated_files.sort(key=lambda x: x[1])

    river_dates = [(river, datetime.strptime(river.acquisition_date, '%Y-%m-%d'))
                   for river in rivers]

    # For each image, find the closest river within tolerance; each river is
    # used at most once (closest image wins)
    candidates = []  # (day_difference, image_index, river)
    for i, (file_path, image_date) in enumerate(dated_files):
        for river, river_date in river_dates:
            diff = abs((image_date - river_date).days)
            if diff <= tolerance_days:
                candidates.append((diff, i, river))
    candidates.sort(key=lambda c: c[0])
    used_images, used_rivers = set(), set()
    pairs = {}
    for diff, i, river in candidates:
        if i in used_images or id(river) in used_rivers:
            continue
        used_images.add(i)
        used_rivers.add(id(river))
        pairs[i] = river

    matched_indices = sorted(pairs.keys())
    matched_rivers = [pairs[i] for i in matched_indices]
    image_files = [dated_files[i][0] for i in matched_indices]
    dates = [dated_files[i][1] for i in matched_indices]

    print(f"Matched {len(matched_rivers)} of {len(rivers)} rivers to "
          f"{len(dated_files)} images (tolerance: {tolerance_days} day(s))")
    return matched_rivers, image_files, dates