Initial commit

README.md

@@ -1,3 +1,61 @@
 # pyhld
 
 Fiber to the home high level autodesign in Python using shapefiles for input/output.
+
+## USAGE
+
+Clone the repo
+
+git clone https://gitea.bomar.cloud/OpticalFlyer/pyhld.git
+
+Make all scripts executable if not already, i.e.
+
+cd ~/scripts
+chmod +x *.py
+
+Add scripts folder to path:
+
+nano ~/.zshrc
+
+Add the following to end of file, save and close:
+
+export PATH="/Users/youruser/path/to/scripts:$PATH"
+
+Apply the changes:
+
+source ~/.zshrc
+
+Create a folder for your HLD project
+Save sites to a shape file called home_points.shp.  Use a projected CRS.
+
+Run centerlines_from_homes.py from your project directory containing home_points.shp
+
+Run drops_split_centerlines.py - need progress
+
+If the design will have aerial path and you have pole data, drop the pole shape file into your project folder and name it poles.shp
+
+Run connect_poles.py if you need to generate aerial spans.  This connects all poles together using a minimum spanning tree.
+
+Run create_aerial_drops.py
+
+Run create_aerial_edges.py
+
+Run create_transitions.py
+
+Run create_nodes.py - need progress
+
+Run associate_drop_points_v2.py
+
+Run cluster_fdh_v2.py
+
+Run create_network_v2.py
+
+***Make Manual Revisions Here***
+
+Run create_mst_clusters.py - need progress (v2 not ready)
+
+Run poles_used.py - add progress indicator
+
+Run report.py
+
+Run create_map.py

associate_drop_points.py

@@ -0,0 +1,51 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import pandas as pd
+from shapely.ops import nearest_points
+import sys
+
+print("Associating home points to drop points...")
+
+# Load the NODES and HOME_POINT shapefiles
+nodes_gdf = gpd.read_file('nodes.shp')
+home_points_gdf = gpd.read_file('home_points.shp')
+
+# Check if 'drop_point' column exists in home_points_gdf, if not add it
+if 'drop_point' not in home_points_gdf.columns:
+    home_points_gdf['drop_point'] = pd.NA
+
+# Check if 'type' column exists in nodes_gdf, if not add it
+if 'type' not in nodes_gdf.columns:
+    nodes_gdf['type'] = None  # Initialize with None (or use pd.NA for pandas >= 1.0)
+
+total_home_points = len(home_points_gdf)
+
+# Iterating through each home point
+for hp_index, home_point in home_points_gdf.iterrows():
+    # Create a temporary DataFrame with distances to each node
+    distances = nodes_gdf.geometry.apply(lambda node: home_point.geometry.distance(node))
+
+    # Find the index of the closest node
+    closest_node_index = distances.idxmin()
+
+    # Get the ID of the closest node
+    closest_node_id = nodes_gdf.iloc[closest_node_index]['id']
+
+    # Update 'drop_point' for the home point
+    home_points_gdf.at[hp_index, 'drop_point'] = closest_node_id
+
+    # Mark the node as associated with a home point by setting its 'type' to 'HP'
+    nodes_gdf.at[closest_node_index, 'type'] = 'HP'
+
+    # Print progress indicator
+    print(f"Processing home point {hp_index + 1}/{total_home_points}...", end='\r')
+    sys.stdout.flush()  # Ensure the progress text is displayed immediately
+
+# Save the updated home_points_gdf back to a shapefile, overwriting if it exists
+home_points_gdf.to_file('home_points.shp', overwrite=True)
+
+# Also, save the updated nodes_gdf with the 'type' attribute
+nodes_gdf.to_file('nodes.shp', overwrite=True)
+
+print("\nDone.")

associate_drop_points_v2.py

@@ -0,0 +1,41 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import pandas as pd
+import numpy as np
+
+print("Associating home points to drop points and updating node types...")
+
+# Load the NODES and HOME_POINT shapefiles
+nodes_gdf = gpd.read_file('nodes.shp')
+home_points_gdf = gpd.read_file('home_points.shp')
+
+# Ensure 'drop_point' and 'type' columns exist
+home_points_gdf['drop_point'] = pd.NA
+nodes_gdf['type'] = None
+
+def find_nearest_node(home_point, nodes):
+    # Calculate distances from the home point to all nodes
+    distances = nodes.geometry.distance(home_point.geometry)
+    # Find the index of the minimum distance
+    closest_node_index = distances.idxmin()
+    # Return the ID of the closest node
+    return nodes.iloc[closest_node_index].id
+
+# Vectorize the find_nearest_node function if possible or use apply as a fallback
+try:
+    # Attempt to use a vectorized approach for performance
+    home_points_gdf['drop_point'] = np.vectorize(find_nearest_node)(home_points_gdf.geometry, nodes_gdf)
+except Exception as e:
+    # Fallback to using apply if the vectorized approach fails
+    home_points_gdf['drop_point'] = home_points_gdf.apply(lambda x: find_nearest_node(x, nodes_gdf), axis=1)
+
+# Update nodes 'type' based on associated home points
+unique_drop_points = home_points_gdf['drop_point'].unique()
+nodes_gdf.loc[nodes_gdf['id'].isin(unique_drop_points), 'type'] = 'HP'
+
+# Save the updated GeoDataFrames back to shapefiles
+home_points_gdf.to_file('home_points.shp', overwrite=True)
+nodes_gdf.to_file('nodes.shp', overwrite=True)
+
+print("Done.")

centerlines_from_homes.py

@@ -0,0 +1,90 @@
+#!/usr/bin/env python3
+
+import osmnx as ox
+import geopandas as gpd
+from shapely.geometry import LineString, box
+
+print("Getting road centerlines...")
+
+# Load home points from a shapefile
+home_points_gdf = gpd.read_file('home_points.shp')
+
+# Check if there's at least one feature in the DataFrame
+if not home_points_gdf.empty:
+    # Check if CRS is projected, if not, prompt to check the CRS
+    if home_points_gdf.crs.is_geographic:
+        print("Warning: The CRS of the shapefile is geographic. Buffering might be inaccurate.")
+    else:
+        print(f"Using projected CRS {home_points_gdf.crs} for buffering.")
+
+    # Determine the units of the CRS
+    crs_units = home_points_gdf.crs.axis_info[0].unit_name
+    print(f"CRS units: {crs_units}")
+
+    # Calculate the total bounds of the home points and extend by 1000 feet on each side
+    extension_distance_feet = 1000  # Distance in feet
+    extension_distance = extension_distance_feet  # Default in feet
+
+    # Convert the extension distance from feet to the CRS units if necessary
+    if crs_units == 'meter':
+        extension_distance = extension_distance_feet * 0.3048  # Convert feet to meters
+
+    minx, miny, maxx, maxy = home_points_gdf.total_bounds
+    extended_bounds = (minx - extension_distance, miny - extension_distance,
+                       maxx + extension_distance, maxy + extension_distance)
+
+    # Create an extended bounding box around the home points
+    bbox = box(*extended_bounds)
+    
+    # Reproject the bounding box to EPSG:4326 for OSMnx
+    bbox_reprojected = gpd.GeoSeries([bbox], crs=home_points_gdf.crs).to_crs('epsg:4326')
+
+    # Use OSMnx to download the street network
+    G = ox.graph_from_polygon(bbox_reprojected.iloc[0], network_type='drive')
+
+    # Convert the graph into GeoDataFrames
+    gdf_nodes, gdf_edges = ox.graph_to_gdfs(G)
+
+    # Check each row in each column and convert non-compatible types to strings, skip geometry column
+    for column in gdf_edges.columns:
+        if column != 'geometry':  # Skip the geometry column
+            for idx, value in gdf_edges[column].items():
+                if isinstance(value, (list, dict, set)):
+                    gdf_edges.at[idx, column] = ', '.join(map(str, value))
+                elif not isinstance(value, (int, float, str, type(None))):
+                    gdf_edges.at[idx, column] = str(value)
+            print(f"Processed column: {column}")
+
+    # Reproject gdf_edges back to the original CRS of home points
+    original_crs = home_points_gdf.crs
+    gdf_edges = gdf_edges.to_crs(original_crs)
+
+    print("Normalizing geometries...")
+
+    # Normalize geometries and deduplicate
+    def normalize_geometry(geometry):
+        if isinstance(geometry, LineString):
+            return LineString(sorted(geometry.coords))
+        return geometry
+
+    gdf_edges['normalized_geom'] = gdf_edges['geometry'].apply(normalize_geometry)
+    gdf_edges_deduped = gdf_edges.drop_duplicates(subset=['osmid', 'normalized_geom'])
+
+    # Drop the temporary 'normalized_geom' column
+    gdf_edges_deduped = gdf_edges_deduped.drop(columns=['normalized_geom'])
+
+    # Reproject gdf_edges_deduped back to the original CRS of home points
+    original_crs = home_points_gdf.crs
+    gdf_edges_deduped = gdf_edges_deduped.to_crs(original_crs)
+
+    # Attempt to save the deduplicated DataFrame to a shapefile
+    try:
+        gdf_edges_deduped.to_file('road_centerlines.shp')
+        print(f"Road centerlines have been saved to road_centerlines.shp in CRS {original_crs}.")
+    except Exception as e:
+        print(f"Error encountered: {e}")
+
+else:
+    print("No points found in the home_points.shp file.")
+
+print("Done.")

check_graph.py

@@ -0,0 +1,32 @@
+#!/usr/bin/env python3
+
+import networkx as nx
+import geopandas as gpd
+
+print("Building the graph...")
+
+# Load shapefiles
+edges_gdf = gpd.read_file('edges.shp')
+nodes_gdf = gpd.read_file('nodes.shp')
+
+# Build the graph
+G = nx.Graph()
+for _, edge in edges_gdf.iterrows():
+    G.add_edge(edge['start_node'], edge['end_node'], weight=edge['cost'], type=edge['type'], length=edge['length'], cost=edge['cost'])
+
+# Get connected components
+components = list(nx.connected_components(G))
+
+# Check if the graph is fully connected
+if len(components) == 1:
+    print("The graph is fully connected.")
+else:
+    num_components = len(components)
+    print(f"The graph is not fully connected. It has {num_components} connected components.")
+    for i, component in enumerate(components):
+        print(f"Connected Component {i + 1}: {component}")
+
+# Print nodes that are not connected
+if len(components) > 1:
+    isolated_nodes = set(G.nodes()) - set.union(*components)
+    print(f"Isolated Nodes (not connected to any component): {isolated_nodes}")

cluster_fdh.py

@@ -0,0 +1,76 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import geopandas as gpd
+from scipy.cluster.hierarchy import fcluster, linkage
+from shapely.geometry import Point
+from shapely.ops import nearest_points
+
+def exclude_outliers(data):
+    """Exclude outliers using the IQR method."""
+    if len(data) == 0:
+        return np.array([], dtype=bool)
+    Q1 = np.percentile(data, 25, axis=0)
+    Q3 = np.percentile(data, 75, axis=0)
+    IQR = Q3 - Q1
+    lower_bound = Q1 - 1.5 * IQR
+    upper_bound = Q3 + 1.5 * IQR
+    # Returns a boolean mask where True indicates the points that are not outliers
+    return (data >= lower_bound) & (data <= upper_bound)
+
+def find_median_center(cluster_points):
+    """Find the median center of a cluster, excluding outliers."""
+    if len(cluster_points) == 0:
+        return None, None
+    outlier_mask = exclude_outliers(cluster_points)
+    # Apply mask to filter both dimensions
+    filtered_points = cluster_points[np.all(outlier_mask, axis=1)]
+    median_x = np.median(filtered_points[:, 0])
+    median_y = np.median(filtered_points[:, 1])
+    return median_x, median_y
+
+def snap_to_road(median_center, roads):
+    """Snap median center to the closest point on the closest road line."""
+    closest_road = roads.geometry.unary_union
+    closest_point = nearest_points(median_center, closest_road)[1]
+    return closest_point
+
+def cluster_homes_and_save(shapefile_path, road_centerlines_path='road_centerlines.shp', home_points_path='home_points_fdh.shp', fdh_path='fdh.shp', max_homes_per_cluster=432):
+    homes = gpd.read_file(shapefile_path)
+    roads = gpd.read_file(road_centerlines_path)
+    coordinates = np.array(list(homes.geometry.apply(lambda x: (x.x, x.y))))
+    Z = linkage(coordinates, method='ward')
+
+    # Attempt to find an appropriate distance threshold dynamically
+    max_distance = Z[-1, 2]  # Maximum distance in the linkage matrix
+    distance_threshold = max_distance / 2  # Start with half of the maximum distance
+    clusters = fcluster(Z, t=distance_threshold, criterion='distance')
+    while np.max(np.bincount(clusters)) > max_homes_per_cluster and distance_threshold > 0:
+        distance_threshold *= 0.95  # Gradually decrease the threshold
+        clusters = fcluster(Z, t=distance_threshold, criterion='distance')
+    
+    if distance_threshold <= 0:
+        print("Unable to find a suitable distance threshold to meet the cluster size constraint.")
+        return
+
+    homes['fdh_id'] = clusters
+    homes.to_file(home_points_path, driver='ESRI Shapefile')
+    print(f"Clustered homes saved to {home_points_path}")
+
+    # Calculate and save median centers
+    median_centers = []
+    for cluster_id in np.unique(homes['fdh_id']):
+        cluster_points = np.array(list(homes.loc[homes['fdh_id'] == cluster_id, 'geometry'].apply(lambda p: (p.x, p.y))))
+        median_x, median_y = find_median_center(cluster_points)
+        if median_x is not None and median_y is not None:
+            median_center = Point(median_x, median_y)
+            snapped_center = snap_to_road(median_center, roads)
+            median_centers.append({'geometry': snapped_center, 'id': cluster_id})
+
+    median_centers_gdf = gpd.GeoDataFrame(median_centers, geometry='geometry')
+    median_centers_gdf.crs = homes.crs
+    median_centers_gdf.to_file(fdh_path, driver='ESRI Shapefile')
+    print(f"Median centers saved to {fdh_path}")
+
+# Adjust the call to cluster_homes_and_save as needed
+cluster_homes_and_save('home_points.shp', home_points_path='home_points.shp')

cluster_fdh_v2.py

@@ -0,0 +1,80 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import geopandas as gpd
+from scipy.cluster.hierarchy import fcluster, linkage
+from shapely.geometry import Point
+from shapely.ops import nearest_points
+
+def exclude_outliers(data):
+    """Exclude outliers using the IQR method."""
+    if len(data) == 0:
+        return np.array([], dtype=bool)
+    Q1 = np.percentile(data, 25, axis=0)
+    Q3 = np.percentile(data, 75, axis=0)
+    IQR = Q3 - Q1
+    lower_bound = Q1 - 1.5 * IQR
+    upper_bound = Q3 + 1.5 * IQR
+    return (data >= lower_bound) & (data <= upper_bound)
+
+def find_median_center(cluster_points):
+    """Find the median center of a cluster, excluding outliers."""
+    if len(cluster_points) == 0:
+        return None, None
+    outlier_mask = exclude_outliers(cluster_points)
+    filtered_points = cluster_points[np.all(outlier_mask, axis=1)]
+    median_x = np.median(filtered_points[:, 0])
+    median_y = np.median(filtered_points[:, 1])
+    return median_x, median_y
+
+def snap_to_nearest_node(median_center, nodes):
+    """Snap median center to the closest node where 'type' does not equal 'HP'."""
+    # Filter nodes to exclude those with type 'HP'
+    eligible_nodes = nodes[nodes['type'] != 'HP']
+    # Find the nearest eligible node
+    closest_node = eligible_nodes.geometry.unary_union
+    closest_point = nearest_points(median_center, closest_node)[1]
+    # Get the ID of the closest node
+    closest_node_id = eligible_nodes.loc[eligible_nodes.geometry == closest_point, 'id'].values[0]
+    return closest_point, closest_node_id
+
+def cluster_homes_and_save(shapefile_path, nodes_path='nodes.shp', home_points_path='home_points.shp', fdh_path='fdh.shp', max_homes_per_cluster=432):
+    homes = gpd.read_file(shapefile_path)
+    nodes = gpd.read_file(nodes_path)
+    coordinates = np.array(list(homes.geometry.apply(lambda x: (x.x, x.y))))
+    Z = linkage(coordinates, method='ward')
+
+    max_distance = Z[-1, 2]
+    distance_threshold = max_distance / 2
+    clusters = fcluster(Z, t=distance_threshold, criterion='distance')
+    while np.max(np.bincount(clusters)) > max_homes_per_cluster and distance_threshold > 0:
+        distance_threshold *= 0.95
+        clusters = fcluster(Z, t=distance_threshold, criterion='distance')
+    
+    if distance_threshold <= 0:
+        print("Unable to find a suitable distance threshold to meet the cluster size constraint.")
+        return
+
+    homes['fdh_id'] = clusters
+    homes.to_file(home_points_path, driver='ESRI Shapefile')
+    print(f"Clustered homes saved to {home_points_path}")
+
+    median_centers = []
+    for cluster_id in np.unique(homes['fdh_id']):
+        cluster_points = np.array(list(homes.loc[homes['fdh_id'] == cluster_id, 'geometry'].apply(lambda p: (p.x, p.y))))
+        median_x, median_y = find_median_center(cluster_points)
+        if median_x is not None and median_y is not None:
+            median_center = Point(median_x, median_y)
+            snapped_center, node_id = snap_to_nearest_node(median_center, nodes)
+            median_centers.append({
+                'geometry': snapped_center,
+                'id': cluster_id,
+                'node_id': node_id
+            })
+
+    median_centers_gdf = gpd.GeoDataFrame(median_centers, columns=['geometry', 'id', 'node_id'], crs=homes.crs)
+    median_centers_gdf.to_file(fdh_path, driver='ESRI Shapefile')
+    print(f"Median centers saved to {fdh_path}")
+
+# Adjust the call to cluster_homes_and_save as needed
+cluster_homes_and_save('home_points.shp')

