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.")