Calculate_CFI.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a48ab391",
   "metadata": {},
   "source": [
    "# Calculate channel form index\n",
    "\n",
    "The following code is directed to a given local path containing 2-D water mask rasters. The code takes the water mask, and start by creating a \"skeleton\" of the mask. It then dilates the tips of the skeleton to improve connection of the channel network, reskeletonizes, and reduces the skeleton to only the identifiable river channels. From the final skeletion, channel links and nodes are created. The links are filtered according to criteria, and a shortest path line or \"main channel\" is extracted, along with a simplified main channel which acts as a valley center line, enabling sinuosity calculations. Finally, cross sections of the river are created and channel count index is calculated across the cross-sections. With sinuosity and channel-count index, the chanel form index can be calculated. These river metrics are provided and exported to a .csv. The processed skeleton, nodes, channel links, main channel, valley center-line, and channel-belt cross-sections are output as shapefiles, and the network is plotted for the user's convenience.\n",
    "\n",
    "Author: James (Huck) Rees; PhD Student, UCSB Geography\n",
    "\n",
    "Date: June 7th, 2024"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b19a039",
   "metadata": {},
   "source": [
    "## Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "facbeca7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import re\n",
    "import logging\n",
    "from glob import glob\n",
    "from itertools import combinations\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "\n",
    "import rasterio\n",
    "from rasterio.plot import show\n",
    "from rasterio.transform import xy\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from skimage.morphology import skeletonize, label\n",
    "from skimage.measure import regionprops\n",
    "\n",
    "from shapefile import Reader, Writer\n",
    "from shapely.geometry import LineString, Point, MultiLineString, MultiPoint\n",
    "from shapely.ops import split, linemerge, snap, nearest_points\n",
    "from shapely.strtree import STRtree\n",
    "\n",
    "from rtree import index as rtree_index\n",
    "\n",
    "import networkx as nx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d84ad93d",
   "metadata": {},
   "source": [
    "## Initialize functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "62ef65ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to load raster data\n",
    "def load_raster(file_path):\n",
    "    \"\"\"\n",
    "    Loads a raster file and returns the data of the first band along with its metadata.\n",
    "\n",
    "    Parameters:\n",
    "    file_path (str): The path to the raster file.\n",
    "\n",
    "    Returns:\n",
    "    tuple: A tuple containing:\n",
    "        - data (numpy.ndarray): The data of the first band of the raster.\n",
    "        - metadata (dict): The metadata of the raster file.\n",
    "    \"\"\"\n",
    "    with rasterio.open(file_path) as dataset:\n",
    "        data = dataset.read(1)  # Read the first band\n",
    "        metadata = dataset.meta\n",
    "    return data, metadata\n",
    "\n",
    "# Function to save a raster file\n",
    "def save_raster(output_path, data, metadata):\n",
    "    \"\"\"\n",
    "    Saves a raster file with the given data and metadata.\n",
    "\n",
    "    Parameters:\n",
    "    output_path (str): The path to save the output raster file.\n",
    "    data (numpy.ndarray): The data to be written to the raster file.\n",
    "    metadata (dict): The metadata of the raster file, including CRS and transform information.\n",
    "\n",
    "    Returns:\n",
    "    None\n",
    "    \"\"\"\n",
    "    with rasterio.open(\n",
    "        output_path, \n",
    "        'w', \n",
    "        driver='GTiff', \n",
    "        height=data.shape[0], \n",
    "        width=data.shape[1], \n",
    "        count=1, \n",
    "        dtype='uint8', \n",
    "        crs=metadata['crs'], \n",
    "        transform=metadata['transform']\n",
    "    ) as dst:\n",
    "        dst.write(data.astype('uint8'), 1)\n",
    "\n",
    "# Function to perform conditional dilation\n",
    "def conditional_dilation(image, radius=5):\n",
    "    \"\"\"\n",
    "    Performs a conditional dilation on a binary image. Pixels with a value of 1 that have \n",
    "    two or fewer neighbors with the same value will cause a dilation within a given radius.\n",
    "\n",
    "    Parameters:\n",
    "    image (numpy.ndarray): The input binary image (2D array) to be processed.\n",
    "    radius (int, optional): The radius for the dilation operation. Default is 5.\n",
    "\n",
    "    Returns:\n",
    "    numpy.ndarray: The dilated image.\n",
    "    \"\"\"\n",
    "    dilated_image = np.copy(image)\n",
    "    for row in range(1, image.shape[0] - 1):\n",
    "        for col in range(1, image.shape[1] - 1):\n",
    "            if image[row, col] == 1:\n",
    "                neighbors = image[row-1:row+2, col-1:col+2]\n",
    "                if np.sum(neighbors) <= 2:  # Include the pixel itself in the count\n",
    "                    dilated_image[max(0, row-radius):min(row+radius+1, image.shape[0]), \n",
    "                                  max(0, col-radius):min(col+radius+1, image.shape[1])] = 1\n",
    "    return dilated_image\n",
    "\n",
    "# Function to keep only the largest connected component\n",
    "def keep_largest_component(image):\n",
    "    \"\"\"\n",
    "    Identifies and retains the largest connected component in a binary image. All other components are removed.\n",
    "\n",
    "    Parameters:\n",
    "    image (numpy.ndarray): The input binary image (2D array).\n",
    "\n",
    "    Returns:\n",
    "    numpy.ndarray: A binary image containing only the largest connected component.\n",
    "    \"\"\"\n",
    "    labeled_image, num_features = label(image, connectivity=2, return_num=True)\n",
    "    if num_features == 0:\n",
    "        return image\n",
    "    regions = regionprops(labeled_image)\n",
    "    largest_region = max(regions, key=lambda r: r.area)\n",
    "    largest_component = (labeled_image == largest_region.label)\n",
    "    return largest_component\n",
    "\n",
    "# Function to create node shapefile and return node points\n",
    "def create_nodes(image, metadata):\n",
    "    \"\"\"\n",
    "    Identifies and creates nodes in a binary image. Nodes are classified as 'endpoint' or 'junction' based on their connectivity.\n",
    "\n",
    "    Parameters:\n",
    "    image (numpy.ndarray): The input binary image (2D array) to be processed.\n",
    "    metadata (dict): The metadata of the raster file, including transform information.\n",
    "\n",
    "    Returns:\n",
    "    geopandas.GeoDataFrame: A GeoDataFrame containing the nodes with their types and geometries.\n",
    "    \"\"\"\n",
    "    node_points = []\n",
    "    transform = metadata['transform']\n",
    "    node_id = 1\n",
    "\n",
    "    adjacent_positions = [\n",
    "        (-1, -1), (-1, 0), (-1, 1),\n",
    "        (0, -1),          (0, 1),\n",
    "        (1, -1),  (1, 0),  (1, 1)\n",
    "    ]\n",
    "\n",
    "    # Iterate over each pixel in the image, excluding the borders\n",
    "    for row in range(1, image.shape[0] - 1):\n",
    "        for col in range(1, image.shape[1] - 1):\n",
    "            if image[row, col] == 1:\n",
    "                # Find neighbors that are part of the segment\n",
    "                neighbors = [\n",
    "                    (row + dr, col + dc)\n",
    "                    for dr, dc in adjacent_positions\n",
    "                    if image[row + dr, col + dc] == 1\n",
    "                ]\n",
    "                count = len(neighbors)\n",
    "                if count == 1:  # Endpoints\n",
    "                    x, y = xy(transform, row, col)\n",
    "                    node_points.append((node_id, 'endpoint', Point(x, y)))\n",
    "                    node_id += 1\n",
    "                elif count == 3:  # Potential junctions\n",
    "                    is_not_on_one_side = not (\n",
    "                        (image[row - 1, col - 1] == 1 and image[row - 1, col] == 1 and image[row - 1, col + 1] == 1) or  # above\n",
    "                        (image[row + 1, col - 1] == 1 and image[row + 1, col] == 1 and image[row + 1, col + 1] == 1) or  # below\n",
    "                        (image[row - 1, col - 1] == 1 and image[row, col - 1] == 1 and image[row + 1, col - 1] == 1) or  # left\n",
    "                        (image[row - 1, col + 1] == 1 and image[row, col + 1] == 1 and image[row + 1, col + 1] == 1)     # right\n",
    "                    )\n",
    "\n",
    "                    if is_not_on_one_side:\n",
    "                        is_junction = True\n",
    "                        for i in range(len(neighbors)):\n",
    "                            for j in range(i + 1, len(neighbors)):\n",
    "                                if (abs(neighbors[i][0] - neighbors[j][0]), abs(neighbors[i][1] - neighbors[j][1])) in [(0, 1), (1, 0)]:\n",
    "                                    is_junction = False\n",
    "                                    break\n",
    "                            if not is_junction:\n",
    "                                break\n",
    "                        if is_junction:\n",
    "                            x, y = xy(transform, row, col)\n",
    "                            node_points.append((node_id, 'junction', Point(x, y)))\n",
    "                            node_id += 1\n",
    "                elif count >= 4:  # Nodes with 4 or more adjacent pixels\n",
    "                    direct_pairs = sum(\n",
    "                        1 for i in range(len(neighbors))\n",
    "                        for j in range(i + 1, len(neighbors))\n",
    "                        if (abs(neighbors[i][0] - neighbors[j][0]), abs(neighbors[i][1] - neighbors[j][1])) in [(0, 1), (1, 0)]\n",
    "                    )\n",
    "\n",
    "                    if direct_pairs < 2:\n",
    "                        x, y = xy(transform, row, col)\n",
    "                        node_points.append((node_id, 'junction', Point(x, y)))\n",
    "                        node_id += 1\n",
    "\n",
    "    # Create a GeoDataFrame\n",
    "    gdf = gpd.GeoDataFrame(node_points, columns=['node_id', 'type', 'geometry'])\n",
    "    gdf.set_crs(epsg=4326, inplace=True)\n",
    "    \n",
    "    return gdf\n",
    "\n",
    "# Function to create links shapefile\n",
    "def create_links(image, metadata):\n",
    "    \"\"\"\n",
    "    Identifies and creates links between adjacent pixels in a binary image. Links are represented as LineStrings.\n",
    "\n",
    "    Parameters:\n",
    "    image (numpy.ndarray): The input binary image (2D array) to be processed.\n",
    "    metadata (dict): The metadata of the raster file, including transform information.\n",
    "\n",
    "    Returns:\n",
    "    geopandas.GeoDataFrame: A GeoDataFrame containing the links as LineStrings.\n",
    "    \"\"\"\n",
    "    links = []\n",
    "    transform = metadata['transform']\n",
    "    link_id = 1\n",
    "\n",
    "    # Iterate over each pixel in the image\n",
    "    for row in range(image.shape[0]):\n",
    "        for col in range(image.shape[1]):\n",
    "            if image[row, col] == 1:  # Check if the pixel is part of a segment\n",
    "                # Identify neighboring pixels that are also part of the segment\n",
    "                neighbors = [\n",
    "                    (row + dr, col + dc) \n",
    "                    for dr, dc in [(-1, 0), (1, 0), (0, -1), (0, 1), (-1, -1), (-1, 1), (1, -1), (1, 1)] \n",
    "                    if 0 <= row + dr < image.shape[0] and 0 <= col + dc < image.shape[1] and image[row + dr, col + dc] == 1\n",
    "                ]\n",
    "                # Create LineString for each neighbor\n",
    "                for nr, nc in neighbors:\n",
    "                    x1, y1 = xy(transform, row, col)  # Convert pixel coordinates to spatial coordinates\n",
    "                    x2, y2 = xy(transform, nr, nc)\n",
    "                    line = LineString([(x1, y1), (x2, y2)])\n",
    "                    links.append((link_id, line))  # Append link to the list\n",
    "                    link_id += 1\n",
    "\n",
    "    # Remove duplicate links by sorting the coordinates of each LineString\n",
    "    unique_links = []\n",
    "    seen = set()\n",
    "    for link in links:\n",
    "        coords = tuple(sorted(link[1].coords))\n",
    "        if coords not in seen:\n",
    "            seen.add(coords)\n",
    "            unique_links.append(link)\n",
    "\n",
    "    # Create a GeoDataFrame from the unique links\n",
    "    gdf = gpd.GeoDataFrame(unique_links, columns=['id', 'geometry'])\n",
    "    \n",
    "    # Set the coordinate reference system (CRS)\n",
    "    gdf.set_crs(epsg=4326, inplace=True)\n",
    "\n",
    "    return gdf\n",
    "\n",
    "# Function to filter links\n",
    "def filter_links(gdf):\n",
    "    \"\"\"\n",
    "    Filters out diagonal links from a GeoDataFrame of line segments, retaining only those\n",
    "    that are not part of an intersection with horizontal and vertical links.\n",
    "\n",
    "    Parameters:\n",
    "    gdf (geopandas.GeoDataFrame): The input GeoDataFrame containing line segments.\n",
    "\n",
    "    Returns:\n",
    "    geopandas.GeoDataFrame: A filtered GeoDataFrame with certain diagonal links removed.\n",
    "    \"\"\"\n",
    "    # Function to categorize the line segments\n",
    "    def categorize_line(row):\n",
    "        if row['start_point'][1] == row['end_point'][1]:\n",
    "            return 'horizontal'\n",
    "        elif row['start_point'][0] == row['end_point'][0]:\n",
    "            return 'vertical'\n",
    "        else:\n",
    "            return 'diagonal'\n",
    "    \n",
    "    # Function to extract start and end coordinates of each line segment\n",
    "    def get_coordinates(geometry):\n",
    "        start_point = geometry.coords[0]\n",
    "        end_point = geometry.coords[1]\n",
    "        return start_point, end_point\n",
    "    \n",
    "    # Apply the function to get coordinates and categorize each segment\n",
    "    gdf[['start_point', 'end_point']] = gdf.apply(lambda row: get_coordinates(row.geometry), axis=1, result_type='expand')\n",
    "    gdf['category'] = gdf.apply(categorize_line, axis=1)\n",
    "    \n",
    "    # Initialize spatial indexes for horizontal and vertical links\n",
    "    idx_horizontal = rtree_index.Index()\n",
    "    idx_vertical = rtree_index.Index()\n",
    "    \n",
    "    for idx, row in gdf.iterrows():\n",
    "        if row['category'] == 'horizontal':\n",
    "            idx_horizontal.insert(idx, row['geometry'].bounds)\n",
    "        elif row['category'] == 'vertical':\n",
    "            idx_vertical.insert(idx, row['geometry'].bounds)\n",
    "    \n",
    "    diagonals_to_remove = set()\n",
    "    \n",
    "    # Loop through each diagonal link\n",
    "    for index, diag_row in gdf[gdf['category'] == 'diagonal'].iterrows():\n",
    "        diag_start = diag_row['start_point']\n",
    "        diag_end = diag_row['end_point']\n",
    "        diag_bounds = diag_row['geometry'].bounds\n",
    "        x_coords = {diag_start[0], diag_end[0]}\n",
    "        y_coords = {diag_start[1], diag_end[1]}\n",
    "        hor = ver = False\n",
    "        \n",
    "        # Find horizontal links intersecting with the diagonal link using spatial index\n",
    "        for hor_idx in idx_horizontal.intersection(diag_bounds):\n",
    "            hor_row = gdf.loc[hor_idx]\n",
    "            hor_start = hor_row['start_point']\n",
    "            hor_end = hor_row['end_point']\n",
    "            if (hor_start[1] in y_coords or hor_end[1] in y_coords) and (hor_start[0] in x_coords and hor_end[0] in x_coords):\n",
    "                hor = True\n",
    "                break\n",
    "        \n",
    "        # Find vertical links intersecting with the diagonal link using spatial index\n",
    "        for ver_idx in idx_vertical.intersection(diag_bounds):\n",
    "            ver_row = gdf.loc[ver_idx]\n",
    "            ver_start = ver_row['start_point']\n",
    "            ver_end = ver_row['end_point']\n",
    "            if (ver_start[0] in x_coords or ver_end[0] in x_coords) and (ver_start[1] in y_coords and ver_end[1] in y_coords):\n",
    "                ver = True\n",
    "                break\n",
    "        \n",
    "        # Mark the diagonal for removal if it satisfies both conditions\n",
    "        if hor and ver:\n",
    "            diagonals_to_remove.add(index)\n",
    "    \n",
    "    # Drop the identified diagonal links\n",
    "    filtered_links = gdf.drop(index=diagonals_to_remove)\n",
    "    \n",
    "    # Drop the unnecessary columns before returning\n",
    "    filtered_links = filtered_links.drop(columns=['start_point', 'end_point', 'category'])\n",
    "    \n",
    "    return filtered_links\n",
    "\n",
    "# Function to find furthest endpoints\n",
    "def find_furthest_endpoints(gdf_points):\n",
    "    \"\"\"\n",
    "    Finds the two furthest nodes in the geodataframe, which may be of type 'endpoint' or 'junction'.\n",
    "\n",
    "    Parameters:\n",
    "    gdf_points (geopandas.GeoDataFrame): The geodataframe of points (nodes).\n",
    "\n",
    "    Returns:\n",
    "    geopandas.GeoDataFrame: A GeoDataFrame containing the two furthest points.\n",
    "    \"\"\"\n",
    "    if len(gdf_points) < 2:\n",
    "        raise ValueError(\"Not enough points to find the furthest pair.\")\n",
    "    \n",
    "    max_distance = 0\n",
    "    furthest_pair = None\n",
    "    for (idx1, point1), (idx2, point2) in combinations(gdf_points.iterrows(), 2):\n",
    "        distance = point1.geometry.distance(point2.geometry)\n",
    "        if distance > max_distance:\n",
    "            max_distance = distance\n",
    "            furthest_pair = (point1, point2)\n",
    "    \n",
    "    furthest_geometries = [furthest_pair[0].geometry, furthest_pair[1].geometry]\n",
    "    start_end_pts = gpd.GeoDataFrame(geometry=furthest_geometries, crs=gdf_points.crs)\n",
    "    return start_end_pts\n",
    "\n",
    "# Function to prune network\n",
    "def prune_network(nodes, filtered_links, start_end_pts):\n",
    "    \"\"\"\n",
    "    Prunes a network by removing spur links that are not part of the main network.\n",
    "\n",
    "    Parameters:\n",
    "    nodes (geopandas.GeoDataFrame): The GeoDataFrame containing nodes with 'endpoint' types.\n",
    "    filtered_links (geopandas.GeoDataFrame): The GeoDataFrame containing filtered links (line segments).\n",
    "    start_end_pts (geopandas.GeoDataFrame): The GeoDataFrame containing start and end points to retain.\n",
    "\n",
    "    Returns:\n",
    "    geopandas.GeoDataFrame: A pruned GeoDataFrame with spur links removed.\n",
    "    \"\"\"\n",
    "    endpoints = nodes[nodes['type'] == 'endpoint']\n",
    "    small_ends = endpoints[~endpoints.geometry.apply(lambda x: any(start_end_pts.geometry.intersects(x)))]\n",
    "    spurs = gpd.GeoDataFrame(columns=['geometry'], geometry='geometry', crs=nodes.crs)\n",
    "    G = nx.Graph()\n",
    "    for idx, row in filtered_links.iterrows():\n",
    "        coords = list(row.geometry.coords)\n",
    "        G.add_edge(coords[0], coords[1], index=idx, geometry=row.geometry)\n",
    "\n",
    "    for idx, p1 in small_ends.iterrows():\n",
    "        pn1 = nodes[nodes.geometry == p1.geometry].iloc[0]\n",
    "        nodes_excluding_pn1 = nodes[nodes.geometry != pn1.geometry]\n",
    "        nearest_node_geom = nearest_points(pn1.geometry, nodes_excluding_pn1.unary_union)[1]\n",
    "        nearest_node = nodes[nodes.geometry == nearest_node_geom].iloc[0]\n",
    "        try:\n",
    "            path = nx.shortest_path(G, source=tuple(pn1.geometry.coords[0]), target=tuple(nearest_node.geometry.coords[0]))\n",
    "            path_geometries = [G.edges[path[i], path[i+1]]['geometry'] for i in range(len(path)-1)]\n",
    "            if len(path_geometries) <= 20:\n",
    "                spur_geometry = LineString([point for line in path_geometries for point in line.coords])\n",
    "                new_spur = gpd.GeoDataFrame([{'geometry': spur_geometry}], geometry='geometry', crs=nodes.crs)\n",
    "                spurs = pd.concat([spurs, new_spur], ignore_index=True)\n",
    "        except nx.NetworkXNoPath:\n",
    "            continue\n",
    "            \n",
    "    pruned_links = filtered_links.overlay(spurs, how='difference')\n",
    "    return pruned_links\n",
    "\n",
    "# Function to find shortest path\n",
    "def find_shortest_path(start_end_pts, filtered_links):\n",
    "    \"\"\"\n",
    "    Finds the shortest path between two points in a network of filtered links.\n",
    "\n",
    "    Parameters:\n",
    "    start_end_pts (geopandas.GeoDataFrame): The GeoDataFrame containing the start and end points.\n",
    "    filtered_links (geopandas.GeoDataFrame): The GeoDataFrame containing the network of links (line segments).\n",
    "\n",
    "    Returns:\n",
    "    geopandas.GeoDataFrame: A GeoDataFrame containing the shortest path as a LineString.\n",
    "    \"\"\"\n",
    "    G = nx.Graph()\n",
    "    for idx, row in filtered_links.iterrows():\n",
    "        line = row.geometry\n",
    "        for i in range(len(line.coords) - 1):\n",
    "            start = Point(line.coords[i])\n",
    "            end = Point(line.coords[i + 1])\n",
    "            distance = start.distance(end)\n",
    "            G.add_edge(tuple(start.coords[0]), tuple(end.coords[0]), weight=distance)\n",
    "    \n",
    "    start_point = tuple(start_end_pts.geometry.iloc[0].coords[0])\n",
    "    end_point = tuple(start_end_pts.geometry.iloc[1].coords[0])\n",
    "    shortest_path = nx.shortest_path(G, source=start_point, target=end_point, weight='weight')\n",
    "    shortest_path_coords = [Point(coord) for coord in shortest_path]\n",
    "    shortest_path_line = LineString(shortest_path_coords)\n",
    "    shortest_path_length = shortest_path_line.length\n",
    "    shortest_path_gdf = gpd.GeoDataFrame({'geometry': [shortest_path_line]}, crs=filtered_links.crs)\n",
    "    return shortest_path_gdf\n",
    "\n",
    "# Function to classify channels\n",
    "def classify_channels(filtered_links, shortest_path):\n",
    "    \"\"\"\n",
    "    Classifies links in a network as 'main_channel' or 'other' based on their relationship to the shortest path.\n",
    "\n",
    "    Parameters:\n",
    "    filtered_links (geopandas.GeoDataFrame): The GeoDataFrame containing the network of links (line segments).\n",
    "    shortest_path (geopandas.GeoDataFrame): The GeoDataFrame containing the shortest path as a LineString.\n",
    "\n",
    "    Returns:\n",
    "    geopandas.GeoDataFrame: A GeoDataFrame with an additional column 'chnl_cat' classifying each link.\n",
    "    \"\"\"\n",
    "    classified_links = filtered_links.copy()\n",
    "    classified_links['chnl_cat'] = 'other'\n",
    "    shortest_path_line = shortest_path.geometry.iloc[0]\n",
    "    \n",
    "    def is_main_channel(link):\n",
    "        return link.within(shortest_path_line)\n",
    "    \n",
    "    classified_links['chnl_cat'] = classified_links.apply(\n",
    "        lambda row: 'main_channel' if is_main_channel(row.geometry) else 'other',\n",
    "        axis=1\n",
    "    )\n",
    "    \n",
    "    return classified_links\n",
    "\n",
    "# Function to simplify shortest path\n",
    "def simplify_shortest_path(shortest_path, num_vertices=10):\n",
    "    \"\"\"\n",
    "    Simplifies the shortest path to a specified number of vertices.\n",
    "\n",
    "    Parameters:\n",
    "    shortest_path (geopandas.GeoDataFrame): The GeoDataFrame containing the shortest path as a LineString.\n",
    "    num_vertices (int, optional): The number of vertices for the simplified path. Default is 10.\n",
    "\n",
    "    Returns:\n",
    "    geopandas.GeoDataFrame: A GeoDataFrame containing the simplified shortest path as a LineString.\n",
    "    \"\"\"\n",
    "    original_line = shortest_path.geometry.iloc[0]\n",
    "    simplified_coords = [\n",
    "        original_line.interpolate(i / (num_vertices - 1), normalized=True).coords[0] \n",
    "        for i in range(num_vertices)\n",
    "    ]\n",
    "    simplified_line = LineString(simplified_coords)\n",
    "    simplified_path_gdf = gpd.GeoDataFrame({'geometry': [simplified_line]}, crs=shortest_path.crs)\n",
    "    return simplified_path_gdf\n",
    "\n",
    "# Function to create perpendicular lines\n",
    "def create_perpendicular_lines(simplified_path, num_lines=10, fraction_length=1/5):\n",
    "    \"\"\"\n",
    "    Creates perpendicular lines along the simplified path at equal intervals.\n",
    "\n",
    "    Parameters:\n",
    "    simplified_path (geopandas.GeoDataFrame): A GeoDataFrame containing the simplified path as a LineString.\n",
    "    num_lines (int): Number of perpendicular lines to create.\n",
    "    fraction_length (float): Fraction of the total path length for the length of each perpendicular line.\n",
    "\n",
    "    Returns:\n",
    "    geopandas.GeoDataFrame: A GeoDataFrame containing the perpendicular lines.\n",
    "    \"\"\"\n",
    "    # Extract the LineString from the GeoDataFrame\n",
    "    line = simplified_path.geometry.iloc[0]\n",
    "    line_length = line.length\n",
    "    \n",
    "    # Calculate spacing between perpendicular lines and half the length of each perpendicular line\n",
    "    spacing = line_length / num_lines\n",
    "    half_length = (line_length * fraction_length) / 2\n",
    "    \n",
    "    # Generate points at equal intervals along the line\n",
    "    points = [line.interpolate(i * spacing, normalized=False) for i in range(num_lines)]\n",
    "    \n",
    "    perpendicular_lines = []\n",
    "    \n",
    "    coords = list(line.coords)\n",
    "    \n",
    "    for idx, point in enumerate(points):\n",
    "        # Find the segment that the point falls on\n",
    "        segment = None\n",
    "        for i in range(len(coords) - 1):\n",
    "            segment_line = LineString([coords[i], coords[i+1]])\n",
    "            if segment_line.project(point) < segment_line.length:\n",
    "                segment = segment_line\n",
    "                break\n",
    "        \n",
    "        if segment is None:\n",
    "            print(f\"No segment found for point {idx}: {point}\")\n",
    "            continue\n",
    "        \n",
    "        # Calculate the perpendicular direction to the segment\n",
    "        dx = segment.coords[1][0] - segment.coords[0][0]\n",
    "        dy = segment.coords[1][1] - segment.coords[0][1]\n",
    "        length = np.sqrt(dx**2 + dy**2)\n",
    "        perpendicular_direction = (-dy / length, dx / length)\n",
    "        \n",
    "        # Calculate the start and end points of the perpendicular line\n",
    "        start_point = Point(point.x + half_length * perpendicular_direction[0],\n",
    "                            point.y + half_length * perpendicular_direction[1])\n",
    "        end_point = Point(point.x - half_length * perpendicular_direction[0],\n",
    "                          point.y - half_length * perpendicular_direction[1])\n",
    "        \n",
    "        # Create the perpendicular line and add it to the list\n",
    "        perpendicular_line = LineString([start_point, end_point])\n",
    "        perpendicular_lines.append(perpendicular_line)\n",
    "    \n",
    "    # Create a GeoDataFrame from the perpendicular lines\n",
    "    channel_belt_cross_sections = gpd.GeoDataFrame({'geometry': perpendicular_lines}, crs=simplified_path.crs)\n",
    "    \n",
    "    return channel_belt_cross_sections\n",
    "\n",
    "# Function to calculate channel count index\n",
    "def calc_channel_count_index(filtered_links, cross_sections):\n",
    "    \"\"\"\n",
    "    Calculates the Channel Count Index (CCI) for a network of links intersecting with cross sections.\n",
    "\n",
    "    Parameters:\n",
    "    filtered_links (geopandas.GeoDataFrame): The GeoDataFrame containing the network of links (line segments) with a 'chnl_cat' classification.\n",
    "    cross_sections (geopandas.GeoDataFrame): The GeoDataFrame containing the cross sections.\n",
    "\n",
    "    Returns:\n",
    "    tuple: A tuple containing:\n",
    "        - cci (float): The Channel Count Index.\n",
    "        - cross_sections (geopandas.GeoDataFrame): The cross sections GeoDataFrame with an additional 'channel_count' column.\n",
    "    \"\"\"\n",
    "    channel_counts = []\n",
    "    for idx, cross_section in cross_sections.iterrows():\n",
    "        cross_section_geom = cross_section.geometry\n",
    "        main_channel_count = filtered_links[(filtered_links['chnl_cat'] == 'main_channel') & (filtered_links.intersects(cross_section_geom))].shape[0]\n",
    "        total_count = filtered_links[filtered_links.intersects(cross_section_geom)].shape[0]\n",
    "        if main_channel_count > 1:\n",
    "            total_count -= (main_channel_count - 1)\n",
    "        channel_counts.append(total_count)\n",
    "    cross_sections['channel_count'] = channel_counts\n",
    "    cci = sum(channel_counts) / len(channel_counts)\n",
    "    print(f\"Channel Count Index (CCI): {cci}\")\n",
    "    return cci, cross_sections\n",
    "\n",
    "# Function to calculate sinuosity\n",
    "def calc_sinuosity(shortest_path, simplified_path):\n",
    "    \"\"\"\n",
    "    Calculates the sinuosity of a path by comparing the lengths of the shortest path and the simplified path.\n",
    "\n",
    "    Parameters:\n",
    "    shortest_path (geopandas.GeoDataFrame): The GeoDataFrame containing the shortest path as a LineString.\n",
    "    simplified_path (geopandas.GeoDataFrame): The GeoDataFrame containing the simplified path as a LineString.\n",
    "\n",
    "    Returns:\n",
    "    float: The sinuosity value, which is the ratio of the shortest path length to the simplified path length.\n",
    "    \"\"\"\n",
    "    shortest_path_line = shortest_path.geometry.iloc[0]\n",
    "    simplified_path_line = simplified_path.geometry.iloc[0]\n",
    "    shortest_path_length = shortest_path_line.length\n",
    "    simplified_path_length = simplified_path_line.length\n",
    "    sinuosity = shortest_path_length / simplified_path_length\n",
    "    print(f\"Sinuosity: {sinuosity}\")\n",
    "    return sinuosity\n",
    "\n",
    "# Function to calculate channel form index\n",
    "def calculate_channel_form_index(sinuosity, cci):\n",
    "    \"\"\"\n",
    "    Calculates the Channel Form Index (CFI) based on sinuosity and Channel Count Index (CCI).\n",
    "\n",
    "    Parameters:\n",
    "    sinuosity (float): The sinuosity of the channel.\n",
    "    cci (float): The Channel Count Index.\n",
    "\n",
    "    Returns:\n",
    "    float: The Channel Form Index (CFI).\n",
    "    \"\"\"\n",
    "    cfi = sinuosity / cci\n",
    "    print(f\"Channel Form Index (CFI): {cfi}\")\n",
    "    return cfi\n",
    "\n",
    "# Main function to process network\n",
    "def process_network_folder(river, \n",
    "                           radius, \n",
    "                           year_range=\"All\", \n",
    "                           reach_range=\"All\", \n",
    "                           num_lines=10, \n",
    "                           num_vertices=10, \n",
    "                           fraction_length=1/5, \n",
    "                           root_input=\"C:/Users/huckr/Desktop/UCSB/Dissertation/Data/RiverMapping/RiverMasks\", \n",
    "                           root_output=\"C:/Users/huckr/Desktop/UCSB/Dissertation/Data/RiverMapping/Channels\"):\n",
    "    \"\"\"\n",
    "    Processes a folder containing water mask rasters to extract river channel networks and calculate metrics.\n",
    "\n",
    "    Parameters:\n",
    "    river (str): Name of the river.\n",
    "    radius (int): Radius for conditional dilation.\n",
    "    year_range (tuple or str): Year range for processing (default is \"All\").\n",
    "    reach_range (tuple or str): Reach range for processing (default is \"All\").\n",
    "    num_lines (int): Number of perpendicular lines (cross-sections) (default is 10).\n",
    "    num_vertices (int): Number of vertices for simplifying the shortest path (default is 10).\n",
    "    fraction_length (float): Fraction length for creating cross-sections (default is 1/5).\n",
    "    root_input (str): Root input directory (default is the specified path).\n",
    "    root_output (str): Root output directory (default is the specified path).\n",
    "\n",
    "    Returns:\n",
    "    None\n",
    "    \"\"\"\n",
    "    input_folder = os.path.join(root_input, river)\n",
    "    output_folder_base = os.path.join(root_output, river)\n",
    "    \n",
    "    if year_range == \"All\":\n",
    "        year_range = (1000, 9999)  # Arbitrary wide range to include all years\n",
    "    if reach_range == \"All\":\n",
    "        reach_range = (1, 9999)  # Arbitrary wide range to include all reaches\n",
    "    \n",
    "    year_start, year_end = year_range\n",
    "    reach_start, reach_end = reach_range\n",
    "    \n",
    "    # Initialize a dictionary to store metrics\n",
    "    metrics = {}\n",
    "\n",
    "    # Process each reach folder\n",
    "    for reach_folder in glob(os.path.join(input_folder, 'reach_*')):\n",
    "        reach_folder_name = os.path.basename(reach_folder)\n",
    "        match_reach = re.match(r\"reach_(\\d+)\", reach_folder_name)\n",
    "        if match_reach:\n",
    "            reach = int(match_reach.group(1))\n",
    "            if reach_start <= reach <= reach_end:\n",
    "                for file_path in glob(os.path.join(reach_folder, '*.tif')):\n",
    "                    file_name = os.path.basename(file_path)\n",
    "                    match_year = re.match(rf\"{river}_reach_{reach}_(\\d{{4}})_.*\\.tif\", file_name)\n",
    "                    if match_year:\n",
    "                        year = int(match_year.group(1))\n",
    "                        if year_start <= year <= year_end:\n",
    "                            output_folder = os.path.join(output_folder_base, f\"reach_{reach}\", str(year))\n",
    "                            os.makedirs(output_folder, exist_ok=True)\n",
    "                            \n",
    "                            try:\n",
    "                                water_mask, metadata = load_raster(file_path)\n",
    "                                skeleton = skeletonize(water_mask > 0)\n",
    "                                dilated_skeleton = conditional_dilation(skeleton, radius)\n",
    "                                reskeletonized = skeletonize(dilated_skeleton > 0)\n",
    "                                largest_component = keep_largest_component(reskeletonized)\n",
    "                                \n",
    "                                largest_component_output_path = os.path.join(output_folder, 'largest_component.tif')\n",
    "                                save_raster(largest_component_output_path, largest_component, metadata)\n",
    "\n",
    "                                node_points = create_nodes(largest_component, metadata)\n",
    "                                initial_links = create_links(largest_component, metadata)\n",
    "                                filtered_links = filter_links(initial_links)\n",
    "                                start_end_pts = find_furthest_endpoints(node_points)\n",
    "                                pruned_links = prune_network(node_points, filtered_links, start_end_pts)\n",
    "                                shortest_path_gdf = find_shortest_path(start_end_pts, pruned_links)\n",
    "                                classified_links = classify_channels(pruned_links, shortest_path_gdf)\n",
    "                                valley_center_line = simplify_shortest_path(shortest_path_gdf, num_vertices)\n",
    "                                channel_belt_cross_sections = create_perpendicular_lines(valley_center_line, num_lines, fraction_length)\n",
    "                                \n",
    "                                sinuosity_value = calc_sinuosity(shortest_path_gdf, valley_center_line)\n",
    "                                cci, updated_cross_sections = calc_channel_count_index(classified_links, channel_belt_cross_sections)\n",
    "                                cfi_value = calculate_channel_form_index(sinuosity_value, cci)\n",
    "                                \n",
    "                                classified_links.to_file(os.path.join(output_folder, 'channel_links.shp'))\n",
    "                                channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n",
    "                                node_points.to_file(os.path.join(output_folder, 'nodes.shp'))\n",
    "                                shortest_path_gdf.to_file(os.path.join(output_folder, 'main_channel.shp'))\n",
    "                                valley_center_line.to_file(os.path.join(output_folder, 'valley_center_line.shp'))\n",
    "                                \n",
    "                                # Store metrics\n",
    "                                reach_key = f\"reach_{reach}\"\n",
    "                                if reach_key not in metrics:\n",
    "                                    metrics[reach_key] = {}\n",
    "                                metrics[reach_key][year] = {\n",
    "                                    'Sinuosity': sinuosity_value,\n",
    "                                    'CCI': cci,\n",
    "                                    'CFI': cfi_value\n",
    "                                }\n",
    "                            except Exception as e:\n",
    "                                logging.error(f\"Error processing file {file_path}: {e}\")\n",
    "                                continue\n",
    "\n",
    "    # Save metrics to an Excel workbook\n",
    "    metrics_output_path = os.path.join(root_output, f'{river}_metrics.xlsx')\n",
    "    with pd.ExcelWriter(metrics_output_path) as writer:\n",
    "        for reach, reach_metrics in metrics.items():\n",
    "            df = pd.DataFrame.from_dict(reach_metrics, orient='index')\n",
    "            df.to_excel(writer, sheet_name=reach)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fe48584",
   "metadata": {},
   "source": [
    "## Execute code for a river, a reach, or specific years"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "df9d3ce3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.2372219337513015\n",
      "Channel Count Index (CCI): 10.6\n",
      "Channel Form Index (CFI): 0.11671905035389638\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.2335158798187458\n",
      "Channel Count Index (CCI): 12.6\n",
      "Channel Form Index (CFI): 0.09789808569990047\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.2271631251943798\n",
      "Channel Count Index (CCI): 7.4\n",
      "Channel Form Index (CFI): 0.16583285475599727\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.215190930048138\n",
      "Channel Count Index (CCI): 11.9\n",
      "Channel Form Index (CFI): 0.10211688487799479\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.2134020148285742\n",
      "Channel Count Index (CCI): 14.8\n",
      "Channel Form Index (CFI): 0.08198662262355232\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.1976217188575813\n",
      "Channel Count Index (CCI): 9.9\n",
      "Channel Form Index (CFI): 0.12097189079369508\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.2085707840791324\n",
      "Channel Count Index (CCI): 13.8\n",
      "Channel Form Index (CFI): 0.08757759304921249\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.2223609987470696\n",
      "Channel Count Index (CCI): 12.1\n",
      "Channel Form Index (CFI): 0.10102157014438592\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.2033042420572442\n",
      "Channel Count Index (CCI): 7.3\n",
      "Channel Form Index (CFI): 0.16483619754208825\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.2290428228541121\n",
      "Channel Count Index (CCI): 10.9\n",
      "Channel Form Index (CFI): 0.11275622228019377\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.2276142973623139\n",
      "Channel Count Index (CCI): 6.7\n",
      "Channel Form Index (CFI): 0.18322601453168863\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.24888923546072\n",
      "Channel Count Index (CCI): 7.7\n",
      "Channel Form Index (CFI): 0.1621934072026909\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.22747542293799\n",
      "Channel Count Index (CCI): 12.1\n",
      "Channel Form Index (CFI): 0.10144424982958596\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.1961433827704657\n",
      "Channel Count Index (CCI): 11.3\n",
      "Channel Form Index (CFI): 0.10585339670535095\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.2311470043903034\n",
      "Channel Count Index (CCI): 8.3\n",
      "Channel Form Index (CFI): 0.1483309643843739\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.210467265652506\n",
      "Channel Count Index (CCI): 12.0\n",
      "Channel Form Index (CFI): 0.10087227213770883\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sinuosity: 1.236727570460706\n",
      "Channel Count Index (CCI): 12.4\n",
      "Channel Form Index (CFI): 0.09973609439199241\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\huckr\\AppData\\Local\\Temp\\ipykernel_27416\\2859813101.py:494: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
      "  channel_belt_cross_sections.to_file(os.path.join(output_folder, 'channel_belt_cross_sections.shp'))\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[3], line 25\u001b[0m\n\u001b[0;32m     22\u001b[0m fraction_length\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m\u001b[38;5;241m/\u001b[39m\u001b[38;5;241m3\u001b[39m\n\u001b[0;32m     24\u001b[0m \u001b[38;5;66;03m# Process network\u001b[39;00m\n\u001b[1;32m---> 25\u001b[0m process_network_folder(river, \n\u001b[0;32m     26\u001b[0m                            radius, \n\u001b[0;32m     27\u001b[0m                            year_range, \n\u001b[0;32m     28\u001b[0m                            reach_range, \n\u001b[0;32m     29\u001b[0m                            num_lines, \n\u001b[0;32m     30\u001b[0m                            num_vertices, \n\u001b[0;32m     31\u001b[0m                            fraction_length, \n\u001b[0;32m     32\u001b[0m                            root_input, \n\u001b[0;32m     33\u001b[0m                            root_output)\n",
      "Cell \u001b[1;32mIn[2], line 483\u001b[0m, in \u001b[0;36mprocess_network_folder\u001b[1;34m(river, radius, year_range, reach_range, num_lines, num_vertices, fraction_length, root_input, root_output)\u001b[0m\n\u001b[0;32m    481\u001b[0m filtered_links \u001b[38;5;241m=\u001b[39m filter_links(initial_links)\n\u001b[0;32m    482\u001b[0m start_end_pts \u001b[38;5;241m=\u001b[39m find_furthest_endpoints(node_points)\n\u001b[1;32m--> 483\u001b[0m pruned_links \u001b[38;5;241m=\u001b[39m prune_network(node_points, filtered_links, start_end_pts)\n\u001b[0;32m    484\u001b[0m shortest_path_gdf \u001b[38;5;241m=\u001b[39m find_shortest_path(start_end_pts, pruned_links)\n\u001b[0;32m    485\u001b[0m classified_links \u001b[38;5;241m=\u001b[39m classify_channels(pruned_links, shortest_path_gdf)\n",
      "Cell \u001b[1;32mIn[2], line 262\u001b[0m, in \u001b[0;36mprune_network\u001b[1;34m(nodes, filtered_links, start_end_pts)\u001b[0m\n\u001b[0;32m    260\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m idx, p1 \u001b[38;5;129;01min\u001b[39;00m small_ends\u001b[38;5;241m.\u001b[39miterrows():\n\u001b[0;32m    261\u001b[0m     pn1 \u001b[38;5;241m=\u001b[39m nodes[nodes\u001b[38;5;241m.\u001b[39mgeometry \u001b[38;5;241m==\u001b[39m p1\u001b[38;5;241m.\u001b[39mgeometry]\u001b[38;5;241m.\u001b[39miloc[\u001b[38;5;241m0\u001b[39m]\n\u001b[1;32m--> 262\u001b[0m     nodes_excluding_pn1 \u001b[38;5;241m=\u001b[39m nodes[nodes\u001b[38;5;241m.\u001b[39mgeometry \u001b[38;5;241m!=\u001b[39m pn1\u001b[38;5;241m.\u001b[39mgeometry]\n\u001b[0;32m    263\u001b[0m     nearest_node_geom \u001b[38;5;241m=\u001b[39m nearest_points(pn1\u001b[38;5;241m.\u001b[39mgeometry, nodes_excluding_pn1\u001b[38;5;241m.\u001b[39munary_union)[\u001b[38;5;241m1\u001b[39m]\n\u001b[0;32m    264\u001b[0m     nearest_node \u001b[38;5;241m=\u001b[39m nodes[nodes\u001b[38;5;241m.\u001b[39mgeometry \u001b[38;5;241m==\u001b[39m nearest_node_geom]\u001b[38;5;241m.\u001b[39miloc[\u001b[38;5;241m0\u001b[39m]\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\pandas\\core\\ops\\common.py:76\u001b[0m, in \u001b[0;36m_unpack_zerodim_and_defer.<locals>.new_method\u001b[1;34m(self, other)\u001b[0m\n\u001b[0;32m     72\u001b[0m             \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mNotImplemented\u001b[39m\n\u001b[0;32m     74\u001b[0m other \u001b[38;5;241m=\u001b[39m item_from_zerodim(other)\n\u001b[1;32m---> 76\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m method(\u001b[38;5;28mself\u001b[39m, other)\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\pandas\\core\\arraylike.py:44\u001b[0m, in \u001b[0;36mOpsMixin.__ne__\u001b[1;34m(self, other)\u001b[0m\n\u001b[0;32m     42\u001b[0m \u001b[38;5;129m@unpack_zerodim_and_defer\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m__ne__\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m     43\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__ne__\u001b[39m(\u001b[38;5;28mself\u001b[39m, other):\n\u001b[1;32m---> 44\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_cmp_method(other, operator\u001b[38;5;241m.\u001b[39mne)\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\pandas\\core\\series.py:6119\u001b[0m, in \u001b[0;36mSeries._cmp_method\u001b[1;34m(self, other, op)\u001b[0m\n\u001b[0;32m   6116\u001b[0m lvalues \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_values\n\u001b[0;32m   6117\u001b[0m rvalues \u001b[38;5;241m=\u001b[39m extract_array(other, extract_numpy\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m, extract_range\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[1;32m-> 6119\u001b[0m res_values \u001b[38;5;241m=\u001b[39m ops\u001b[38;5;241m.\u001b[39mcomparison_op(lvalues, rvalues, op)\n\u001b[0;32m   6121\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_construct_result(res_values, name\u001b[38;5;241m=\u001b[39mres_name)\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\pandas\\core\\ops\\array_ops.py:330\u001b[0m, in \u001b[0;36mcomparison_op\u001b[1;34m(left, right, op)\u001b[0m\n\u001b[0;32m    321\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[0;32m    322\u001b[0m             \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mLengths must match to compare\u001b[39m\u001b[38;5;124m\"\u001b[39m, lvalues\u001b[38;5;241m.\u001b[39mshape, rvalues\u001b[38;5;241m.\u001b[39mshape\n\u001b[0;32m    323\u001b[0m         )\n\u001b[0;32m    325\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m should_extension_dispatch(lvalues, rvalues) \u001b[38;5;129;01mor\u001b[39;00m (\n\u001b[0;32m    326\u001b[0m     (\u001b[38;5;28misinstance\u001b[39m(rvalues, (Timedelta, BaseOffset, Timestamp)) \u001b[38;5;129;01mor\u001b[39;00m right \u001b[38;5;129;01mis\u001b[39;00m NaT)\n\u001b[0;32m    327\u001b[0m     \u001b[38;5;129;01mand\u001b[39;00m lvalues\u001b[38;5;241m.\u001b[39mdtype \u001b[38;5;241m!=\u001b[39m \u001b[38;5;28mobject\u001b[39m\n\u001b[0;32m    328\u001b[0m ):\n\u001b[0;32m    329\u001b[0m     \u001b[38;5;66;03m# Call the method on lvalues\u001b[39;00m\n\u001b[1;32m--> 330\u001b[0m     res_values \u001b[38;5;241m=\u001b[39m op(lvalues, rvalues)\n\u001b[0;32m    332\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m is_scalar(rvalues) \u001b[38;5;129;01mand\u001b[39;00m isna(rvalues):  \u001b[38;5;66;03m# TODO: but not pd.NA?\u001b[39;00m\n\u001b[0;32m    333\u001b[0m     \u001b[38;5;66;03m# numpy does not like comparisons vs None\u001b[39;00m\n\u001b[0;32m    334\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m op \u001b[38;5;129;01mis\u001b[39;00m operator\u001b[38;5;241m.\u001b[39mne:\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\geopandas\\array.py:1513\u001b[0m, in \u001b[0;36mGeometryArray.__ne__\u001b[1;34m(self, other)\u001b[0m\n\u001b[0;32m   1512\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__ne__\u001b[39m(\u001b[38;5;28mself\u001b[39m, other):\n\u001b[1;32m-> 1513\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_binop(other, operator\u001b[38;5;241m.\u001b[39mne)\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\geopandas\\array.py:1504\u001b[0m, in \u001b[0;36mGeometryArray._binop\u001b[1;34m(self, other, op)\u001b[0m\n\u001b[0;32m   1500\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mLengths must match to compare\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m   1502\u001b[0m \u001b[38;5;66;03m# If the operator is not defined for the underlying objects,\u001b[39;00m\n\u001b[0;32m   1503\u001b[0m \u001b[38;5;66;03m# a TypeError should be raised\u001b[39;00m\n\u001b[1;32m-> 1504\u001b[0m res \u001b[38;5;241m=\u001b[39m [op(a, b) \u001b[38;5;28;01mfor\u001b[39;00m (a, b) \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mzip\u001b[39m(lvalues, rvalues)]\n\u001b[0;32m   1506\u001b[0m res \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39masarray(res, dtype\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mbool\u001b[39m)\n\u001b[0;32m   1507\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m res\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\geopandas\\array.py:1504\u001b[0m, in \u001b[0;36m<listcomp>\u001b[1;34m(.0)\u001b[0m\n\u001b[0;32m   1500\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mLengths must match to compare\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m   1502\u001b[0m \u001b[38;5;66;03m# If the operator is not defined for the underlying objects,\u001b[39;00m\n\u001b[0;32m   1503\u001b[0m \u001b[38;5;66;03m# a TypeError should be raised\u001b[39;00m\n\u001b[1;32m-> 1504\u001b[0m res \u001b[38;5;241m=\u001b[39m [op(a, b) \u001b[38;5;28;01mfor\u001b[39;00m (a, b) \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mzip\u001b[39m(lvalues, rvalues)]\n\u001b[0;32m   1506\u001b[0m res \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39masarray(res, dtype\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mbool\u001b[39m)\n\u001b[0;32m   1507\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m res\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\shapely\\geometry\\base.py:212\u001b[0m, in \u001b[0;36mBaseGeometry.__ne__\u001b[1;34m(self, other)\u001b[0m\n\u001b[0;32m    210\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(other, BaseGeometry):\n\u001b[0;32m    211\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mNotImplemented\u001b[39m\n\u001b[1;32m--> 212\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__eq__\u001b[39m(other)\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\shapely\\geometry\\base.py:206\u001b[0m, in \u001b[0;36mBaseGeometry.__eq__\u001b[1;34m(self, other)\u001b[0m\n\u001b[0;32m    202\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mNotImplemented\u001b[39m\n\u001b[0;32m    203\u001b[0m \u001b[38;5;66;03m# equal_nan=False is the default, but not yet available for older numpy\u001b[39;00m\n\u001b[0;32m    204\u001b[0m \u001b[38;5;66;03m# TODO updated once we require numpy >= 1.19\u001b[39;00m\n\u001b[0;32m    205\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mtype\u001b[39m(other) \u001b[38;5;241m==\u001b[39m \u001b[38;5;28mtype\u001b[39m(\u001b[38;5;28mself\u001b[39m) \u001b[38;5;129;01mand\u001b[39;00m np\u001b[38;5;241m.\u001b[39marray_equal(\n\u001b[1;32m--> 206\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcoords, other\u001b[38;5;241m.\u001b[39mcoords  \u001b[38;5;66;03m# , equal_nan=False\u001b[39;00m\n\u001b[0;32m    207\u001b[0m )\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\shapely\\geometry\\base.py:223\u001b[0m, in \u001b[0;36mBaseGeometry.coords\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    220\u001b[0m \u001b[38;5;129m@property\u001b[39m\n\u001b[0;32m    221\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mcoords\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[0;32m    222\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Access to geometry's coordinates (CoordinateSequence)\"\"\"\u001b[39;00m\n\u001b[1;32m--> 223\u001b[0m     coords_array \u001b[38;5;241m=\u001b[39m shapely\u001b[38;5;241m.\u001b[39mget_coordinates(\u001b[38;5;28mself\u001b[39m, include_z\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhas_z)\n\u001b[0;32m    224\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m CoordinateSequence(coords_array)\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\shapely\\geometry\\base.py:624\u001b[0m, in \u001b[0;36mBaseGeometry.has_z\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    620\u001b[0m \u001b[38;5;129m@property\u001b[39m\n\u001b[0;32m    621\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mhas_z\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[0;32m    622\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"True if the geometry's coordinate sequence(s) have z values (are\u001b[39;00m\n\u001b[0;32m    623\u001b[0m \u001b[38;5;124;03m    3-dimensional)\"\"\"\u001b[39;00m\n\u001b[1;32m--> 624\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mbool\u001b[39m(shapely\u001b[38;5;241m.\u001b[39mhas_z(\u001b[38;5;28mself\u001b[39m))\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\shapely\\decorators.py:77\u001b[0m, in \u001b[0;36mmultithreading_enabled.<locals>.wrapped\u001b[1;34m(*args, **kwargs)\u001b[0m\n\u001b[0;32m     75\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m arr \u001b[38;5;129;01min\u001b[39;00m array_args:\n\u001b[0;32m     76\u001b[0m         arr\u001b[38;5;241m.\u001b[39mflags\u001b[38;5;241m.\u001b[39mwriteable \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mFalse\u001b[39;00m\n\u001b[1;32m---> 77\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m func(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n\u001b[0;32m     78\u001b[0m \u001b[38;5;28;01mfinally\u001b[39;00m:\n\u001b[0;32m     79\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m arr, old_flag \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mzip\u001b[39m(array_args, old_flags):\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\shapely\\predicates.py:68\u001b[0m, in \u001b[0;36mhas_z\u001b[1;34m(geometry, **kwargs)\u001b[0m\n\u001b[0;32m     41\u001b[0m \u001b[38;5;129m@multithreading_enabled\u001b[39m\n\u001b[0;32m     42\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mhas_z\u001b[39m(geometry, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[0;32m     43\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Returns True if a geometry has a Z coordinate.\u001b[39;00m\n\u001b[0;32m     44\u001b[0m \n\u001b[0;32m     45\u001b[0m \u001b[38;5;124;03m    Note that this function returns False if the (first) Z coordinate equals NaN or\u001b[39;00m\n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m     66\u001b[0m \u001b[38;5;124;03m    False\u001b[39;00m\n\u001b[0;32m     67\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[1;32m---> 68\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m lib\u001b[38;5;241m.\u001b[39mhas_z(geometry, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# Required user inputs\n",
    "\n",
    "# river is the name of the river. Needs to match the name of the folder in which the mask resides\n",
    "river = \"Brahmaputra\"\n",
    "\n",
    "# Specify the range in years to process. The user can enter either a range of years, a specific year, or \"All\" to process\n",
    "# all years for the given reach/reaches. Example inputs are (1997, 2017) for a range of years, 2017 for a single year, or \n",
    "# \"All\" for all years\n",
    "year_range = \"All\"\n",
    "\n",
    "# Specify reach/reaches to process. The user can enter either a range of reaches, a single reach, or \"All\" to process all\n",
    "# reaches for the given river. Example inputs are (1, 40) for a range of reaches, 8 for a single reach, or \n",
    "# \"All\" for all reaches\n",
    "reach_range = (25, 28)\n",
    "\n",
    "# Optional user inputs\n",
    "root_input = \"C:/Users/huckr/Desktop/UCSB/Dissertation/Data/RiverMapping/RiverMasks\"\n",
    "root_output = \"C:/Users/huckr/Desktop/UCSB/Dissertation/Data/RiverMapping/Channels\"\n",
    "radius = 3\n",
    "num_lines=10\n",
    "num_vertices=5\n",
    "fraction_length=1/3\n",
    "\n",
    "# Process network\n",
    "process_network_folder(river, \n",
    "                           radius, \n",
    "                           year_range, \n",
    "                           reach_range, \n",
    "                           num_lines, \n",
    "                           num_vertices, \n",
    "                           fraction_length, \n",
    "                           root_input, \n",
    "                           root_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0728eab1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}