connect_poles.py

@@ -0,0 +1,45 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import networkx as nx
+from scipy.spatial import distance_matrix
+from shapely.geometry import LineString
+import numpy as np
+import pandas as pd
+
+# Load the poles from a shapefile
+poles_gdf = gpd.read_file('poles.shp')
+
+# Extract the points coordinates for the distance matrix
+points = np.array(list(zip(poles_gdf.geometry.x, poles_gdf.geometry.y)))
+
+# Create a complete distance matrix
+dist_matrix = distance_matrix(points, points)
+
+# Create a complete graph from the distance matrix
+G = nx.complete_graph(len(points))
+
+# Add edges to the graph with distances as weights
+for i, node in enumerate(G.nodes()):
+    for j, neighbor in enumerate(G.nodes()):
+        if i != j:  # skip if same node
+            G[node][neighbor]['weight'] = dist_matrix[i][j]
+
+# Compute the minimum spanning tree of the complete graph
+mst = nx.minimum_spanning_tree(G)
+
+# Generate the lines for the MST
+mst_lines = []
+for edge in mst.edges(data=False):
+    point1 = points[edge[0]]
+    point2 = points[edge[1]]
+    line = LineString([point1, point2])
+    mst_lines.append(line)
+
+# Create a GeoDataFrame from the lines
+lines_gdf = gpd.GeoDataFrame(geometry=mst_lines, crs=poles_gdf.crs)
+
+# Save the lines to a new shapefile
+lines_gdf.to_file('pole_lines.shp')
+
+print("Pole lines have been saved to pole_lines.shp.")

create_aerial_drops.py

@@ -0,0 +1,76 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import pandas as pd
+from shapely.geometry import LineString
+import numpy as np
+from scipy.spatial import cKDTree
+
+# Constants
+SEARCH_RADIUS = 200  # Search radius in feet
+COST_PER_FOOT = 1.5  # Cost per foot
+
+# Load the shapefiles
+home_points_gdf = gpd.read_file('home_points.shp')
+poles_gdf = gpd.read_file('poles.shp')
+edges_gdf = gpd.read_file('edges.shp')
+
+# Convert search radius to the CRS units of the poles layer (assuming it's in feet if in a projected CRS)
+# If the CRS is in meters, use a conversion factor from feet to meters
+conversion_factor = 1 # if poles_gdf.crs.axis_info[0].unit_name == 'foot' else 0.3048
+search_radius_in_crs_units = SEARCH_RADIUS * conversion_factor
+
+# Check for invalid geometries
+invalid_geometries = poles_gdf[~poles_gdf.is_valid]
+if not invalid_geometries.empty:
+    print(f"Found {len(invalid_geometries)} invalid geometries. Removing them.")
+    poles_gdf = poles_gdf[poles_gdf.is_valid]
+
+# Ensure all coordinates are finite
+finite_coords_check = np.isfinite(poles_gdf.geometry.x) & np.isfinite(poles_gdf.geometry.y)
+if not finite_coords_check.all():
+    print("Found non-finite coordinates. Removing corresponding entries.")
+    poles_gdf = poles_gdf[finite_coords_check]
+
+# Create KDTree for poles for efficient nearest neighbor search
+poles_coords = np.array(list(zip(poles_gdf.geometry.x, poles_gdf.geometry.y)))
+tree = cKDTree(poles_coords)
+
+# List to store new drop edges
+new_edges = []
+
+# Iterate over every home point
+for idx, home in home_points_gdf.iterrows():
+    # Query the nearest pole within the search radius
+    distance, index = tree.query(home.geometry.coords[0], k=1, distance_upper_bound=search_radius_in_crs_units)
+    if distance != np.inf:
+        # If a pole is found within the search radius, create the new edge
+        nearest_pole = poles_gdf.iloc[index].geometry
+        #line = LineString([home.geometry, nearest_pole])
+        home_coords = (home.geometry.x, home.geometry.y)
+        nearest_pole_coords = (nearest_pole.x, nearest_pole.y)
+        line = LineString([home_coords, nearest_pole_coords])
+        
+        # Calculate cost
+        length = line.length  # Length in CRS units
+        cost = length * conversion_factor * COST_PER_FOOT  # Convert length to feet and calculate cost
+        
+        # Create a new edge feature with the specified attributes
+        new_edge = {
+            'type': 'Aerial Drop',
+            'length': length,
+            'cost': cost,
+            'geometry': line
+        }
+        new_edges.append(new_edge)
+
+# Append new edges to the existing edges GeoDataFrame
+new_edges_gdf = gpd.GeoDataFrame(new_edges, columns=['type', 'length', 'cost', 'geometry'], crs=edges_gdf.crs)
+
+# Concatenate new edges with the existing edges GeoDataFrame
+edges_gdf = pd.concat([edges_gdf, new_edges_gdf], ignore_index=True)
+
+# Save the updated edges GeoDataFrame back to the shapefile
+edges_gdf.to_file('edges.shp')
+
+print("Aerial drops have been added to the edges shapefile.")

create_aerial_edges.py

@@ -0,0 +1,63 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+from shapely.geometry import LineString
+import pandas as pd
+
+# Constants
+COST_PER_UNIT = 2.5  # Assuming the unit is in the same as your CRS
+
+# Load the shapefiles
+pole_lines_gdf = gpd.read_file('pole_lines.shp')
+edges_gdf = gpd.read_file('edges.shp')
+
+# Ensure CRS match if necessary
+# pole_lines_gdf = pole_lines_gdf.to_crs(edges_gdf.crs)
+
+# List to store new edges
+new_edges = []
+
+# Iterate over every feature of the pole_lines_gdf GeoDataFrame
+for idx, pole_line in pole_lines_gdf.iterrows():
+    # Get the geometry of the feature
+    geom = pole_line.geometry
+    
+    # Assume each feature is a LineString (not MultiLineString)
+    if isinstance(geom, LineString):
+        vertices = list(geom.coords)
+        
+        for i in range(1, len(vertices)):
+            # Get two consecutive vertices
+            pt1 = vertices[i - 1]
+            pt2 = vertices[i]
+            
+            # Create a LineString from the vertices
+            line = LineString([pt1, pt2])
+            
+            # Calculate the length of the line
+            length = line.length  # Length is in CRS units
+            
+            # Calculate the cost
+            cost = length * COST_PER_UNIT
+            
+            # Create a new edge feature with attributes and geometry
+            new_edge = {
+                'type': 'Aerial',
+                'length': length,
+                'cost': cost,
+                'geometry': line
+            }
+            
+            # Append the new edge to the list
+            new_edges.append(new_edge)
+
+# Create a GeoDataFrame from the list of new edges
+new_edges_gdf = gpd.GeoDataFrame(new_edges, crs=edges_gdf.crs)
+
+# Concatenate new edges with the existing edges GeoDataFrame
+edges_gdf = pd.concat([edges_gdf, new_edges_gdf], ignore_index=True)
+
+# Save the updated edges GeoDataFrame back to the shapefile
+edges_gdf.to_file('edges.shp')
+
+print("Aerial pole lines have been added to the edges shapefile.")

create_map.py

@@ -0,0 +1,140 @@
+#!/usr/bin/env python3
+
+import folium
+import geopandas as gpd
+from shapely.geometry import box
+import matplotlib.pyplot as plt
+import matplotlib.colors
+import numpy as np
+import os
+
+# Function to check if a file exists and is not empty
+def check_file_exists(file_path):
+    return os.path.exists(file_path) and os.path.getsize(file_path) > 0
+
+# Load the shapefiles
+network_gdf = gpd.read_file('network.shp')
+home_points_gdf = gpd.read_file('home_points.shp')
+if check_file_exists('poles_used.shp'):
+    poles_used_gdf = gpd.read_file('poles_used.shp').to_crs(epsg=4326)
+else:
+    poles_used_gdf = gpd.GeoDataFrame()  # Create an empty GeoDataFrame if the file does not exist
+#poles_used_gdf = gpd.read_file('poles_used.shp')
+mst_gdf = gpd.read_file('mst.shp')
+if check_file_exists('headend.shp'):
+    headend_gdf = gpd.read_file('headend.shp').to_crs(epsg=4326)
+else:
+    headend_gdf = gpd.GeoDataFrame()  # Create an empty GeoDataFrame if the file does not exist
+#headend_gdf = gpd.read_file('headend.shp')
+fdh_gdf = gpd.read_file('fdh.shp').to_crs(epsg=4326)
+
+# Reproject to EPSG:4326 for Folium compatibility
+network_gdf = network_gdf.to_crs(epsg=4326)
+home_points_gdf = home_points_gdf.to_crs(epsg=4326)
+if not poles_used_gdf.empty:
+    poles_used_gdf = poles_used_gdf.to_crs(epsg=4326)
+mst_gdf = mst_gdf.to_crs(epsg=4326)
+if not headend_gdf.empty:
+    headend_gdf = headend_gdf.to_crs(epsg=4326)
+
+# Calculate the center based on the bounding box of network geometries
+center = [network_gdf.total_bounds[1] + (network_gdf.total_bounds[3] - network_gdf.total_bounds[1]) / 2,
+          network_gdf.total_bounds[0] + (network_gdf.total_bounds[2] - network_gdf.total_bounds[0]) / 2]
+
+# Initialize a Folium map centered on the calculated center
+m = folium.Map(location=center, zoom_start=14)
+
+# Function to add lines to the map with specific styles
+def add_lines(gdf, line_color, line_weight, line_dash_array=None):
+    for _, row in gdf.iterrows():
+        points = [[point[1], point[0]] for point in row.geometry.coords]
+        folium.PolyLine(points, color=line_color, weight=line_weight, dash_array=line_dash_array).add_to(m)
+
+# Add lines from network.shp with different styles based on type
+add_lines(network_gdf[network_gdf['type'] == 'Underground'], 'green', 3)
+add_lines(network_gdf[network_gdf['type'] == 'Aerial'], 'blue', 3)
+add_lines(network_gdf[network_gdf['type'] == 'Aerial Drop'], 'cyan', 1)
+add_lines(network_gdf[network_gdf['type'] == 'Buried Drop'], 'orange', 1)
+add_lines(network_gdf[network_gdf['type'] == 'Transition'], 'green', 3, '5,5')
+
+# Function to add filled circles to the map
+def add_filled_circles(gdf, color, fill_color, radius):
+    for _, row in gdf.iterrows():
+        folium.Circle(
+            location=[row.geometry.y, row.geometry.x], 
+            radius=radius,
+            color=color,
+            fill=True,
+            fill_color=fill_color,
+            fill_opacity=1.0  # Ensure the fill is fully opaque
+        ).add_to(m)
+
+def add_fdh_markers(gdf, icon_color, icon_icon):
+    for _, row in gdf.iterrows():
+        folium.Marker(
+            location=[row.geometry.y, row.geometry.x],
+            icon=folium.Icon(color=icon_color, icon=icon_icon),
+            popup=f"FDH Cabinet {row['id']}"  # Optional: Add a popup label to each FDH marker
+        ).add_to(m)
+
+# Add home points as orange filled circles
+#add_filled_circles(home_points_gdf, 'orange', 'orange', 5)
+        
+def add_home_points_by_fdh(gdf, radius):
+    # Generate a color palette with enough colors for each fdh_id
+    unique_fdh_ids = gdf['fdh_id'].unique()
+    # Update: Use recommended method for accessing colormaps in newer versions of Matplotlib
+    color_palette = plt.colormaps['hsv'](np.linspace(0, 1, len(unique_fdh_ids)))
+    
+    # Create a dictionary to map each fdh_id to a color
+    fdh_id_to_color = {fdh_id: color_palette[i] for i, fdh_id in enumerate(unique_fdh_ids)}
+    
+    # Add home points with colors based on their fdh_id
+    for _, row in gdf.iterrows():
+        fdh_id = row['fdh_id']
+        color = matplotlib.colors.to_hex(fdh_id_to_color[fdh_id])  # Convert color to hex format for Folium
+        folium.Circle(
+            location=[row.geometry.y, row.geometry.x],
+            radius=radius,
+            color=color,
+            fill=True,
+            fill_color=color,
+            fill_opacity=1.0
+        ).add_to(m)
+
+# Replace the original call to add home points with the new function
+add_home_points_by_fdh(home_points_gdf, 5)
+
+# Only add poles if poles_used_gdf is not empty
+if not poles_used_gdf.empty:
+    add_filled_circles(poles_used_gdf, 'black', 'white', 3)
+
+# Add used poles as white filled circles, make them smaller
+#add_filled_circles(poles_used_gdf, 'black', 'white', 3)
+
+add_filled_circles(mst_gdf, 'black', 'red', 4)
+
+'''
+for _, row in headend_gdf.iterrows():
+        folium.Marker(
+            location=[row.geometry.y, row.geometry.x],
+            icon=folium.Icon(color='darkblue', icon='cloud'),
+            popup="Headend"
+        ).add_to(m)
+'''
+
+# Only add headend markers if headend_gdf is not empty
+if not headend_gdf.empty:
+    for _, row in headend_gdf.iterrows():
+        folium.Marker(
+            location=[row.geometry.y, row.geometry.x],
+            icon=folium.Icon(color='darkblue', icon='cloud'),
+            popup="Headend"
+        ).add_to(m)
+
+add_fdh_markers(fdh_gdf, 'darkpurple', 'glyphicon-tower')
+
+# Save the map to an HTML file
+m.save('network_map.html')
+
+print("Map has been saved to network_map.html.")

create_map_buried.py

