diff --git a/README.md b/README.md index 3370b60..0c9349d 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,61 @@ # pyhld -Fiber to the home high level autodesign in Python using shapefiles for input/output. \ No newline at end of file +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 \ No newline at end of file diff --git a/associate_drop_points.py b/associate_drop_points.py new file mode 100755 index 0000000..789eccd --- /dev/null +++ b/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.") diff --git a/associate_drop_points_v2.py b/associate_drop_points_v2.py new file mode 100755 index 0000000..8ae5ab3 --- /dev/null +++ b/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.") diff --git a/centerlines_from_homes.py b/centerlines_from_homes.py new file mode 100755 index 0000000..d92259c --- /dev/null +++ b/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.") diff --git a/check_graph.py b/check_graph.py new file mode 100755 index 0000000..12a06b6 --- /dev/null +++ b/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}") diff --git a/cluster_fdh.py b/cluster_fdh.py new file mode 100755 index 0000000..04a2d64 --- /dev/null +++ b/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') \ No newline at end of file diff --git a/cluster_fdh_v2.py b/cluster_fdh_v2.py new file mode 100755 index 0000000..c07c4b7 --- /dev/null +++ b/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') diff --git a/connect_poles.py b/connect_poles.py new file mode 100755 index 0000000..82354d8 --- /dev/null +++ b/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.") diff --git a/create_aerial_drops.py b/create_aerial_drops.py new file mode 100755 index 0000000..a5b3905 --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/create_aerial_edges.py b/create_aerial_edges.py new file mode 100755 index 0000000..bc2d195 --- /dev/null +++ b/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.") diff --git a/create_map.py b/create_map.py new file mode 100755 index 0000000..5c91a66 --- /dev/null +++ b/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.") diff --git a/create_map_buried.py b/create_map_buried.py new file mode 100755 index 0000000..5ddd308 --- /dev/null +++ b/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.") diff --git a/create_mst_clusters.py b/create_mst_clusters.py new file mode 100755 index 0000000..d111d67 --- /dev/null +++ b/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.") diff --git a/create_mst_clusters_v2.py b/create_mst_clusters_v2.py new file mode 100755 index 0000000..8fe529d --- /dev/null +++ b/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') \ No newline at end of file diff --git a/create_network.py b/create_network.py new file mode 100755 index 0000000..800ce27 --- /dev/null +++ b/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') diff --git a/create_network_v2.1.py b/create_network_v2.1.py new file mode 100755 index 0000000..2b8e26e --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/create_network_v2.2.py b/create_network_v2.2.py new file mode 100755 index 0000000..7365baf --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/create_network_v2.py b/create_network_v2.py new file mode 100755 index 0000000..08ceef5 --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/create_network_v3.py b/create_network_v3.py new file mode 100755 index 0000000..306f8ab --- /dev/null +++ b/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.") diff --git a/create_network_v4.py b/create_network_v4.py new file mode 100755 index 0000000..5f0f4c8 --- /dev/null +++ b/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.") diff --git a/create_nodes.py b/create_nodes.py new file mode 100755 index 0000000..1499317 --- /dev/null +++ b/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.") diff --git a/create_transitions.py b/create_transitions.py new file mode 100755 index 0000000..b5fc786 --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/design_from_roads.sh b/design_from_roads.sh new file mode 100755 index 0000000..65726ae --- /dev/null +++ b/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 \ No newline at end of file diff --git a/design_from_start.sh b/design_from_start.sh new file mode 100755 index 0000000..d6251e8 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/design_underground_only.sh b/design_underground_only.sh new file mode 100755 index 0000000..adb7b86 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/drops_split_centerlines.py b/drops_split_centerlines.py new file mode 100755 index 0000000..f15cf0a --- /dev/null +++ b/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.") diff --git a/drops_split_centerlines_v2.py b/drops_split_centerlines_v2.py new file mode 100755 index 0000000..022a4b7 --- /dev/null +++ b/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.") diff --git a/exchange_boundary.py b/exchange_boundary.py new file mode 100755 index 0000000..0313c65 --- /dev/null +++ b/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 ") + sys.exit(1) + + input_shapefile = sys.argv[1] + output_shapefile = 'exchange_boundary.shp' + create_concave_hull_shapefile(input_shapefile, output_shapefile) diff --git a/open_all_kml.py b/open_all_kml.py new file mode 100755 index 0000000..9232086 --- /dev/null +++ b/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) diff --git a/poles_used.py b/poles_used.py new file mode 100755 index 0000000..83f91ec --- /dev/null +++ b/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.") + diff --git a/qc/drop_endpoints.py b/qc/drop_endpoints.py new file mode 100644 index 0000000..d9ca1e9 --- /dev/null +++ b/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() diff --git a/qc/drop_endpoints_folium.py b/qc/drop_endpoints_folium.py new file mode 100644 index 0000000..2969225 --- /dev/null +++ b/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') diff --git a/qgis_scripts/associate_drop_points.py b/qgis_scripts/associate_drop_points.py new file mode 100644 index 0000000..c6e1219 --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/qgis_scripts/centerlines_from_homes.py b/qgis_scripts/centerlines_from_homes.py new file mode 100644 index 0000000..25a9e28 --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/qgis_scripts/create_aerial_drops.py b/qgis_scripts/create_aerial_drops.py new file mode 100644 index 0000000..1adb4a8 --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/qgis_scripts/create_aerial_edges_v2.py b/qgis_scripts/create_aerial_edges_v2.py new file mode 100644 index 0000000..3026a72 --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/qgis_scripts/create_network_v3.py b/qgis_scripts/create_network_v3.py new file mode 100644 index 0000000..6ce74eb --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/qgis_scripts/create_nodes.py b/qgis_scripts/create_nodes.py new file mode 100644 index 0000000..f0883a9 --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/qgis_scripts/create_poles_from_edges.py b/qgis_scripts/create_poles_from_edges.py new file mode 100644 index 0000000..f6dba07 --- /dev/null +++ b/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) diff --git a/qgis_scripts/create_transitions.py b/qgis_scripts/create_transitions.py new file mode 100644 index 0000000..492ccd2 --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/qgis_scripts/drops_split_centerlines.py b/qgis_scripts/drops_split_centerlines.py new file mode 100644 index 0000000..cdbf7f7 --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/report.py b/report.py new file mode 100755 index 0000000..5aa3da7 --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/reproject.py b/reproject.py new file mode 100755 index 0000000..ec39929 --- /dev/null +++ b/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 ") + 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) diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..74e48e3 --- /dev/null +++ b/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