@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+
+import folium
+import geopandas as gpd
+from shapely.geometry import box
+import matplotlib.pyplot as plt
+import matplotlib.colors
+import numpy as np
+
+# Load the shapefiles
+network_gdf = gpd.read_file('network.shp')
+home_points_gdf = gpd.read_file('home_points.shp')
+#poles_used_gdf = gpd.read_file('poles_used.shp')
+mst_gdf = gpd.read_file('mst.shp')
+#headend_gdf = gpd.read_file('headend.shp')
+fdh_gdf = gpd.read_file('fdh.shp').to_crs(epsg=4326)
+
+# Reproject to EPSG:4326 for Folium compatibility
+network_gdf = network_gdf.to_crs(epsg=4326)
+home_points_gdf = home_points_gdf.to_crs(epsg=4326)
+#poles_used_gdf = poles_used_gdf.to_crs(epsg=4326)
+mst_gdf = mst_gdf.to_crs(epsg=4326)
+#headend_gdf = headend_gdf.to_crs(epsg=4326)
+
+# Calculate the center based on the bounding box of network geometries
+center = [network_gdf.total_bounds[1] + (network_gdf.total_bounds[3] - network_gdf.total_bounds[1]) / 2,
+          network_gdf.total_bounds[0] + (network_gdf.total_bounds[2] - network_gdf.total_bounds[0]) / 2]
+
+# Initialize a Folium map centered on the calculated center
+m = folium.Map(location=center, zoom_start=14)
+
+# Function to add lines to the map with specific styles
+def add_lines(gdf, line_color, line_weight, line_dash_array=None):
+    for _, row in gdf.iterrows():
+        points = [[point[1], point[0]] for point in row.geometry.coords]
+        folium.PolyLine(points, color=line_color, weight=line_weight, dash_array=line_dash_array).add_to(m)
+
+# Add lines from network.shp with different styles based on type
+add_lines(network_gdf[network_gdf['type'] == 'Underground'], 'green', 3)
+add_lines(network_gdf[network_gdf['type'] == 'Aerial'], 'blue', 3)
+add_lines(network_gdf[network_gdf['type'] == 'Aerial Drop'], 'cyan', 1)
+add_lines(network_gdf[network_gdf['type'] == 'Buried Drop'], 'orange', 1)
+add_lines(network_gdf[network_gdf['type'] == 'Transition'], 'green', 3, '5,5')
+
+# Function to add filled circles to the map
+def add_filled_circles(gdf, color, fill_color, radius):
+    for _, row in gdf.iterrows():
+        folium.Circle(
+            location=[row.geometry.y, row.geometry.x], 
+            radius=radius,
+            color=color,
+            fill=True,
+            fill_color=fill_color,
+            fill_opacity=1.0  # Ensure the fill is fully opaque
+        ).add_to(m)
+
+def add_fdh_markers(gdf, icon_color, icon_icon):
+    for _, row in gdf.iterrows():
+        folium.Marker(
+            location=[row.geometry.y, row.geometry.x],
+            icon=folium.Icon(color=icon_color, icon=icon_icon),
+            popup=f"FDH Cabinet {row['id']}"  # Optional: Add a popup label to each FDH marker
+        ).add_to(m)
+
+# Add home points as orange filled circles
+#add_filled_circles(home_points_gdf, 'orange', 'orange', 5)
+        
+def add_home_points_by_fdh(gdf, radius):
+    # Generate a color palette with enough colors for each fdh_id
+    unique_fdh_ids = gdf['fdh_id'].unique()
+    # Update: Use recommended method for accessing colormaps in newer versions of Matplotlib
+    color_palette = plt.colormaps['hsv'](np.linspace(0, 1, len(unique_fdh_ids)))
+    
+    # Create a dictionary to map each fdh_id to a color
+    fdh_id_to_color = {fdh_id: color_palette[i] for i, fdh_id in enumerate(unique_fdh_ids)}
+    
+    # Add home points with colors based on their fdh_id
+    for _, row in gdf.iterrows():
+        fdh_id = row['fdh_id']
+        color = matplotlib.colors.to_hex(fdh_id_to_color[fdh_id])  # Convert color to hex format for Folium
+        folium.Circle(
+            location=[row.geometry.y, row.geometry.x],
+            radius=radius,
+            color=color,
+            fill=True,
+            fill_color=color,
+            fill_opacity=1.0
+        ).add_to(m)
+
+# Replace the original call to add home points with the new function
+add_home_points_by_fdh(home_points_gdf, 5)
+
+# Add used poles as white filled circles, make them smaller
+#add_filled_circles(poles_used_gdf, 'black', 'white', 3)
+
+add_filled_circles(mst_gdf, 'black', 'red', 4)
+'''
+for _, row in headend_gdf.iterrows():
+        folium.Marker(
+            location=[row.geometry.y, row.geometry.x],
+            icon=folium.Icon(color='darkblue', icon='cloud'),
+            popup="Headend"
+        ).add_to(m)
+'''
+add_fdh_markers(fdh_gdf, 'darkpurple', 'glyphicon-tower')
+
+# Save the map to an HTML file
+m.save('network_map.html')
+
+print("Map has been saved to network_map.html.")

create_mst_clusters.py

@@ -0,0 +1,89 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+from scipy.spatial import KDTree
+from shapely.geometry import Point, LineString
+from shapely.ops import nearest_points
+import numpy as np
+
+# Load home points from a shapefile
+home_points_gdf = gpd.read_file('home_points.shp')
+
+# Load network lines from a shapefile
+network_gdf = gpd.read_file('network.shp')
+
+# Filter network lines by type
+network_gdf = network_gdf[(network_gdf['type'] == 'Underground') | (network_gdf['type'] == 'Aerial')]
+
+# Extract coordinates as a NumPy array
+coords = np.array(list(zip(home_points_gdf.geometry.x, home_points_gdf.geometry.y)))
+
+# Create a KDTree for fast nearest neighbor search
+tree = KDTree(coords)
+
+# Parameters for clustering
+distance_threshold = 700  # feet
+max_homes_per_cluster = 9
+
+# Initialize clusters and a flag for each home indicating if it has been assigned to a cluster
+clusters = []
+is_clustered = np.zeros(len(home_points_gdf), dtype=bool)
+
+# Greedy clustering
+for i, coord in enumerate(coords):
+    if is_clustered[i]:
+        continue
+
+    # Find all points within the distance threshold
+    indices = tree.query_ball_point(coord, r=distance_threshold)
+    
+    # Remove points that are already clustered
+    indices = [idx for idx in indices if not is_clustered[idx]]
+    
+    # If more points than max cluster size, keep only the nearest ones
+    if len(indices) > max_homes_per_cluster:
+        distances, nearest_indices = tree.query(coord, k=max_homes_per_cluster)
+        indices = nearest_indices.tolist()
+    
+    # Create the cluster and update the flags
+    clusters.append(indices)
+    is_clustered[indices] = True
+
+# Function to snap point to nearest line
+def snap_to_nearest_line(point, lines_gdf):
+    nearest_line = None
+    min_distance = float('inf')
+
+    for line in lines_gdf.geometry:
+        if line.distance(point) < min_distance:
+            nearest_line = line
+            min_distance = line.distance(point)
+
+    nearest_point_on_line = nearest_points(point, nearest_line)[1]
+    return nearest_point_on_line
+
+# Calculate the centroids for each cluster and snap to nearest line
+cluster_centroids = []
+for cluster_indices in clusters:
+    cluster_coords = coords[cluster_indices]
+    centroid = Point(np.mean(cluster_coords[:, 0]), np.mean(cluster_coords[:, 1]))
+    snapped_point = snap_to_nearest_line(centroid, network_gdf)
+    cluster_centroids.append({'geometry': snapped_point})
+
+# Create a GeoDataFrame for cluster centroids
+clusters_gdf = gpd.GeoDataFrame(cluster_centroids, crs=home_points_gdf.crs)
+
+# Save the cluster centroids to mst.shp
+clusters_gdf.to_file('mst.shp')
+
+# Add a new attribute 'mst' to the home points GeoDataFrame
+home_points_gdf['mst'] = None
+
+# Assign the 'mst' id from the clusters to the corresponding home points
+for idx, cluster_indices in enumerate(clusters):
+    home_points_gdf.loc[cluster_indices, 'mst'] = idx
+
+# Save the home points with the 'mst' attribute to a new shapefile (home_points_clustered.shp)
+home_points_gdf.to_file('home_points_clustered.shp')
+
+print(f"{len(clusters)} clusters have been created and saved to mst.shp.")

create_mst_clusters_v2.py

@@ -0,0 +1,72 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import networkx as nx
+from shapely.geometry import Point, box
+
+# Load shapefiles
+edges_gdf = gpd.read_file('edges.shp')
+nodes_gdf = gpd.read_file('nodes.shp')
+home_points_gdf = gpd.read_file('home_points.shp')
+
+# Parameters for drop types
+drop_length = {'Buried Drop': 1000, 'Aerial Drop': 1000}
+
+# Initialize MST info, MST locations, and a structure to track homes to MST assignments
+mst_info = {}
+mst_locations = {}
+mst_id_counter = 1
+mst_to_homes = {}  # Correctly initialize mst_to_homes here
+
+# Build initial graph
+G = nx.Graph()
+for _, edge in edges_gdf.iterrows():
+    G.add_edge(edge['start_node'], edge['end_node'], weight=edge['cost'], type=edge['type'])
+
+# Spatial indexing on nodes
+nodes_gdf['geometry'] = nodes_gdf.apply(lambda row: Point(row['geometry'].x, row['geometry'].y), axis=1)
+nodes_sindex = nodes_gdf.sindex
+
+#Initialize a structure to hold MSTs and their assigned homes
+mst_assignments = {}
+
+# Function to determine the closest MST with available capacity
+def find_closest_available_mst(home_point, max_length):
+    closest_mst = None
+    min_distance = float('inf')
+    for mst_id, mst_data in mst_assignments.items():
+        if len(mst_data['homes']) < 9:
+            distance = home_point.distance(mst_data['geometry'])
+            if distance <= max_length and distance < min_distance:
+                min_distance = distance
+                closest_mst = mst_id
+    return closest_mst
+
+# Build MST assignments
+for _, home in home_points_gdf.iterrows():
+    home_point = Point(home['geometry'].x, home['geometry'].y)
+    max_length = 1000  # Adjusted allowable length
+
+    # Attempt to find an existing MST within the allowable length with capacity
+    closest_mst = find_closest_available_mst(home_point, max_length)
+    
+    if closest_mst:
+        # Assign home to the closest available MST
+        mst_assignments[closest_mst]['homes'].append(home['drop_point'])
+    else:
+        # Create a new MST for this home
+        new_mst_id = f'MST_{len(mst_assignments) + 1}'
+        mst_assignments[new_mst_id] = {'geometry': home_point, 'homes': [home['drop_point']]}
+
+# Update home_points_gdf with assigned MST IDs
+for mst_id, mst_data in mst_assignments.items():
+    for home_id in mst_data['homes']:
+        home_points_gdf.loc[home_points_gdf['drop_point'] == home_id, 'mst'] = mst_id
+
+# Prepare and save the MST locations as a GeoDataFrame
+mst_data = [{'mst_id': mst_id, 'geometry': mst_data['geometry']} for mst_id, mst_data in mst_assignments.items()]
+mst_gdf = gpd.GeoDataFrame(mst_data, geometry='geometry', crs=home_points_gdf.crs)
+mst_gdf.to_file('mstv2.shp')
+
+# Save updated home points with MST ID
+home_points_gdf.to_file('updated_home_points_with_mst.shp')

create_network.py

@@ -0,0 +1,78 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import networkx as nx
+import time
+
+print("Building the graph...")
+
+# Load shapefiles
+edges_gdf = gpd.read_file('edges.shp')
+nodes_gdf = gpd.read_file('nodes.shp')
+home_points_gdf = gpd.read_file('home_points.shp')
+
+# Build the graph
+G = nx.Graph()
+for _, edge in edges_gdf.iterrows():
+    G.add_edge(edge['start_node'], edge['end_node'], weight=edge['cost'], type=edge['type'], length=edge['length'], cost=edge['cost'])
+
+# Identify home nodes
+home_nodes = set(home_points_gdf['drop_point'].dropna())
+
+# Create a new graph that only contains the home nodes
+home_graph = nx.Graph()
+
+total = len(home_nodes)
+counter = 0
+
+target_node = '202'
+
+# Store the edges connected to home nodes
+home_node_edges = {}
+for node in home_nodes:
+    home_node_edges[node] = list(G.edges(node, data=True))
+
+# Remove all home nodes from the graph
+G_temp = G.copy()
+G_temp.remove_nodes_from(home_nodes)
+
+# Start the timer
+start_time = time.time()
+
+# Find the minimum cost path from each home node to the target node
+for start_node in home_nodes:
+    # Add the current start node back to the graph along with its edges
+    G_temp.add_node(start_node)
+    G_temp.add_edges_from(home_node_edges[start_node])
+
+    try:
+        path = nx.shortest_path(G_temp, start_node, target_node, weight='cost')
+        for i in range(len(path) - 1):
+            home_graph.add_edge(path[i], path[i+1])
+    except nx.NetworkXNoPath:
+        pass
+    finally:
+        # Remove the start node again
+        G_temp.remove_node(start_node)
+
+    counter += 1
+    print(f'Progress: {counter}/{total}', end='\r')
+
+# Stop the timer and print the elapsed time
+end_time = time.time()
+elapsed_time = end_time - start_time
+print(f"\nDone in {elapsed_time:.2f} seconds.")
+
+# Create a GeoDataFrame for the network
+network_data = []
+
+# Iterate over the edges in the home_graph and use edge geometry from edges.shp
+for edge in home_graph.edges():
+    start_node, end_node = edge
+    edge_data = edges_gdf[((edges_gdf['start_node'] == start_node) & (edges_gdf['end_node'] == end_node)) | ((edges_gdf['start_node'] == end_node) & (edges_gdf['end_node'] == start_node))]
+    if not edge_data.empty:
+        edge_geom = edge_data.iloc[0]['geometry']
+        network_data.append({'geometry': edge_geom, 'type': edge_data.iloc[0]['type'], 'length': edge_data.iloc[0]['length'], 'cost': edge_data.iloc[0]['cost']})
+
+network_gdf = gpd.GeoDataFrame(network_data, crs=edges_gdf.crs)
+network_gdf.to_file('network.shp')

create_network_v2.1.py

@@ -0,0 +1,140 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import networkx as nx
+import time
+
+# Load shapefiles
+edges_gdf = gpd.read_file('edges.shp')
+nodes_gdf = gpd.read_file('nodes.shp')
+home_points_gdf = gpd.read_file('home_points.shp')
+fdh_gdf = gpd.read_file('fdh.shp')
+
+print("Grouping homes by FDH and sorting by distance to FDH...")
+
+# Build the graph
+G = nx.Graph()
+for _, edge in edges_gdf.iterrows():
+    G.add_edge(edge['start_node'], edge['end_node'], weight=edge['cost'], type=edge['type'], length=edge['length'], cost=edge['cost'])
+
+# Create a mapping from fdh_id to node_id for quick lookup
+fdh_to_node = fdh_gdf.set_index('id')['node_id'].to_dict()
+
+# Identify home nodes
+home_nodes = set(home_points_gdf['drop_point'].dropna())
+
+# Store the edges connected to home nodes for re-adding later
+home_node_edges = {}
+for node in home_nodes:
+    home_node_edges[node] = list(G.edges(node, data=True))
+
+# Group home points by FDH
+fdh_to_homes = {}
+for _, home_point in home_points_gdf.iterrows():
+    fdh_id = home_point['fdh_id']
+    if fdh_id not in fdh_to_homes:
+        fdh_to_homes[fdh_id] = []
+    fdh_to_homes[fdh_id].append(home_point['drop_point'])
+
+# Before starting the sorting, ensure the graph is ready
+G_temp = G.copy()
+G_temp.remove_nodes_from(home_nodes)  # Initially remove all home nodes
+
+# Sort homes within each FDH group by proximity to FDH
+for fdh_id, homes in fdh_to_homes.items():
+    target_node = fdh_to_node.get(fdh_id)
+    if target_node:
+        homes_with_distance = []
+        for home in homes:
+            if home in home_node_edges:
+                # Temporarily add back the home node and its edges for distance calculation
+                G_temp.add_node(home)
+                G_temp.add_edges_from(home_node_edges[home])
+
+                try:
+                    # Calculate distance only if both nodes are present
+                    if G_temp.has_node(home) and G_temp.has_node(target_node):
+                        distance = nx.shortest_path_length(G_temp, source=home, target=target_node, weight='length')
+                        homes_with_distance.append((home, distance))
+                except nx.NetworkXNoPath:
+                    print(f"No path found from home node {home} to FDH node {target_node}.")
+                finally:
+                    # Remove the home node again to ensure it's not included in the next home's calculation
+                    G_temp.remove_node(home)
+
+        # Sort homes by calculated distance
+        homes_sorted = sorted(homes_with_distance, key=lambda x: x[1])
+        fdh_to_homes[fdh_id] = [home for home, _ in homes_sorted]  # Update with sorted homes
+
+# Create a new graph to store paths
+home_graph = nx.Graph()
+
+# Start the timer
+print("Building the network...")
+start_time = time.time()
+
+# Process each FDH group
+for fdh_id, homes in fdh_to_homes.items():
+    print(f"Processing FDH {fdh_id} with {len(homes)} homes...")
+    target_node = fdh_to_node.get(fdh_id)
+    if not target_node:
+        continue  # Skip if FDH node is not in the graph
+
+    # Temporary copy of G to modify for each pathfinding operation, excluding home nodes initially
+    G_temp = G.copy()
+    G_temp.remove_nodes_from(home_nodes)
+
+    for start_node in homes:
+        if start_node in home_node_edges:
+            # Temporarily add the start node and its edges back to G_temp for pathfinding
+            G_temp.add_node(start_node)
+            G_temp.add_edges_from(home_node_edges[start_node])
+
+            try:
+                if G_temp.has_node(target_node):  # Ensure the target FDH node exists in the graph
+                    path = nx.shortest_path(G_temp, start_node, target_node, weight='cost')
+                    for i in range(len(path) - 1):
+                        # Include fdh_id as an edge attribute
+                        home_graph.add_edge(path[i], path[i+1], weight=G[path[i]][path[i+1]]['cost'], fdh_id=fdh_id)
+            except nx.NetworkXNoPath:
+                print(f"No path found from home node {start_node} to FDH node {target_node}.")
+            finally:
+                # Remove the start node again to reset G_temp for the next iteration
+                G_temp.remove_node(start_node)
+
+# Stop the timer and print the elapsed time
+end_time = time.time()
+elapsed_time = end_time - start_time
+print(f"\nNetwork complete in {elapsed_time:.2f} seconds.")
+
+print("Saving the network...")
+
+# Extract edges data from the home_graph to create a GeoDataFrame
+network_data = []
+
+edge_counter = 0
+total_edges = len(home_graph.edges(data=True))
+
+for edge in home_graph.edges(data=True):
+    start_node, end_node, edge_attrs = edge
+    edge_data = edges_gdf[((edges_gdf['start_node'] == start_node) & (edges_gdf['end_node'] == end_node)) |
+                           ((edges_gdf['start_node'] == end_node) & (edges_gdf['end_node'] == start_node))]
+    if not edge_data.empty:
+        fdh_id = edge_attrs['fdh_id']
+        edge_geom = edge_data.iloc[0]['geometry']
+        network_data.append({
+            'geometry': edge_geom,
+            'type': edge_data.iloc[0]['type'],
+            'length': edge_data.iloc[0]['length'],
+            'cost': edge_data.iloc[0]['cost'],
+            'fdh_id': fdh_id
+        })
+
+        # Update the progress indicator
+        edge_counter += 1
+        print(f'Processing edge {edge_counter}/{total_edges}', end='\r')
+
+network_gdf = gpd.GeoDataFrame(network_data, crs=edges_gdf.crs)
+network_gdf.to_file('network.shp')
+
+print("\nDone.")

create_network_v2.2.py

@@ -0,0 +1,167 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import networkx as nx
+import time
+from shapely.geometry import box
+
+# Load shapefiles
+edges_gdf = gpd.read_file('edges.shp')
+nodes_gdf = gpd.read_file('nodes.shp')
+home_points_gdf = gpd.read_file('home_points.shp')
+fdh_gdf = gpd.read_file('fdh.shp')
+
+print("Grouping homes by FDH and sorting by distance to FDH...")
+
+# Build the graph
+G = nx.Graph()
+for _, edge in edges_gdf.iterrows():
+    G.add_edge(edge['start_node'], edge['end_node'], weight=edge['cost'], type=edge['type'], length=edge['length'], cost=edge['cost'])
+
+# Create a mapping from fdh_id to node_id for quick lookup
+fdh_to_node = fdh_gdf.set_index('id')['node_id'].to_dict()
+
+# Identify home nodes
+home_nodes = set(home_points_gdf['drop_point'].dropna())
+
+# Store the edges connected to home nodes for re-adding later
+home_node_edges = {}
+for node in home_nodes:
+    home_node_edges[node] = list(G.edges(node, data=True))
+
+# Group home points by FDH
+fdh_to_homes = {}
+for _, home_point in home_points_gdf.iterrows():
+    fdh_id = home_point['fdh_id']
+    if fdh_id not in fdh_to_homes:
+        fdh_to_homes[fdh_id] = []
+    fdh_to_homes[fdh_id].append(home_point['drop_point'])
+
+# Before starting the sorting, ensure the graph is ready
+G_temp = G.copy()
+G_temp.remove_nodes_from(home_nodes)  # Initially remove all home nodes
+
+# Sort homes within each FDH group by proximity to FDH
+for fdh_id, homes in fdh_to_homes.items():
+    target_node = fdh_to_node.get(fdh_id)
+    if target_node:
+        homes_with_distance = []
+        for home in homes:
+            if home in home_node_edges:
+                # Temporarily add back the home node and its edges for distance calculation
+                G_temp.add_node(home)
+                G_temp.add_edges_from(home_node_edges[home])
+
+                try:
+                    # Calculate distance only if both nodes are present
+                    if G_temp.has_node(home) and G_temp.has_node(target_node):
+                        distance = nx.shortest_path_length(G_temp, source=home, target=target_node, weight='length')
+                        homes_with_distance.append((home, distance))
+                except nx.NetworkXNoPath:
+                    print(f"No path found from home node {home} to FDH node {target_node}.")
+                finally:
+                    # Remove the home node again to ensure it's not included in the next home's calculation
+                    G_temp.remove_node(home)
+
+        # Sort homes by calculated distance
+        homes_sorted = sorted(homes_with_distance, key=lambda x: x[1])
+        fdh_to_homes[fdh_id] = [home for home, _ in homes_sorted]  # Update with sorted homes
+
+# Create a new graph to store paths
+home_graph = nx.Graph()
+
+# Start the timer
+print("Building the network...")
+start_time = time.time()
+
+# Process each FDH group
+# Ensure nodes_gdf is indexed by node ID for efficient lookup
+nodes_gdf = nodes_gdf.set_index('id')
+nodes_sindex = nodes_gdf.sindex
+
+for fdh_id, homes in fdh_to_homes.items():
+    print(f"Processing FDH {fdh_id} with {len(homes)} homes...")
+    target_node = fdh_to_node.get(fdh_id)
+    if not target_node:
+        continue
+
+    G_temp = G.copy()
+    G_temp.remove_nodes_from(home_nodes)  # Exclude home nodes
+
+    for start_node in homes:
+        if start_node in home_node_edges and start_node in nodes_gdf.index:
+            start_geom = nodes_gdf.loc[start_node, 'geometry']
+            
+            G_temp.add_node(start_node)
+            G_temp.add_edges_from(home_node_edges[start_node])
+
+            try:
+                path_to_fdh = nx.shortest_path(G_temp, start_node, target_node, weight='cost')
+                path_to_fdh_cost = sum(G_temp[path_to_fdh[i]][path_to_fdh[i+1]]['cost'] for i in range(len(path_to_fdh)-1))
+
+                alternate_path_cost = float('inf')
+                alternate_path = []
+
+                nearest_items = list(nodes_sindex.nearest(start_geom, 1))  # Using start_geom directly
+                if nearest_items:
+                    nearest_item_index = nearest_items[0]
+                    nearest_node_id = nodes_gdf.iloc[nearest_item_index].index  # Correctly access the index of the GeoDataFrame
+                    
+                    if nearest_node_id != start_node and nearest_node_id in G_temp:
+                        path = nx.shortest_path(G_temp, start_node, nearest_node_id, weight='cost')
+                        path_cost = sum(G_temp[path[i]][path[i+1]]['cost'] for i in range(len(path)-1))
+                        if path_cost < alternate_path_cost:
+                            alternate_path_cost = path_cost
+                            alternate_path = path
+                    else:
+                        print(f"Nearest node {nearest_node_id} is not a valid node in the graph.")
+                else:
+                    print(f"No nearest node found for home node {start_node}.")
+
+                final_path = path_to_fdh if path_to_fdh_cost <= alternate_path_cost else alternate_path
+
+                for i in range(len(final_path)-1):
+                    home_graph.add_edge(final_path[i], final_path[i+1], weight=G[final_path[i]][final_path[i+1]]['cost'], fdh_id=fdh_id)
+            
+            except nx.NetworkXNoPath:
+                print(f"No path found from home node {start_node} to FDH node {target_node}.")
+            finally:
+                G_temp.remove_node(start_node)
+
+
+# Stop the timer and print the elapsed time
+end_time = time.time()
+elapsed_time = end_time - start_time
+print(f"\nNetwork complete in {elapsed_time:.2f} seconds.")
+
+print("Saving the network...")
+
+# Extract edges data from the home_graph to create a GeoDataFrame
+network_data = []
+
+edge_counter = 0
+total_edges = len(home_graph.edges(data=True))
+
+for edge in home_graph.edges(data=True):
+    start_node, end_node, edge_attrs = edge
+    edge_data = edges_gdf[((edges_gdf['start_node'] == start_node) & (edges_gdf['end_node'] == end_node)) |
+                           ((edges_gdf['start_node'] == end_node) & (edges_gdf['end_node'] == start_node))]
+    if not edge_data.empty:
+        fdh_id = edge_attrs['fdh_id']
+        edge_geom = edge_data.iloc[0]['geometry']
+        network_data.append({
+            'geometry': edge_geom,
+            'type': edge_data.iloc[0]['type'],
+            'length': edge_data.iloc[0]['length'],
+            'cost': edge_data.iloc[0]['cost'],
+            'fdh_id': fdh_id
+        })
+
+        # Update the progress indicator
+        edge_counter += 1
+        print(f'Processing edge {edge_counter}/{total_edges}', end='\r')
+
+network_gdf = gpd.GeoDataFrame(network_data, crs=edges_gdf.crs)
+network_gdf.to_file('network.shp')
+
+print("\nDone.")

create_network_v2.py

@@ -0,0 +1,106 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import networkx as nx
+import time
+
+print("Building the graph...")
+
+# Load shapefiles
+edges_gdf = gpd.read_file('edges.shp')
+nodes_gdf = gpd.read_file('nodes.shp')
+home_points_gdf = gpd.read_file('home_points.shp')
+fdh_gdf = gpd.read_file('fdh.shp')
+
+# Build the graph
+G = nx.Graph()
+for _, edge in edges_gdf.iterrows():
+    G.add_edge(edge['start_node'], edge['end_node'], weight=edge['cost'], type=edge['type'], length=edge['length'], cost=edge['cost'])
+
+# Create a mapping from fdh_id to node_id for quick lookup
+fdh_to_node = fdh_gdf.set_index('id')['node_id'].to_dict()
+
+# Identify home nodes
+home_nodes = set(home_points_gdf['drop_point'].dropna())
+
+# Store the edges connected to home nodes for re-adding later
+home_node_edges = {}
+for node in home_nodes:
+    home_node_edges[node] = list(G.edges(node, data=True))
+
+# Create a new graph to store paths
+home_graph = nx.Graph()
+
+total = len(home_points_gdf)
+counter = 0
+
+# Start the timer
+start_time = time.time()
+
+# Temporary copy of G to modify for each pathfinding operation, excluding home nodes initially
+G_temp = G.copy()
+G_temp.remove_nodes_from(home_nodes)
+
+# Find the minimum cost path from each home node to its associated FDH node
+for _, home_point in home_points_gdf.iterrows():
+    start_node = home_point['drop_point']
+    fdh_id = home_point['fdh_id']
+    if fdh_id in fdh_to_node and start_node in home_node_edges:
+        target_node = fdh_to_node[fdh_id]  # Lookup the target_node using fdh_id
+        
+        # Temporarily add the start node and its edges back to G_temp for pathfinding
+        G_temp.add_node(start_node)
+        G_temp.add_edges_from(home_node_edges[start_node])
+        
+        try:
+            if G_temp.has_node(target_node):  # Ensure the target FDH node exists in the graph
+                path = nx.shortest_path(G_temp, start_node, target_node, weight='cost')
+                for i in range(len(path) - 1):
+                    #home_graph.add_edge(path[i], path[i+1], weight=G[path[i]][path[i+1]]['cost'])
+                    # Include fdh_id as an edge attribute
+                    home_graph.add_edge(path[i], path[i+1], weight=G[path[i]][path[i+1]]['cost'], fdh_id=fdh_id)
+        except nx.NetworkXNoPath:
+            print(f"No path found from home node {start_node} to FDH node {target_node}.")
+        finally:
+            # Remove the start node again to reset G_temp for the next iteration
+            G_temp.remove_node(start_node)
+
+    counter += 1
+    print(f'Progress: {counter}/{total}', end='\r')
+
+# Stop the timer and print the elapsed time
+end_time = time.time()
+elapsed_time = end_time - start_time
+print(f"\nGraph complete in {elapsed_time:.2f} seconds.")
+
+print("Saving the network...")
+
+# Extract edges data from the home_graph to create a GeoDataFrame
+network_data = []
+
+edge_counter = 0
+total_edges = len(home_graph.edges(data=True))
+
+for edge in home_graph.edges(data=True):
+    start_node, end_node, edge_attrs = edge
+    edge_data = edges_gdf[((edges_gdf['start_node'] == start_node) & (edges_gdf['end_node'] == end_node)) |
+                           ((edges_gdf['start_node'] == end_node) & (edges_gdf['end_node'] == start_node))]
+    if not edge_data.empty:
+        fdh_id = edge_attrs['fdh_id']
+        edge_geom = edge_data.iloc[0]['geometry']
+        network_data.append({
+            'geometry': edge_geom,
+            'type': edge_data.iloc[0]['type'],
+            'length': edge_data.iloc[0]['length'],
+            'cost': edge_data.iloc[0]['cost'],
+            'fdh_id': fdh_id
+        })
+
+        # Update the progress indicator
+        edge_counter += 1
+        print(f'Processing edge {edge_counter}/{total_edges}', end='\r')
+
+network_gdf = gpd.GeoDataFrame(network_data, crs=edges_gdf.crs)
+network_gdf.to_file('network.shp')
+
+print("\nDone.")

create_network_v3.py

@@ -0,0 +1,91 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import networkx as nx
+import time
+from shapely.geometry import LineString
+
+print("Building the graph...")
+
+# Load shapefiles
+edges_gdf = gpd.read_file('edges.shp')
+nodes_gdf = gpd.read_file('nodes.shp').set_index('id')
+home_points_gdf = gpd.read_file('home_points.shp')
+fdh_gdf = gpd.read_file('fdh.shp')  # Load fdh.shp
+
+# Build the graph
+G = nx.Graph()
+for _, edge in edges_gdf.iterrows():
+    G.add_edge(edge['start_node'], edge['end_node'], weight=edge['cost'], type=edge['type'], length=edge['length'], cost=edge['cost'])
+
+# Create a mapping from fdh_id to node_id for quick lookup
+fdh_to_node = fdh_gdf.set_index('id')['node_id'].to_dict()
+
+# Identify home nodes
+home_nodes = set(home_points_gdf['drop_point'].dropna())
+
+# Store the edges connected to home nodes for re-adding later
+home_node_edges = {}
+for node in home_nodes:
+    home_node_edges[node] = list(G.edges(node, data=True))
+
+# Create a new graph to store paths and a list for circuits
+home_graph = nx.Graph()
+circuits_data = []
+
+total = len(home_points_gdf)
+counter = 0
+
+# Start the timer
+start_time = time.time()
+
+# Temporary copy of G to modify for each pathfinding operation, excluding home nodes initially
+G_temp = G.copy()
+G_temp.remove_nodes_from(home_nodes)
+
+# Find the minimum cost path from each home node to its associated FDH node and construct circuits
+for _, home_point in home_points_gdf.iterrows():
+    start_node = home_point['drop_point']
+    fdh_id = home_point['fdh_id']
+    if fdh_id in fdh_to_node and start_node in home_node_edges:
+        target_node = fdh_to_node[fdh_id]  # Lookup the target_node using fdh_id
+        
+        # Temporarily add the start node and its edges back to G_temp for pathfinding
+        G_temp.add_node(start_node)
+        G_temp.add_edges_from(home_node_edges[start_node])
+        
+        try:
+            if G_temp.has_node(target_node):  # Ensure the target FDH node exists in the graph
+                path = nx.shortest_path(G_temp, start_node, target_node, weight='cost')
+                # Calculate total cost for the path
+                path_cost = sum(G_temp[path[i]][path[i+1]]['cost'] for i in range(len(path)-1))
+                path_geoms = [nodes_gdf.loc[node, 'geometry'] for node in path if node in nodes_gdf.index]
+                if len(path_geoms) > 1:  # Ensure there are at least two points to form a line
+                    linestring = LineString(path_geoms)
+                    circuits_data.append({
+                        'geometry': linestring,
+                        'home_node': start_node,
+                        'fdh_id': fdh_id,
+                        'cost': path_cost  # Include the total path cost
+                    })
+                for i in range(len(path) - 1):
+                    home_graph.add_edge(path[i], path[i+1], weight=G[path[i]][path[i+1]]['cost'], fdh_id=fdh_id)
+        except nx.NetworkXNoPath:
+            print(f"No path found from home node {start_node} to FDH node {target_node}.")
+        finally:
+            # Remove the start node again to reset G_temp for the next iteration
+            G_temp.remove_node(start_node)
+
+    counter += 1
+    print(f'Progress: {counter}/{total}', end='\r')
+
+# Stop the timer and print the elapsed time
+end_time = time.time()
+elapsed_time = end_time - start_time
+print(f"\nGraph complete in {elapsed_time:.2f} seconds.")
+
+# Create GeoDataFrame for circuits and save as shapefile
+circuits_gdf = gpd.GeoDataFrame(circuits_data, crs=edges_gdf.crs)
+circuits_gdf.to_file('circuits.shp')
+
+print("\nDone.")

create_network_v4.py

@@ -0,0 +1,160 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import networkx as nx
+import numpy as np
+from scipy.spatial import KDTree, distance
+import pandas as pd
+import time
+
+print("Building the graph...")
+
+# Load shapefiles
+edges_gdf = gpd.read_file('edges.shp')
+nodes_gdf = gpd.read_file('nodes.shp')
+home_points_gdf = gpd.read_file('home_points.shp')
+fdh_gdf = gpd.read_file('fdh.shp')
+
+# Build the graph
+G = nx.Graph()
+for _, edge in edges_gdf.iterrows():
+    G.add_edge(edge['start_node'], edge['end_node'], weight=edge['cost'], type=edge['type'], length=edge['length'])
+
+# Create KDTree for efficient spatial queries (for nearest node lookups)
+node_coords = np.array(list(zip(nodes_gdf.geometry.x, nodes_gdf.geometry.y)))
+tree = KDTree(node_coords)
+node_id_to_idx = {node_id: idx for idx, node_id in enumerate(nodes_gdf['id'])}
+
+# Function to find nearest node in the graph to a given point
+def find_nearest_node(point):
+    _, idx = tree.query(point, k=1)
+    return nodes_gdf.iloc[idx]['id']
+
+# Group homes by FDH
+fdh_to_homes = home_points_gdf.groupby('fdh_id')['drop_point'].apply(list).to_dict()
+
+# KDTree for FDH locations to sort homes by distance to FDH
+fdh_coords = np.array(list(zip(fdh_gdf.geometry.x, fdh_gdf.geometry.y)))
+fdh_tree = KDTree(fdh_coords)
+fdh_id_to_idx = {fdh_id: idx for idx, fdh_id in enumerate(fdh_gdf['id'])}
+
+# Assuming fdh_gdf has a 'geometry' column and each FDH can be mapped to the nearest network node
+fdh_to_node = {}
+for _, fdh_row in fdh_gdf.iterrows():
+    fdh_point = fdh_row.geometry
+    nearest_node_id = find_nearest_node((fdh_point.x, fdh_point.y))
+    fdh_to_node[fdh_row['id']] = nearest_node_id
+
+# Function to find the nearest node in the optimized graph to a given start node
+def find_nearest_node_in_optimized_graph(optimized_graph, G, start_node, nodes_gdf):
+    if optimized_graph.number_of_nodes() == 0:
+        # If the optimized graph has no nodes, return None values
+        return None, [], float('inf')
+    
+    # Extract coordinates for the start node and all nodes in the optimized graph
+    start_coords = nodes_gdf.loc[nodes_gdf['id'] == start_node, ['geometry']].values[0][0]
+    optimized_nodes_coords = {node: nodes_gdf.loc[nodes_gdf['id'] == node, ['geometry']].values[0][0] for node in optimized_graph.nodes}
+
+    # Calculate distances from start node to each node in the optimized graph
+    distances = {node: distance.euclidean((start_coords.x, start_coords.y), (coords.x, coords.y)) for node, coords in optimized_nodes_coords.items()}
+    
+    # Find the nearest node and its distance
+    nearest_node = min(distances, key=distances.get)
+    nearest_distance = distances[nearest_node]
+    
+    # Calculate the shortest path from start node to nearest node in the original graph G
+    try:
+        path = nx.shortest_path(G, source=start_node, target=nearest_node, weight='weight')
+        path_cost = sum(G[u][v]['weight'] for u, v in zip(path[:-1], path[1:]))
+        return nearest_node, path, path_cost
+    except nx.NetworkXNoPath:
+        # If no path exists, return None values
+        return None, [], float('inf')
+
+# Function to sort homes by distance to FDH
+def sort_homes_by_distance(fdh_id, homes):
+    fdh_idx = fdh_id_to_idx[fdh_id]
+    fdh_point = fdh_coords[fdh_idx]
+    home_points = np.array([nodes_gdf.loc[nodes_gdf['id'] == h].geometry.apply(lambda x: (x.x, x.y)).values[0] for h in homes])
+    distances = np.linalg.norm(home_points - fdh_point, axis=1)
+    sorted_idx = np.argsort(distances)
+    return [homes[i] for i in sorted_idx]
+
+# Start the optimization process
+start_time = time.time()
+
+optimized_graph = nx.Graph()
+
+path_cache = {}
+
+def calculate_direct_path_cost(G, start_node, target_node):
+    cache_key = (start_node, target_node)
+    if cache_key in path_cache:
+        return path_cache[cache_key]
+
+    try:
+        path = nx.shortest_path(G, source=start_node, target=target_node, weight='weight')
+        cost = sum(G[u][v]['weight'] for u, v in zip(path[:-1], path[1:]))
+        path_cache[cache_key] = (cost, path)
+        return cost, path
+    except nx.NetworkXNoPath:
+        path_cache[cache_key] = (float('inf'), [])
+        return float('inf'), []
+
+# Placeholder for a function that updates the optimized graph with a new path
+def update_network_with_path(optimized_graph, path, G):
+    for start, end in zip(path[:-1], path[1:]):
+        if not optimized_graph.has_edge(start, end):
+            optimized_graph.add_edge(start, end, weight=G[start][end]['weight'])
+
+# Before starting the main loop, initialize variables to track progress
+total_homes = sum(len(homes) for homes in fdh_to_homes.values())
+processed_homes = 0
+
+# Main loop for pathfinding and updating the optimized graph
+for fdh_id, homes in fdh_to_homes.items():
+    target_node = fdh_to_node[fdh_id]  # Assuming fdh_to_node is correctly defined earlier
+    sorted_homes = sort_homes_by_distance(fdh_id, homes)
+    for home in sorted_homes:
+        start_node = home  # Assuming 'home' is already a node ID in G
+
+        # Calculate direct path cost from home to FDH
+        direct_cost, direct_path = calculate_direct_path_cost(G, start_node, target_node)
+
+        # Find the nearest node in the optimized graph to the home
+        nearest_optimized_node, nearest_optimized_path, nearest_optimized_cost = find_nearest_node_in_optimized_graph(optimized_graph, G, start_node, nodes_gdf)
+
+        # Calculate path cost from the nearest optimized node to FDH through the existing network
+        if nearest_optimized_node is not None:
+            _, path_from_nearest = calculate_direct_path_cost(optimized_graph, nearest_optimized_node, target_node)
+            total_optimized_cost = nearest_optimized_cost + sum(optimized_graph[u][v]['weight'] for u, v in zip(path_from_nearest[:-1], path_from_nearest[1:]))
+        else:
+            total_optimized_cost = float('inf')
+
+        # Compare and update the optimized graph with the cheaper path
+        if direct_cost <= total_optimized_cost:
+            update_network_with_path(optimized_graph, direct_path, G)
+        else:
+            # Update with path to nearest node then to FDH
+            update_network_with_path(optimized_graph, nearest_optimized_path + path_from_nearest[1:], G)
+
+        # Update and print the progress
+        processed_homes += 1
+        progress_percentage = (processed_homes / total_homes) * 100
+        print(f"\rProgress: {processed_homes}/{total_homes} homes processed ({progress_percentage:.2f}%)", end='')
+
+# Print a newline character at the end to ensure any following output starts on a new line
+print("\nOptimization complete.")
+
+# End the optimization process
+end_time = time.time()
+elapsed_time = end_time - start_time
+print(f"Optimization complete in {elapsed_time:.2f} seconds.")
+
+# TODO: Add the implementation for direct and existing network path calculations and updates
+
+# Saving the optimized network (this is a placeholder, adapt to your actual data structure)
+# network_gdf = gpd.GeoDataFrame([...], crs=edges_gdf.crs)
+# network_gdf.to_file('optimized_network.shp')
+
+print("Optimization done. Network saved.")

create_nodes.py

@@ -0,0 +1,61 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import pandas as pd
+import itertools
+from shapely.geometry import Point
+
+print("Creating nodes...")
+
+# Load the EDGES shapefile
+edges_gdf = gpd.read_file('edges.shp')
+
+# Prepare a GeoDataFrame for the NODES
+nodes_gdf = gpd.GeoDataFrame(columns=['id', 'geometry'], crs=edges_gdf.crs)
+
+# Prepare fields for start_node and end_node in the EDGES GeoDataFrame
+edges_gdf['start_node'] = None
+edges_gdf['end_node'] = None
+
+# To store the unique nodes
+unique_nodes = {}
+
+# Using itertools.count to generate unique IDs
+id_counter = itertools.count(start=1)
+
+# Define a tolerance threshold for coordinate comparison
+epsilon = 2
+
+# Function to handle node creation and ID assignment
+def handle_node(point, unique_nodes, nodes_gdf):
+    for existing_point, node_id in unique_nodes.items():
+        if abs(point[0] - existing_point[0]) < epsilon and abs(point[1] - existing_point[1]) < epsilon:
+            return unique_nodes, nodes_gdf, node_id
+    
+    node_id = next(id_counter)
+    unique_nodes[point] = node_id
+    node_df = gpd.GeoDataFrame({'id': [node_id], 'geometry': [Point(point)]}, crs=edges_gdf.crs)
+    nodes_gdf = pd.concat([nodes_gdf, node_df], ignore_index=True)
+    return unique_nodes, nodes_gdf, node_id
+
+for index, edge in edges_gdf.iterrows():
+    start_point, end_point = edge.geometry.coords[0], edge.geometry.coords[-1]
+    
+    # Handle start node
+    unique_nodes, nodes_gdf, start_node_id = handle_node(start_point, unique_nodes, nodes_gdf)
+
+    # Handle end node
+    unique_nodes, nodes_gdf, end_node_id = handle_node(end_point, unique_nodes, nodes_gdf)
+
+    # Populate the start_node and end_node fields for the EDGES GeoDataFrame
+    edges_gdf.at[index, 'start_node'] = start_node_id
+    edges_gdf.at[index, 'end_node'] = end_node_id
+
+# Save the NODES GeoDataFrame to a new shapefile
+nodes_gdf.to_file('nodes.shp')
+print("Nodes have been saved to nodes.shp.")
+
+# Optionally, save the updated EDGES GeoDataFrame with start_node and end_node attributes
+edges_gdf.to_file('edges.shp')
+print("Edges with node attributes have been saved to edges.shp.")
+print("Done.")

create_transitions.py

@@ -0,0 +1,135 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+from shapely.geometry import LineString, Point, GeometryCollection
+from shapely.ops import nearest_points, split
+import pandas as pd
+
+# Constants
+SEARCH_RADIUS = 50  # feet, assuming your spatial data is in a feet-based coordinate system
+COST_PER_FOOT = 1000  # cost per foot
+
+# Load the shapefiles
+poles_gdf = gpd.read_file('poles.shp')
+edges_gdf = gpd.read_file('edges.shp')
+
+# Separating 'Underground' edges from other edge types
+underground_gdf = edges_gdf[edges_gdf['type'] == 'Underground']
+other_edges_gdf = edges_gdf[edges_gdf['type'] != 'Underground']
+
+# Creating spatial index for underground_gdf
+underground_gdf_sindex = underground_gdf.sindex
+
+# Set to keep track of indices of split edges
+split_edges_indices = set()
+
+# Function to split a line at a given point and return the segments
+def split_line_at_point(line, point):
+    """Split a line at a given point and return the segments."""
+    # Create a temporary tiny line at the point to use for splitting
+    buffer = Point(point).buffer(0.0001)  # Adjust buffer size if needed for precision
+    split_line = split(line, buffer)
+
+    # Check if the result is a GeometryCollection
+    if isinstance(split_line, GeometryCollection):
+        # Iterate through the geometries in the GeometryCollection
+        return [geom for geom in split_line.geoms if isinstance(geom, LineString)]
+    else:
+        return [split_line]
+    
+def find_closest_edge(pole, underground_gdf, underground_gdf_sindex):
+    # Using spatial index to find nearby edges
+    possible_matches_index = list(underground_gdf_sindex.intersection(pole.geometry.buffer(SEARCH_RADIUS).bounds))
+    possible_matches = underground_gdf.iloc[possible_matches_index]
+
+    closest_edge = None
+    closest_point = None
+    min_dist = float('inf')
+
+    # Iterate through potential nearby edges
+    for edge_idx, edge in possible_matches.iterrows():
+        # Find the closest point on the current edge to the pole
+        point_on_edge = nearest_points(pole.geometry, edge.geometry)[1]
+        dist = pole.geometry.distance(point_on_edge)
+
+        # Update the closest edge if this one is closer
+        if dist < min_dist:
+            min_dist = dist
+            closest_edge = edge
+            closest_point = point_on_edge
+
+    return closest_edge, closest_point, min_dist
+
+# List to store new edges and transitions
+new_edges = []
+new_transitions = []
+
+# Before the loop, filter out poles without valid geometries
+valid_poles_gdf = poles_gdf[poles_gdf.geometry.notnull()]
+
+# Progress indicator setup
+total = len(valid_poles_gdf)
+counter = 0
+
+# Iterate over every pole
+for idx, pole in valid_poles_gdf.iterrows():
+    counter += 1
+    print(f'Processing pole {counter}/{total}', end='\r')
+
+    # Update the spatial index for the current state of underground_gdf
+    underground_gdf_sindex = underground_gdf.sindex
+
+    # Find the closest edge to the current pole
+    closest_edge, closest_point, min_dist = find_closest_edge(pole, underground_gdf, underground_gdf_sindex)
+
+    # Check if the closest edge is within the search radius
+    if closest_edge is not None and min_dist <= SEARCH_RADIUS:
+        # Create a transition edge to the closest point
+        #transition_line = LineString([pole.geometry, closest_point])
+        # Extract coordinates from Point objects before creating the LineString
+        transition_line = LineString([(pole.geometry.x, pole.geometry.y), (closest_point.x, closest_point.y)])
+        transition_length = transition_line.length
+        transition_cost = transition_length * COST_PER_FOOT
+        transition = {
+            'type': 'Transition',
+            'length': transition_length,
+            'cost': transition_cost,
+            'geometry': transition_line
+        }
+        new_transitions.append(transition)
+
+        # Split the closest underground edge
+        if closest_point.coords[:] not in closest_edge.geometry.coords[:]:
+            split_segments = split_line_at_point(closest_edge.geometry, closest_point.coords[:])
+            
+            # Remove the split edge from underground_gdf
+            underground_gdf = underground_gdf.drop(closest_edge.name)
+
+            # Add new split segments to underground_gdf
+            for segment in split_segments:
+                new_edge = {
+                    'type': 'Underground',
+                    'length': segment.length,
+                    'cost': segment.length * COST_PER_FOOT,
+                    'geometry': segment
+                }
+                # Create a GeoDataFrame for the new edge and append it
+                new_edge_gdf = gpd.GeoDataFrame([new_edge], crs=poles_gdf.crs)
+                underground_gdf = pd.concat([underground_gdf, new_edge_gdf], ignore_index=True)
+
+            # Update spatial index after modification
+            underground_gdf_sindex = underground_gdf.sindex
+
+# Filter out the split 'Underground' edges
+unsplit_underground_gdf = underground_gdf[~underground_gdf.index.isin(split_edges_indices)]
+
+# Create GeoDataFrames for the new edges and transitions
+transitions_gdf = gpd.GeoDataFrame(new_transitions, crs=poles_gdf.crs)
+
+# Concatenate the unsplit underground edges, the new (updated) underground edges, transitions, and other edge types
+combined_edges_gdf = pd.concat([other_edges_gdf, underground_gdf, transitions_gdf], ignore_index=True)
+
+# Save the updated edges GeoDataFrame back to the shapefile
+combined_edges_gdf.to_file('edges.shp')
+
+print("\nProcessing complete.")

design_from_roads.sh

@@ -0,0 +1,10 @@
+#!/bin/bash
+
+drops_split_centerlines.py
+create_nodes.py
+associate_drop_points_v2.py
+cluster_fdh_v2.py
+create_network_v2.py
+create_mst_clusters.py
+report.py
+create_map_buried.py

design_from_start.sh

@@ -0,0 +1,15 @@
+#!/bin/bash
+
+centerlines_from_homes.py
+drops_split_centerlines.py
+create_aerial_drops.py
+create_aerial_edges.py
+create_transitions.py
+create_nodes.py
+associate_drop_points_v2.py
+cluster_fdh_v2.py
+create_network_v2.py
+create_mst_clusters.py
+poles_used.py
+report.py
+create_map.py

design_underground_only.sh

@@ -0,0 +1,11 @@
+#!/bin/bash
+
+centerlines_from_homes.py
+drops_split_centerlines.py
+create_nodes.py
+associate_drop_points_v2.py
+cluster_fdh_v2.py
+create_network_v2.py
+create_mst_clusters.py
+report.py
+create_map_buried.py

drops_split_centerlines.py

@@ -0,0 +1,92 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+from shapely.geometry import LineString, Point
+from shapely.affinity import translate
+from shapely.ops import split, unary_union, nearest_points
+from geopy.distance import great_circle
+import math
+import pandas as pd
+
+# Define constants
+underground_cpf = 1000.00
+buried_drop_cpf = 200.00
+
+# Function definitions
+def extend_line(line, extension_length):
+    """Extend a line by a given length."""
+    if isinstance(line, LineString) and len(line.coords) >= 2:
+        start_point = line.coords[0]
+        end_point = line.coords[-1]
+        dx = end_point[0] - start_point[0]
+        dy = end_point[1] - start_point[1]
+        length = math.sqrt(dx**2 + dy**2)
+        extension_dx = dx / length * extension_length
+        extension_dy = dy / length * extension_length
+        extended_line = LineString([start_point, translate(Point(end_point), extension_dx, extension_dy).coords[0]])
+        return extended_line
+    return line
+
+def closest_point_on_line(point, lines):
+    """Find the closest point on the closest line to the given point."""
+    nearest_line = min(lines, key=lambda line: line.distance(point))
+    return nearest_points(point, nearest_line)[1]
+
+# Load data from shapefiles
+gdf_homes = gpd.read_file('home_points.shp')
+gdf_centerlines = gpd.read_file('road_centerlines.shp')
+
+# Ensure CRS match
+gdf_homes = gdf_homes.to_crs(gdf_centerlines.crs)
+
+print("Processing drops...")
+
+# Process each home point
+drops = []
+drops_ext = []
+for idx, home in gdf_homes.iterrows():
+    nearest_point = closest_point_on_line(home.geometry, gdf_centerlines.geometry)
+    line = LineString([home.geometry, nearest_point])
+    drops.append({'geometry': line})
+    line_ext = extend_line(line, line.length * 0.05)
+    drops_ext.append({'geometry': line_ext})
+
+# Convert drops to GeoDataFrames and save to shapefiles
+gdf_drops = gpd.GeoDataFrame(drops, crs=gdf_centerlines.crs)
+#gdf_drops.to_file('drops.shp')
+
+gdf_drops_ext = gpd.GeoDataFrame(drops_ext, crs=gdf_centerlines.crs)
+#gdf_drops_ext.to_file('drops_ext.shp')
+
+print("Splitting centerlines...")
+
+# Split the centerlines with extended drops
+split_lines = []
+for line in gdf_centerlines.geometry:
+    split_line = split(line, unary_union(gdf_drops_ext.geometry))
+    for geom in split_line.geoms:  # Iterate through the geometries in the GeometryCollection
+        split_lines.append(geom)
+
+gdf_split_roads = gpd.GeoDataFrame(geometry=split_lines, crs=gdf_centerlines.crs)
+gdf_split_roads = gdf_split_roads.drop_duplicates(subset=['geometry'])
+
+# Calculate additional attributes (length, cost, type)
+gdf_split_roads['type'] = 'Underground'
+gdf_split_roads['length'] = gdf_split_roads.geometry.length  # Length in CRS units
+gdf_split_roads['cost'] = gdf_split_roads['length'] * underground_cpf
+
+# Save outputs to shapefiles
+#gdf_edges.to_file('split_roads.shp')
+
+# Prepare drops for concatenation with edges
+gdf_drops['type'] = 'Buried Drop'
+gdf_drops['length'] = gdf_drops.geometry.length  # Length in CRS units
+gdf_drops['cost'] = gdf_drops['length'] * buried_drop_cpf
+
+# Combine drops with edges
+gdf_combined = pd.concat([gdf_split_roads, gdf_drops], ignore_index=True)
+
+# Save the combined GeoDataFrame to a shapefile
+gdf_combined.to_file('edges.shp')
+
+print("Processing complete.")

drops_split_centerlines_v2.py

@@ -0,0 +1,95 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+from shapely.geometry import LineString, Point
+from shapely.affinity import translate
+from shapely.ops import split, unary_union, nearest_points
+from geopy.distance import great_circle
+import math
+import pandas as pd
+
+# Define constants
+underground_cpf = 1000.00
+buried_drop_cpf = 200.00
+
+# Function definitions
+def extend_line(line, extension_length):
+    """Extend a line by a given length."""
+    if isinstance(line, LineString) and len(line.coords) >= 2:
+        start_point = line.coords[0]
+        end_point = line.coords[-1]
+        dx = end_point[0] - start_point[0]
+        dy = end_point[1] - start_point[1]
+        length = math.sqrt(dx**2 + dy**2)
+        extension_dx = dx / length * extension_length
+        extension_dy = dy / length * extension_length
+        extended_line = LineString([start_point, translate(Point(end_point), extension_dx, extension_dy).coords[0]])
+        return extended_line
+    return line
+
+def closest_point_on_line(point, lines):
+    """Find the closest point on the closest line to the given point."""
+    nearest_line = min(lines, key=lambda line: line.distance(point))
+    return nearest_points(point, nearest_line)[1]
+
+# Load data from shapefiles
+gdf_homes = gpd.read_file('home_points.shp')
+gdf_centerlines = gpd.read_file('road_centerlines.shp')
+
+# Ensure CRS match
+gdf_homes = gdf_homes.to_crs(gdf_centerlines.crs)
+
+print("Processing drops...")
+
+# Process each home point
+drops = []
+drops_ext = []
+for idx, home in gdf_homes.iterrows():
+    # Convert home.geometry to a 2D point by ignoring the Z coordinate
+    home_point_2d = Point(home.geometry.x, home.geometry.y)
+    nearest_point = closest_point_on_line(home.geometry, gdf_centerlines.geometry)
+    # Now, both points are 2D, and you can create a LineString without dimensionality issues
+    line = LineString([home_point_2d, nearest_point])
+    drops.append({'geometry': line})
+    line_ext = extend_line(line, line.length * 0.05)
+    drops_ext.append({'geometry': line_ext})
+
+# Convert drops to GeoDataFrames and save to shapefiles
+gdf_drops = gpd.GeoDataFrame(drops, crs=gdf_centerlines.crs)
+#gdf_drops.to_file('drops.shp')
+
+gdf_drops_ext = gpd.GeoDataFrame(drops_ext, crs=gdf_centerlines.crs)
+#gdf_drops_ext.to_file('drops_ext.shp')
+
+print("Splitting centerlines...")
+
+# Split the centerlines with extended drops
+split_lines = []
+for line in gdf_centerlines.geometry:
+    split_line = split(line, unary_union(gdf_drops_ext.geometry))
+    for geom in split_line.geoms:  # Iterate through the geometries in the GeometryCollection
+        split_lines.append(geom)
+
+gdf_split_roads = gpd.GeoDataFrame(geometry=split_lines, crs=gdf_centerlines.crs)
+gdf_split_roads = gdf_split_roads.drop_duplicates(subset=['geometry'])
+
+# Calculate additional attributes (length, cost, type)
+gdf_split_roads['type'] = 'Underground'
+gdf_split_roads['length'] = gdf_split_roads.geometry.length  # Length in CRS units
+gdf_split_roads['cost'] = gdf_split_roads['length'] * underground_cpf
+
+# Save outputs to shapefiles
+#gdf_edges.to_file('split_roads.shp')
+
+# Prepare drops for concatenation with edges
+gdf_drops['type'] = 'Buried Drop'
+gdf_drops['length'] = gdf_drops.geometry.length  # Length in CRS units
+gdf_drops['cost'] = gdf_drops['length'] * buried_drop_cpf
+
+# Combine drops with edges
+gdf_combined = pd.concat([gdf_split_roads, gdf_drops], ignore_index=True)
+
+# Save the combined GeoDataFrame to a shapefile
+gdf_combined.to_file('edges.shp')
+
+print("Processing complete.")

exchange_boundary.py

@@ -0,0 +1,74 @@
+#!/usr/bin/env python3
+
+import sys
+import geopandas as gpd
+import numpy as np
+from shapely.geometry import MultiPoint, MultiPolygon, Polygon
+from scipy.spatial import Delaunay
+
+def alpha_shape(points, alpha):
+    if len(points) < 4:
+        return MultiPoint(list(points)).convex_hull
+
+    def add_edge(edges, edge_points, coords, i, j):
+        """Add a line between the i-th and j-th points, if not in the list already"""
+        if (i, j) in edges or (j, i) in edges:
+            # already added
+            return
+        edges.add((i, j))
+        edge_points.append(coords[i])  # First point of the edge
+        edge_points.append(coords[j])  # Second point of the edge
+
+    coords = np.array([point.coords[0] for point in points])
+    tri = Delaunay(coords)
+    edges = set()
+    edge_points = []
+    for ia, ib, ic in tri.simplices:
+        pa = coords[ia]
+        pb = coords[ib]
+        pc = coords[ic]
+        a = np.linalg.norm(pa - pb)
+        b = np.linalg.norm(pb - pc)
+        c = np.linalg.norm(pc - pa)
+        s = (a + b + c) / 2.0
+        area = np.sqrt(s * (s - a) * (s - b) * (s - c))
+        circum_r = a * b * c / (4.0 * area)
+        if circum_r < 1.0 / alpha:
+            add_edge(edges, edge_points, coords, ia, ib)
+            add_edge(edges, edge_points, coords, ib, ic)
+            add_edge(edges, edge_points, coords, ic, ia)
+
+    # Make sure edge_points are tuples, as Shapely expects
+    edge_points = [tuple(pt) for pt in edge_points]
+    
+    # Create a MultiPoint object from edge_points
+    m = MultiPoint(edge_points).convex_hull
+    if m.geom_type in ['Polygon', 'MultiPolygon']:
+        return m
+    elif m.geom_type == 'GeometryCollection':
+        # Correct iteration over geometries in GeometryCollection
+        polygons = [geom for geom in m.geoms if geom.geom_type in ['Polygon', 'MultiPolygon']]
+        if polygons:
+            return MultiPolygon(polygons)
+        else:
+            # Fallback to convex hull if no polygonal geometries are found
+            print("Warning: No polygonal geometries found in the geometry collection. Falling back to convex hull.")
+            return MultiPoint(list(points)).convex_hull
+    else:
+        raise ValueError(f"Unsupported geometry type: {m.geom_type}")
+
+def create_concave_hull_shapefile(input_shapefile, output_shapefile, alpha=0.0000003):
+    points_df = gpd.read_file(input_shapefile)
+    concave_hull_polygon = alpha_shape(points_df.geometry, alpha)
+    polygon_gdf = gpd.GeoDataFrame([1], geometry=[concave_hull_polygon], columns=['dummy'])
+    polygon_gdf.crs = points_df.crs
+    polygon_gdf.to_file(output_shapefile)
+
+if __name__ == "__main__":
+    if len(sys.argv) != 2:
+        print("Usage: python script.py <path/to/your/point_shapefile.shp>")
+        sys.exit(1)
+
+    input_shapefile = sys.argv[1]
+    output_shapefile = 'exchange_boundary.shp'
+    create_concave_hull_shapefile(input_shapefile, output_shapefile)

open_all_kml.py

@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+
+import os
+import subprocess
+
+def open_with_google_earth(directory):
+    # Check if the directory exists
+    if not os.path.isdir(directory):
+        print(f"Directory {directory} does not exist.")
+        return
+
+    # Walk through all files and folders within the directory
+    for root, dirs, files in os.walk(directory):
+        for file in files:
+            if file.endswith('.kml') or file.endswith('.kmz'):
+                # Construct the full file path
+                file_path = os.path.join(root, file)
+                print(f"Opening {file_path} with Google Earth...")
+                try:
+                    # Attempt to open the file with the default application
+                    subprocess.run(['open', file_path], check=True)
+                except subprocess.CalledProcessError as e:
+                    print(f"Failed to open {file_path}: {e}")
+
+# Use the current directory as the starting point
+directory_path = os.getcwd()
+open_with_google_earth(directory_path)

poles_used.py

@@ -0,0 +1,40 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+
+# Load the shapefiles
+network_gdf = gpd.read_file('network.shp')
+poles_gdf = gpd.read_file('poles.shp')
+
+# Check for invalid geometries
+invalid_geometries = poles_gdf[~poles_gdf.is_valid]
+if not invalid_geometries.empty:
+    print(f"Found {len(invalid_geometries)} invalid geometries. Removing them.")
+    poles_gdf = poles_gdf[poles_gdf.is_valid]
+
+# Buffer value in feet (adjust as needed for your data's accuracy)
+buffer_distance = 1  # 1 foot buffer
+
+# List to collect poles that are used
+poles_used = []
+
+# Check each pole for intersection with any line in the network
+for _, pole in poles_gdf.iterrows():
+    # Apply buffer to the pole point for intersection check
+    buffered_pole = pole.geometry.buffer(buffer_distance)
+    
+    # Check if the buffered pole intersects with any line in the network
+    if any(buffered_pole.intersects(line) for line in network_gdf.geometry):
+        # Add pole to the list
+        poles_used.append(pole)
+
+# Create a GeoDataFrame from the list of used poles
+poles_used_gdf = gpd.GeoDataFrame(poles_used, columns=poles_gdf.columns, crs=poles_gdf.crs)
+
+# Save the resulting GeoDataFrame to a new shapefile
+poles_used_gdf.to_file('poles_used.shp')
+
+# Output the number of poles saved
+num_poles_saved = len(poles_used_gdf)
+print(f"{num_poles_saved} poles used have been saved to poles_used.shp.")
+

qc/drop_endpoints.py

@@ -0,0 +1,32 @@
+import geopandas as gpd
+import matplotlib.pyplot as plt
+from shapely.geometry import MultiPoint
+
+def extract_endpoints(geometry):
+    """
+    Extracts the start and end points of a LineString or MultiLineString.
+    """
+    if geometry.type == 'LineString':
+        return MultiPoint([geometry.coords[0], geometry.coords[-1]])
+    elif geometry.type == 'MultiLineString':
+        points = []
+        for line in geometry:
+            points.append(line.coords[0])
+            points.append(line.coords[-1])
+        return MultiPoint(points)
+
+# Load the shapefile
+drop_cable_gdf = gpd.read_file('drop_cable.shp')
+
+# Extract endpoints
+drop_cable_gdf['endpoints'] = drop_cable_gdf['geometry'].apply(extract_endpoints)
+
+# Create a GeoDataFrame for endpoints
+endpoints_gdf = gpd.GeoDataFrame(geometry=drop_cable_gdf['endpoints'].explode())
+
+# Plotting
+fig, ax = plt.subplots()
+drop_cable_gdf.plot(ax=ax, color='blue', label='Drop Cables')
+endpoints_gdf.plot(ax=ax, color='red', marker='o', label='Endpoints', markersize=5)
+plt.legend()
+plt.show()

qc/drop_endpoints_folium.py

@@ -0,0 +1,49 @@
+import geopandas as gpd
+import folium
+from shapely.geometry import Point
+
+def extract_endpoints(geometry):
+    """
+    Extracts the start and end points of a LineString or MultiLineString.
+    """
+    points = []
+    if geometry.type == 'LineString':
+        points.append(Point(geometry.coords[0]))
+        points.append(Point(geometry.coords[-1]))
+    elif geometry.type == 'MultiLineString':
+        for line in geometry:
+            points.append(Point(line.coords[0]))
+            points.append(Point(line.coords[-1]))
+    return points
+
+# Load the shapefile
+drop_cable_gdf = gpd.read_file('drop_cable.shp')
+
+# Check and transform CRS to EPSG:4326 if needed
+if drop_cable_gdf.crs is not None and drop_cable_gdf.crs.to_epsg() != 4326:
+    drop_cable_gdf = drop_cable_gdf.to_crs(epsg=4326)
+
+# Extract endpoints and create a list of points
+endpoint_list = []
+for geom in drop_cable_gdf.geometry:
+    endpoint_list.extend(extract_endpoints(geom))
+
+# Create a GeoDataFrame for endpoints
+endpoints_gdf = gpd.GeoDataFrame(geometry=endpoint_list, crs="EPSG:4326")
+
+# Convert GeoDataFrames to GeoJSON
+drop_cable_geojson = drop_cable_gdf.to_json()
+endpoints_geojson = endpoints_gdf.to_json()
+
+# Create a Folium map
+m = folium.Map(location=[drop_cable_gdf.geometry.iloc[0].centroid.y, drop_cable_gdf.geometry.iloc[0].centroid.x], zoom_start=12)
+
+# Add GeoJSON layers to the map
+folium.GeoJson(drop_cable_geojson, name="Drop Cables", style_function=lambda x: {'color': 'blue', 'weight': 2}).add_to(m)
+folium.GeoJson(endpoints_geojson, name="Endpoints", style_function=lambda x: {'color': 'red', 'fillColor': 'red'}).add_to(m)
+
+# Add layer control to toggle layers
+folium.LayerControl().add_to(m)
+
+# Save the map to an HTML file
+m.save('map.html')

qgis_scripts/associate_drop_points.py

@@ -0,0 +1,34 @@
+from qgis.core import (QgsFeature, QgsField, QgsVectorLayer, QgsProject, QgsSpatialIndex)
+
+print("Associating home points to drop points...")
+
+# Load the NODES and HOME_POINT layers
+nodes_layer = QgsProject.instance().mapLayersByName('NODES')[0]
+home_points_layer = QgsProject.instance().mapLayersByName('HOME_POINTS')[0]
+
+# Prepare the drop_point_id field for the HOME_POINT layer
+home_points_layer.startEditing()
+home_points_prov = home_points_layer.dataProvider()
+home_points_prov.addAttributes([QgsField("drop_point_id", QVariant.Int)])
+home_points_layer.updateFields()
+
+# Build a spatial index for the NODES layer
+index = QgsSpatialIndex()
+node_features = {f.id(): f for f in nodes_layer.getFeatures()}
+index.addFeatures(node_features.values())
+
+# Find the nearest node for each home point and update the drop_point_id attribute
+for hp in home_points_layer.getFeatures():
+    nearest_nodes = index.nearestNeighbor(hp.geometry().asPoint(), 1)
+    if nearest_nodes:
+        nearest_node_id = nearest_nodes[0]
+        nearest_node = node_features[nearest_node_id]
+        node_id = nearest_node["id"]
+        home_points_layer.changeAttributeValue(hp.id(), 
+                                               home_points_layer.fields().indexOf("drop_point_id"), 
+                                               node_id)
+
+# Stop editing the HOME_POINT layer and save changes
+home_points_layer.commitChanges()
+
+print("Done.")

qgis_scripts/centerlines_from_homes.py

@@ -0,0 +1,77 @@
+import osmnx as ox
+import geopandas as gpd
+from shapely import wkt
+from shapely.geometry import Point
+from qgis.core import QgsProject, QgsVectorLayer, QgsFeature, QgsGeometry, QgsWkbTypes, QgsField
+
+print("Getting road centerlines...")
+
+# Get the HOME_POINTS layer
+home_points_layer = QgsProject.instance().mapLayersByName('HOME_POINTS')[0]
+
+# Check if there's at least one feature in the layer
+if home_points_layer.featureCount() > 0:
+    # Generate a buffer around each point
+    buffer_distance = 0.01  # change to your desired buffer distance
+    buffers = []
+    for feature in home_points_layer.getFeatures():
+        point = wkt.loads(feature.geometry().asWkt())
+        buffer_polygon = point.buffer(buffer_distance)
+        buffers.append(buffer_polygon)
+
+    # Combine all buffers into a single polygon
+    union_polygon = gpd.GeoSeries(buffers).unary_union
+
+    # Create PROJECTAREA layer in memory
+    project_area_layer = QgsVectorLayer("Polygon?crs=epsg:4326", "PROJECTAREA", "memory")
+    pr = project_area_layer.dataProvider()
+
+    # Start editing the layer
+    project_area_layer.startEditing()
+
+    # Create new feature for the layer
+    f = QgsFeature()
+    f.setGeometry(QgsGeometry.fromWkt(union_polygon.wkt))
+    pr.addFeature(f)
+
+    # Commit changes
+    project_area_layer.commitChanges()
+
+    # Add the layer to the Layers panel in QGIS
+    project_area_layer.setOpacity(0.3)
+    QgsProject.instance().addMapLayer(project_area_layer)
+
+    # Use OSMnx to download the street network
+    G = ox.graph_from_polygon(union_polygon, network_type='drive')
+
+    # Convert the graph into GeoDataFrames
+    gdf_nodes, gdf_edges = ox.graph_to_gdfs(G)
+
+    # Convert gdf_edges (GeoDataFrame) to WKT format for QGIS
+    wkt_list = gdf_edges['geometry'].to_wkt().tolist()
+
+    # Create a new layer for the street centerlines
+    vl = QgsVectorLayer("LineString?crs=epsg:4326", "CENTERLINES", "memory")
+
+    # Start editing the layer
+    vl.startEditing()
+
+    # Get the data provider
+    pr = vl.dataProvider()
+
+    # Create new features for the layer
+    for wkt in wkt_list:
+        f = QgsFeature()
+        f.setGeometry(QgsGeometry.fromWkt(wkt))
+        pr.addFeature(f)
+
+    # Commit changes
+    vl.commitChanges()
+
+    # Add the layer to the Layers panel in QGIS
+    QgsProject.instance().addMapLayer(vl)
+
+else:
+    print("No points found in the HOME_POINTS layer.")
+
+print("Done.")

qgis_scripts/create_aerial_drops.py

@@ -0,0 +1,50 @@
+from qgis.PyQt.QtCore import QVariant
+from qgis.core import (QgsFeature, QgsField, QgsGeometry, QgsPointXY, QgsVectorLayer, QgsProject)
+from geopy.distance import great_circle
+
+# Constants
+SEARCH_RADIUS = 300 / 3.28084  # converting feet to meters
+COST_PER_METER = 1.5
+
+# Locate the existing layers
+home_points_layer = QgsProject.instance().mapLayersByName('HOME_POINTS')[0]
+poles_layer = QgsProject.instance().mapLayersByName('POLES')[0]
+edges_layer = QgsProject.instance().mapLayersByName('EDGES')[0]
+
+# Define the data provider
+pr = edges_layer.dataProvider()
+
+# Iterate over every home point
+for home_feature in home_points_layer.getFeatures():
+    home_point = home_feature.geometry().asPoint()
+
+    # Initialize nearest neighbor search
+    nearest_neighbor = None
+    nearest_distance = None
+
+    # Search for the nearest pole
+    for pole_feature in poles_layer.getFeatures():
+        pole_point = pole_feature.geometry().asPoint()
+        distance = great_circle((home_point.y(), home_point.x()), (pole_point.y(), pole_point.x())).meters
+
+        # If the pole is within 300 feet and it's closer than the previous nearest neighbor
+        if distance < SEARCH_RADIUS and (nearest_neighbor is None or distance < nearest_distance):
+            nearest_neighbor = pole_feature
+            nearest_distance = distance
+
+    # If a nearest neighbor was found within 300 feet
+    if nearest_neighbor is not None:
+        # Create a new edge feature
+        edge = QgsFeature()
+        edge.setGeometry(QgsGeometry.fromPolylineXY([home_point, nearest_neighbor.geometry().asPoint()]))
+        
+        # Calculate cost
+        cost = nearest_distance * COST_PER_METER
+
+        # Set attributes
+        edge.setAttributes(['Aerial Drop', nearest_distance, cost])
+        pr.addFeature(edge)
+
+# Update the layer's extent when new features have been added
+edges_layer.updateExtents()
+print("Done.")

qgis_scripts/create_aerial_edges_v2.py

@@ -0,0 +1,70 @@
+from geopy.distance import great_circle
+from qgis.core import (QgsFeature, QgsField, QgsGeometry, QgsPoint, QgsPointXY, QgsVectorLayer, QgsProject, QgsSpatialIndex)
+
+# Locate the aerial path layer
+aerial_path_layer = QgsProject.instance().mapLayersByName('AERIAL_PATH')[0]
+
+# Locate the edges layer
+edges_layer = QgsProject.instance().mapLayersByName('EDGES')[0]
+
+# Define the data provider
+pr = edges_layer.dataProvider()
+
+# Initialize spatial index
+index = QgsSpatialIndex()
+
+# List to hold new features and their IDs
+new_edges = []
+features_by_id = {}
+
+# Iterate over every feature of the AERIAL_PATH layer
+for feature in aerial_path_layer.getFeatures():
+    # Get the geometry of the feature
+    geom = feature.geometry()
+
+    # Handle both 2D and 3D geometries
+    if geom.isMultipart():
+        vertices = [v for part in geom.asMultiPolyline() for v in part]
+    else:
+        vertices = geom.asPolyline()
+
+    for i in range(1, len(vertices)):
+        # Calculate the distance between the two points in feet
+        pt1 = vertices[i - 1]
+        pt2 = vertices[i]
+        length = great_circle((pt1.y(), pt1.x()), (pt2.y(), pt2.x())).feet
+
+        # Create a new feature in the EDGES layer for each vertex pair
+        edge = QgsFeature()
+        new_geom = QgsGeometry.fromPolyline([QgsPoint(pt1), QgsPoint(pt2)])
+        edge.setGeometry(new_geom)
+        edge.setAttributes(['Aerial', length, length * 2.5])
+
+        new_edges.append(edge)
+        pr.addFeature(edge)
+        index.insertFeature(edge)
+        features_by_id[edge.id()] = edge
+
+# Start editing the edges layer
+edges_layer.startEditing()
+
+# Adjust intersecting line endpoints to the average intersection point
+for edge_id, edge in features_by_id.items():
+    endpoints = [QgsPoint(point) for point in edge.geometry().asPolyline()]
+    for i in [0, -1]:  # Check both endpoints
+        # Find nearby endpoints (within a small threshold)
+        search_radius = 0.00005
+        nearby_ids = index.intersects(QgsGeometry.fromPointXY(QgsPointXY(endpoints[i])).buffer(search_radius, 8).boundingBox())
+        nearby_endpoints = [QgsPoint(point) for id in nearby_ids for point in features_by_id[id].geometry().asPolyline() if QgsPoint(point).distance(endpoints[i]) < search_radius]
+        if len(nearby_endpoints) > 2:
+            # Update endpoint to the average nearby endpoint
+            avg_x = sum(point.x() for point in nearby_endpoints) / len(nearby_endpoints)
+            avg_y = sum(point.y() for point in nearby_endpoints) / len(nearby_endpoints)
+            endpoints[i] = QgsPoint(avg_x, avg_y)  # Create new QgsPoint object
+    # Update geometry in the layer
+    edges_layer.changeGeometry(edge.id(), QgsGeometry.fromPolyline(endpoints))
+
+# Commit the changes and update the layer's extent when new features have been added
+edges_layer.commitChanges()
+edges_layer.updateExtents()
+print("Done.")

qgis_scripts/create_network_v3.py

@@ -0,0 +1,95 @@
+# Import required libraries
+import networkx as nx
+import time
+from qgis.PyQt.QtCore import QVariant
+from qgis.core import QgsProject, QgsVectorLayer, QgsFeature, QgsField, QgsGeometry, QgsPoint
+import heapq
+
+# Get layers
+edges_layer = QgsProject.instance().mapLayersByName('EDGES')[0]
+nodes_layer = QgsProject.instance().mapLayersByName('NODES')[0]
+home_points_layer = QgsProject.instance().mapLayersByName('HOME_POINTS')[0]
+
+# Build the graph
+G = nx.Graph()
+nodes_data = {feature['id']: feature.geometry() for feature in nodes_layer.getFeatures()}  # store node geometries by id
+edges_data = {(feature['start_node'], feature['end_node']): (feature.geometry(), feature['length'], feature['cost']) for feature in edges_layer.getFeatures()}  # store edge geometries by node pairs
+
+for feature in nodes_layer.getFeatures():
+    G.add_node(feature['id'])
+
+for feature in edges_layer.getFeatures():
+    start_node = feature['start_node']  # assuming edges layer has start_node & end_node fields
+    end_node = feature['end_node']
+    G.add_edge(start_node, end_node, weight=feature['cost'], type=feature['type'], length=feature['length'], cost=feature['cost'])
+
+# Identify home nodes
+home_nodes = {feature['drop_point_id'] for feature in home_points_layer.getFeatures() if feature['drop_point_id'] in G}
+
+# Create a subgraph of G that only includes the home nodes and any nodes that connect them
+subgraph = nx.Graph()
+visited = set()
+
+# Start at any 'home' node
+start = next(iter(home_nodes))
+visited.add(start)
+
+counter = 0
+total = len(home_nodes)
+process_times = []
+
+while home_nodes.difference(visited):
+    start_time = time.time()
+
+    # Find the nearest 'home' node that is not in our tree yet
+    next_node = min(
+        (node for node in home_nodes if node not in visited),
+        key=lambda node: nx.shortest_path_length(G, start, node, weight='cost')
+    )
+
+    # Create a subgraph excluding other home nodes except start and next_node
+    G_sub = G.copy()
+    remove_nodes = home_nodes - {start} - {next_node}
+    G_sub.remove_nodes_from(remove_nodes)
+
+    # Add the shortest path to this node to our tree
+    shortest_path = nx.shortest_path(G_sub, start, next_node, weight='cost')
+    subgraph.add_edges_from(zip(shortest_path, shortest_path[1:]))
+    visited.add(next_node)
+    start = next_node
+
+    end_time = time.time()
+    process_times.append(end_time - start_time)
+    counter += 1
+    avg_time = sum(process_times) / len(process_times)
+    remaining = (total - counter) * avg_time
+    hours, remainder = divmod(remaining, 3600)
+    minutes, seconds = divmod(remainder, 60)
+    print(f'Processing home node {counter} of {total}.\nEstimated time remaining: {int(hours)} hours, {int(minutes)} minutes, and {int(seconds)} seconds')
+    
+# Create a minimum spanning tree from the subgraph
+mst = nx.minimum_spanning_tree(subgraph, weight='cost')
+
+# Prepare the NETWORK layer
+network_layer = QgsVectorLayer("LineString", "NETWORK", "memory")
+network_prov = network_layer.dataProvider()
+
+# Add 'id', 'type', 'length', and 'cost' fields to the NETWORK layer
+network_prov.addAttributes([QgsField("id", QVariant.Int), QgsField("type", QVariant.String), QgsField("length", QVariant.Double), QgsField("cost", QVariant.Double)])
+network_layer.updateFields()
+
+# Add edges from the minimum spanning tree to the NETWORK layer
+for edge in mst.edges():
+    edge_geometry, edge_length, edge_cost = edges_data[(edge[0], edge[1])] if (edge[0], edge[1]) in edges_data else edges_data[(edge[1], edge[0])]  # fetch geometry, length, cost by node pair
+    edge_type = G.get_edge_data(edge[0], edge[1])['type']  # get edge 'type' from the original graph G
+    network_feature = QgsFeature()
+    network_feature.setGeometry(edge_geometry)
+    network_feature.setAttributes([edge[0], edge_type, edge_length, edge_cost])  # set the type, length, cost attribute from edge data
+    network_prov.addFeature(network_feature)
+
+# Update NETWORK layer
+network_layer.updateExtents()
+
+# Add NETWORK layer to Layers panel
+QgsProject.instance().addMapLayer(network_layer)
+print("Done.")

qgis_scripts/create_nodes.py

@@ -0,0 +1,78 @@
+import itertools
+
+from qgis.PyQt.QtCore import QVariant
+from qgis.core import (QgsFeature, QgsField, QgsGeometry, QgsPointXY, 
+                       QgsVectorLayer, QgsProject, QgsVectorDataProvider)
+
+print("Creating nodes...")
+
+# Load the EDGES layer
+edges_layer = QgsProject.instance().mapLayersByName('EDGES')[0]
+
+# Prepare the NODES layer
+nodes_layer = QgsVectorLayer("Point", "NODES", "memory")
+prov = nodes_layer.dataProvider()
+
+# Add 'id' field to the NODES layer
+prov.addAttributes([QgsField("id", QVariant.Int)])
+nodes_layer.updateFields()
+
+# Prepare the start_node and end_node fields for the EDGES layer
+edges_layer.startEditing()
+edges_prov = edges_layer.dataProvider()
+edges_prov.addAttributes([QgsField("start_node", QVariant.Int),
+                          QgsField("end_node", QVariant.Int)])
+edges_layer.updateFields()
+
+# To store the unique nodes
+unique_nodes = []
+
+id_counter = itertools.count(start=1)
+
+node_id_map = {}
+
+total_features = edges_layer.featureCount()
+count = 0
+
+for feature in edges_layer.getFeatures():
+    count += 1
+    print(f'Processing feature {count} of {total_features} ({(count/total_features)*100:.2f}%)')
+    
+    # Get the start and end points
+    start_point = feature.geometry().asPolyline()[0]
+    end_point = feature.geometry().asPolyline()[-1]
+    
+    # Check if nodes are unique and add them to NODES layer
+    if start_point not in unique_nodes:
+        unique_nodes.append(start_point)
+        node_id = next(id_counter)
+        node_feature = QgsFeature()
+        node_feature.setGeometry(QgsGeometry.fromPointXY(start_point))
+        node_feature.setAttributes([node_id])
+        prov.addFeature(node_feature)
+        node_id_map[start_point] = node_id
+        
+    if end_point not in unique_nodes:
+        unique_nodes.append(end_point)
+        node_id = next(id_counter)
+        node_feature = QgsFeature()
+        node_feature.setGeometry(QgsGeometry.fromPointXY(end_point))
+        node_feature.setAttributes([node_id])
+        prov.addFeature(node_feature)
+        node_id_map[end_point] = node_id
+
+    # Populate the start_node and end_node fields for the EDGES layer
+    edges_layer.changeAttributeValue(feature.id(), 
+                                     edges_layer.fields().indexOf("start_node"), 
+                                     node_id_map[start_point])
+    edges_layer.changeAttributeValue(feature.id(), 
+                                     edges_layer.fields().indexOf("end_node"), 
+                                     node_id_map[end_point])
+
+# Stop editing the EDGES layer and save changes
+edges_layer.commitChanges()
+
+# Add the layer to the Layers panel
+QgsProject.instance().addMapLayer(nodes_layer)
+
+print("Done.")

qgis_scripts/create_poles_from_edges.py

@@ -0,0 +1,52 @@
+from qgis.PyQt.QtCore import QVariant
+from qgis.core import (QgsFeature, QgsField, QgsGeometry, QgsPointXY, QgsVectorLayer, QgsProject)
+
+# Locate the edges layer
+edges_layer = QgsProject.instance().mapLayersByName('EDGES')[0]
+
+# Create a new temporary layer for the poles
+poles_layer = QgsVectorLayer('Point?crs=epsg:4326', 'POLES', 'memory')
+
+# Define the data provider
+pr = poles_layer.dataProvider()
+
+# Add a new field to the new layer
+pr.addAttributes([QgsField('id', QVariant.Int)])
+poles_layer.updateFields()
+
+# Store coordinates of existing poles to avoid duplicates
+existing_poles = set()
+
+# Iterate over every feature of the EDGES layer
+for feature in edges_layer.getFeatures():
+    # Only consider edges of type 'Aerial'
+    if feature['type'] == 'Aerial':
+        # Get the geometry of the feature
+        geom = feature.geometry()
+
+        # Handle both 2D and 3D geometries
+        if geom.isMultipart():
+            vertices = [v for part in geom.asMultiPolyline() for v in part]
+        else:
+            vertices = geom.asPolyline()
+
+        for i, vertex in enumerate(vertices):
+            # Do not create a new pole if one already exists at this location
+            vertex_coordinates = (vertex[0], vertex[1])
+            if vertex_coordinates in existing_poles:
+                continue
+
+            # Create a new feature in the POLES layer for each vertex
+            pole = QgsFeature()
+            pole.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(*vertex_coordinates)))
+            pole.setAttributes([i])
+            pr.addFeature(pole)
+
+            # Remember this location so we don't create a duplicate pole
+            existing_poles.add(vertex_coordinates)
+
+# Update the layer's extent when new features have been added
+poles_layer.updateExtents()
+
+# Add the new layer to the project
+QgsProject.instance().addMapLayer(poles_layer)

qgis_scripts/create_transitions.py

@@ -0,0 +1,55 @@
+from qgis.PyQt.QtCore import QVariant
+from qgis.core import (QgsFeature, QgsField, QgsGeometry, QgsPointXY, QgsVectorLayer, QgsProject)
+from geopy.distance import great_circle
+
+# Constants
+SEARCH_RADIUS = 50
+COST_PER_METER = 19.10
+
+# Locate the existing layers
+poles_layer = QgsProject.instance().mapLayersByName('POLES')[0]
+edges_layer = QgsProject.instance().mapLayersByName('EDGES')[0]
+
+# Define the data provider
+pr = edges_layer.dataProvider()
+
+# Iterate over every pole
+for pole_feature in poles_layer.getFeatures():
+    pole_point = pole_feature.geometry().asPoint()
+
+    # Initialize nearest neighbor search
+    nearest_vertex = None
+    nearest_distance = None
+
+    # Search for the nearest vertex in underground edges
+    for edge_feature in edges_layer.getFeatures():
+        if edge_feature['type'] == 'Underground':
+            # Extract vertices
+            geom = edge_feature.geometry()
+            vertices = geom.asPolyline() if not geom.isMultipart() else [v for part in geom.asMultiPolyline() for v in part]
+
+            for vertex in vertices:
+                vertex_point = QgsPointXY(vertex[0], vertex[1])
+                distance = great_circle((pole_point.y(), pole_point.x()), (vertex_point.y(), vertex_point.x())).meters
+
+                # If the vertex is within 50 feet and it's closer than the previous nearest neighbor
+                if distance < SEARCH_RADIUS and (nearest_vertex is None or distance < nearest_distance):
+                    nearest_vertex = vertex_point
+                    nearest_distance = distance
+
+    # If a nearest vertex was found within 50 feet
+    if nearest_vertex is not None:
+        # Create a new edge feature
+        edge = QgsFeature()
+        edge.setGeometry(QgsGeometry.fromPolylineXY([pole_point, nearest_vertex]))
+        
+        # Calculate cost
+        cost = nearest_distance * COST_PER_METER
+
+        # Set attributes
+        edge.setAttributes(['Transition', nearest_distance, cost])
+        pr.addFeature(edge)
+
+# Update the layer's extent when new features have been added
+edges_layer.updateExtents()
+print("Done.")

qgis_scripts/drops_split_centerlines.py

@@ -0,0 +1,203 @@
+import geopandas as gpd
+from shapely.geometry import LineString, Point
+from shapely.affinity import translate
+from geopy.distance import great_circle
+from qgis.core import (
+    QgsProject, QgsVectorLayer, QgsFeature, QgsGeometry, 
+    QgsWkbTypes, QgsField, QgsFields
+)
+
+underground_cpf = 19.00
+buried_drop_cpf = 3.00
+
+def extend_line(line, extension_length):
+    """Extend a line by a given length."""
+    start_point = line.coords[0]
+    end_point = line.coords[-1]
+
+    # Calculate the line's direction
+    dx = end_point[0] - start_point[0]
+    dy = end_point[1] - start_point[1]
+
+    # Calculate the extension length along each axis
+    length = math.sqrt(dx**2 + dy**2)
+    extension_dx = dx / length * extension_length
+    extension_dy = dy / length * extension_length
+
+    # Create the extended line by translating the end point
+    extended_line = LineString([start_point, translate(line, extension_dx, extension_dy).coords[-1]])
+
+    return extended_line
+
+def closest_point_on_line(point, lines):
+    """Find closest point on the closest line to the given point."""
+    closest_line = min(lines, key=point.distance)
+    return closest_line.interpolate(closest_line.project(point))
+
+# Load your layers
+homes_layer = QgsProject.instance().mapLayersByName('HOME_POINTS')[0]
+centerlines_layer = QgsProject.instance().mapLayersByName('CENTERLINES')[0]
+
+# Convert layers to GeoDataFrames
+gdf_homes = gpd.GeoDataFrame.from_features([f for f in homes_layer.getFeatures()], crs=homes_layer.crs().toWkt())
+gdf_centerlines = gpd.GeoDataFrame.from_features([f for f in centerlines_layer.getFeatures()], crs=centerlines_layer.crs().toWkt())
+
+num_homes = len(gdf_homes)
+print("Number of home points:", num_homes)
+
+# Create an in-memory layer to store the connecting lines
+vl = QgsVectorLayer("LineString?crs=epsg:4326", "DROPS", "memory")
+vl_ext = QgsVectorLayer("LineString?crs=epsg:4326", "DROPS_EXT", "memory")
+pr = vl.dataProvider()
+pr_ext = vl_ext.dataProvider()
+
+print("Creating drops...")
+
+# For each home, find the nearest point on the centerlines and create a line
+for idx, home in gdf_homes.iterrows():
+    nearest_point = closest_point_on_line(home.geometry, gdf_centerlines.geometry)
+    line = LineString([(home.geometry.x, home.geometry.y), (nearest_point.x, nearest_point.y)])
+    
+    # Add the line to the "DROPS" layer
+    feat = QgsFeature()
+    feat.setGeometry(QgsGeometry.fromWkt(line.wkt))
+    pr.addFeature(feat)
+    
+    # Extend the line by 5% of its length beyond the nearest point and add it to the "DROPS_EXT" layer
+    line_ext = extend_line(line, line.length*0.05)
+    feat_ext = QgsFeature()
+    feat_ext.setGeometry(QgsGeometry.fromWkt(line_ext.wkt))
+    pr_ext.addFeature(feat_ext)
+    
+    # After processing each home, print the progress
+    progress = (idx + 1) / num_homes * 100  # Progress as percentage
+    print(f'\rProgress: {progress:.2f}%', end='')  # Print the progress
+
+# Update the layer's extent and add it to the project
+vl.updateExtents()
+QgsProject.instance().addMapLayer(vl)
+
+print("\nSplitting centerlines at drops...")
+
+# Define the splitting function
+def split_with_lines(input_layer, split_layer, output):
+    """
+    This function will use native:splitwithlines processing algorithm 
+    from QGIS's processing toolbox to split features in one layer using lines 
+    from another layer.
+    """
+    params = {
+        'INPUT': input_layer,
+        'LINES': split_layer,
+        'OUTPUT': output
+    }
+
+    res = processing.run("native:splitwithlines", params)
+    return res
+
+# Use the 'DROPS_EXT' layer to split features in the 'CENTERLINES' layer
+centerlines_split = split_with_lines(centerlines_layer, vl_ext, 'memory:')
+# Add the split layer to the project instance
+split_layer = centerlines_split['OUTPUT']
+#split_layer.setName('EDGES')
+#QgsProject.instance().addMapLayer(split_layer)
+
+print("Deleting duplicate edges...")
+
+# Define the delete duplicates function
+def delete_duplicates(input_layer, output):
+    """
+    This function will use qgis:deleteduplicategeometries processing algorithm 
+    from QGIS's processing toolbox to remove duplicate geometries in a layer.
+    """
+    params = {
+        'INPUT': input_layer,
+        'OUTPUT': output
+    }
+
+    res = processing.run("qgis:deleteduplicategeometries", params)
+    return res
+
+# Use the delete_duplicates function on the split layer
+split_layer_no_duplicates = delete_duplicates(split_layer, 'memory:')
+
+# Add the split layer without duplicates to the project instance
+split_layer_nd = split_layer_no_duplicates['OUTPUT']
+
+# Get the data provider
+dp = split_layer_nd.dataProvider()
+
+# Add new fields
+dp.addAttributes([
+    QgsField("type", QVariant.String),
+    QgsField("length", QVariant.Double),
+    QgsField("cost", QVariant.Double)
+])
+
+# Update to apply the changes
+split_layer_nd.updateFields()
+
+# Now you have to update all features to have 'type' as 'Underground'
+# Start editing
+split_layer_nd.startEditing()
+
+for feature in split_layer_nd.getFeatures():
+    # set type to 'Underground'
+    feature.setAttribute('type', 'Underground')
+    
+    # set length to the great circle distance of the geometry
+    # convert coordinates to tuples (latitude, longitude)
+    start_point = feature.geometry().asPolyline()[0]
+    end_point = feature.geometry().asPolyline()[-1]
+    start_coordinates = (start_point.y(), start_point.x())
+    end_coordinates = (end_point.y(), end_point.x())
+    
+    # calculate great circle distance in feet
+    length_in_feet = great_circle(start_coordinates, end_coordinates).feet
+    feature.setAttribute('length', length_in_feet)
+    
+    # set cost to a computed value, for example length * underground_cpf
+    feature.setAttribute('cost', length_in_feet * underground_cpf)
+    
+    # update the feature in the layer
+    split_layer_nd.updateFeature(feature)
+    
+# Convert 'DROPS' layer to GeoDataFrame
+gdf_drops = gpd.GeoDataFrame.from_features([f for f in vl.getFeatures()], crs=vl.crs().toWkt())
+
+# Add the 'DROPS' to the 'EDGES'
+for idx, drop in gdf_drops.iterrows():
+    # Create a new feature for the drop
+    drop_feature = QgsFeature(dp.fields())
+
+    # Convert drop's Shapely geometry to QgsGeometry
+    qgis_geom = QgsGeometry.fromWkt(drop.geometry.wkt)
+    drop_feature.setGeometry(qgis_geom)
+    
+    # Add the drop feature to the 'EDGES' layer
+    dp.addFeature(drop_feature)
+
+    # Set 'type' to 'Buried Drop'
+    drop_feature.setAttribute('type', 'Buried Drop')
+
+    # Calculate the great circle length of the drop in feet
+    start_point = drop.geometry.coords[0]
+    end_point = drop.geometry.coords[-1]
+    start_coordinates = (start_point[1], start_point[0])  # Note the reverse order (y, x) for geopy
+    end_coordinates = (end_point[1], end_point[0])  # Note the reverse order (y, x) for geopy
+    length_in_feet = great_circle(start_coordinates, end_coordinates).feet
+    drop_feature.setAttribute('length', length_in_feet)
+
+    # Calculate cost as length * buried_drop_cpf
+    drop_feature.setAttribute('cost', length_in_feet * buried_drop_cpf)
+
+    # Update the feature in the layer
+    split_layer_nd.updateFeature(drop_feature)
+
+# Commit changes
+split_layer_nd.commitChanges()
+
+split_layer_nd.setName('EDGES')
+QgsProject.instance().addMapLayer(split_layer_nd)
+
+print("Done.")

report.py

@@ -0,0 +1,76 @@
+#!/usr/bin/env python3
+
+import geopandas as gpd
+import pandas as pd
+
+# Load the network shapefile
+network_gdf = gpd.read_file('network.shp')
+
+# Load the home points shapefile
+home_points_gdf = gpd.read_file('home_points.shp')
+
+# Initialize a dictionary to hold the aggregated network data
+report_data = {}
+
+# Check if 'unit_count' exists in home_points_gdf
+include_unit_count = 'unit_count' in home_points_gdf.columns
+if include_unit_count:
+    print("unit_count column found in home_points.shp.")
+else:
+    print("unit_count column not found in home_points.shp.")
+
+# Count homes per fdh_id from home_points.shp
+home_counts = home_points_gdf['fdh_id'].value_counts().to_dict()
+
+# Aggregate Unit_count per fdh_id if available
+unit_counts = home_points_gdf.groupby('fdh_id')['unit_count'].sum().to_dict() if include_unit_count else {}
+
+for _, row in network_gdf.iterrows():
+    fdh_id = row['fdh_id']
+    type = row['type']
+    length = row.geometry.length
+
+    # Initialize the fdh_id entry if not already present
+    if fdh_id not in report_data:
+        report_data[fdh_id] = {'FDH ID': fdh_id, 'HP': 0, 'HHP': 0, 'Aerial Drop': 0, 'Buried Drop': 0, 'Aerial': 0, 'Underground': 0}
+        # Add home count if available
+        report_data[fdh_id]['HP'] = home_counts.get(fdh_id, 0)
+        # Add unit count if available
+        if include_unit_count:
+            report_data[fdh_id]['HHP'] = unit_counts.get(fdh_id, 0)
+
+    # Process network data
+    if type in ['Aerial Drop', 'Buried Drop']:
+        report_data[fdh_id][type] += 1  # Increment the count for drop types
+    elif type in ['Aerial', 'Underground', 'Transition']:
+        adjusted_type = 'Underground' if type == 'Transition' else type
+        report_data[fdh_id][adjusted_type] += length  # Add length, including 'Transition' to 'Underground'
+
+# Calculate additional columns and round lengths
+for fdh_id, data in report_data.items():
+    aerial = data['Aerial']
+    underground = data['Underground']
+    total_length = aerial + underground
+    
+    # Determine the divisor for FPP calculation based on the availability of HHP or HP
+    if include_unit_count and data['HHP'] > 0:  # If HHP is available and greater than 0
+        divisor = data['HHP']
+    else:  # If HHP is not available or 0, use HP
+        divisor = data['HP']
+
+    data['Aerial'] = round(aerial)
+    data['Underground'] = round(underground)
+    data['% Aerial'] = round((aerial / total_length * 100) if total_length > 0 else 0, 2)  # Calculate % Aerial
+    data['FPP'] = round((total_length / divisor) if divisor > 0 else 0, 2)  # Calculate Feet per Home Point
+
+# Convert the data to a pandas DataFrame and sort by FDH ID
+report_df = pd.DataFrame(list(report_data.values())).sort_values(by='FDH ID')
+
+# Specify the column order, including the new 'HHP' column if applicable
+columns_order = ['FDH ID', 'HP', 'HHP', 'Aerial Drop', 'Buried Drop', 'Aerial', 'Underground', '% Aerial', 'FPP'] if include_unit_count else ['FDH ID', 'HP', 'Aerial Drop', 'Buried Drop', 'Aerial', 'Underground', '% Aerial', 'FPP']
+report_df = report_df[columns_order]
+
+# Write the DataFrame to an Excel file
+report_df.to_excel('network_report.xlsx', index=False, engine='openpyxl')
+
+print("Sorted report with additional metrics has been saved to network_report.xlsx.")

reproject.py

@@ -0,0 +1,22 @@
+#!/usr/bin/env python3
+
+import sys
+import geopandas as gpd
+
+def reproject_shapefile(input_shapefile, output_shapefile, target_epsg):
+    gdf = gpd.read_file(input_shapefile)
+    gdf = gdf.to_crs(epsg=target_epsg)
+    gdf.to_file(output_shapefile)
+
+    print(f"Shapefile has been reprojected to EPSG:{target_epsg} and saved as {output_shapefile}")
+
+if __name__ == '__main__':
+    if len(sys.argv) != 4:
+        print("Usage: python reproject.py <input_shapefile> <output_shapefile> <target_epsg_code>")
+        sys.exit(1)
+    
+    input_shapefile = sys.argv[1]
+    output_shapefile = sys.argv[2]
+    target_epsg = sys.argv[3]  # EPSG code as a string is fine here, GeoPandas handles the conversion
+
+    reproject_shapefile(input_shapefile, output_shapefile, target_epsg)

requirements.txt

@@ -0,0 +1,40 @@
+attrs==23.2.0
+branca==0.7.1
+certifi==2024.2.2
+charset-normalizer==3.3.2
+click==8.1.7
+click-plugins==1.1.1
+cligj==0.7.2
+contourpy==1.2.0
+cycler==0.12.1
+et-xmlfile==1.1.0
+fiona==1.9.5
+folium==0.16.0
+fonttools==4.49.0
+geographiclib==2.0
+geopandas==0.14.3
+geopy==2.4.1
+idna==3.6
+Jinja2==3.1.3
+kiwisolver==1.4.5
+MarkupSafe==2.1.5
+matplotlib==3.8.3
+networkx==3.2.1
+numpy==1.26.4
+openpyxl==3.1.2
+osmnx==1.9.1
+packaging==23.2
+pandas==2.2.1
+pillow==10.2.0
+pyparsing==3.1.1
+pyproj==3.6.1
+python-dateutil==2.9.0
+pytz==2024.1
+requests==2.31.0
+scipy==1.12.0
+setuptools==69.1.1
+shapely==2.0.3
+six==1.16.0
+tzdata==2024.1
+urllib3==2.2.1
+xyzservices==2023.10.1