modify this script so that the lines that are used for routing are drawn on either side of the road favouring either one side or the other. Both sides can be drawn but only as a last resort in order to get to all the centroids. Centroids can not be connected to these lines by crossing over the roads line. these represent trench lines and should be closer to the erven polygons. import geopandas as gpd
import osmnx as ox
import numpy as np
from shapely.geometry import Point, LineString
from shapely.ops import unary_union
import networkx as nx
from fastkml import kml

# Define verbose logging function
def log(message):
    print(message)

# Load the input shapefile
input_layer_path = r"D:\Piketberg Erven.shp"
log("Loading shapefile...")
erven_gdf = gpd.read_file(input_layer_path)

# Generate centroids for all the polygons
log("Calculating centroids...")
erven_gdf['centroid'] = erven_gdf.geometry.centroid
centroids_gdf = gpd.GeoDataFrame(geometry=erven_gdf['centroid'], crs=erven_gdf.crs)

# Create a buffered area around the polygons
log("Creating buffered area...")
buffered_gdf = erven_gdf.buffer(0.00005)
dissolved_buffer = gpd.GeoDataFrame(geometry=[unary_union(buffered_gdf.geometry)], crs=erven_gdf.crs)

# Get the bounds for OSM data download
bounds = dissolved_buffer.total_bounds  # (west, south, east, north)
west, south, east, north = bounds

# Download road network data within the buffered area
try:
    log("Downloading road network data from OSM...")
    G = ox.graph_from_bbox(north, south, east, west, network_type='drive')
    log("OSM data successfully loaded.")
except ox._errors.InsufficientResponseError:
    log("No data found for the specified area. Please check the bounding box or network type.")
    exit(1)

# Convert graph to GeoDataFrames
node_gdf, edge_gdf = ox.graph_to_gdfs(G)

# Get dissolved polygon edges as network lines
log("Extracting dissolved polygon edges...")
dissolved_lines = edge_gdf[edge_gdf.geometry.intersects(dissolved_buffer.geometry[0])]

# Create aggregation nodes (boundary boxes)
log("Creating boundary boxes for aggregation nodes...")
boundary_boxes = []
for idx in range(0, len(erven_gdf), 2):  # One box for every two erven polygons
    box_geometry = erven_gdf.geometry.iloc[idx:idx+2].unary_union.envelope
    boundary_boxes.append(box_geometry)

boundary_boxes_gdf = gpd.GeoDataFrame(geometry=boundary_boxes, crs=erven_gdf.crs)

# Place a central node
log("Calculating central point...")
central_point = Point(boundary_boxes_gdf.geometry.unary_union.centroid.x, 
                      boundary_boxes_gdf.geometry.unary_union.centroid.y)

# Create an aggregation node GeoDataFrame
node_gdf = gpd.GeoDataFrame(geometry=[central_point], crs=erven_gdf.crs)
node_gdf['label'] = 'Central Node'

# Snap centroids and the central point to the nearest nodes in the graph
def snap_to_nearest_node(point, graph):
    nearest_node = ox.distance.nearest_nodes(graph, point.x, point.y)
    return nearest_node

# Create connections to nearest boundary boxes and number them
log("Creating connections to nearest boundary boxes...")
connections = []
boundary_box_labels = []

for i, boundary in enumerate(boundary_boxes_gdf.geometry):
    label = f'BB{str(i+1).zfill(3)}'  # Format as BB001, BB002, etc.
    boundary_box_labels.append(label)
    
    for centroid in centroids_gdf.geometry:
        connections.append(LineString([centroid, boundary.centroid]))

# Create a GeoDataFrame for connections
connections_gdf = gpd.GeoDataFrame(geometry=connections, crs=erven_gdf.crs)

# Find the nearest nodes for each centroid and the central point
log("Finding nearest nodes for centroids...")
centroid_nodes = [snap_to_nearest_node(centroid, G) for centroid in centroids_gdf.geometry]
central_node = snap_to_nearest_node(central_point, G)

# Check if nodes are valid
valid_centroid_nodes = [node for node in centroid_nodes if node in G.nodes()]
if not valid_centroid_nodes:
    log("No valid centroid nodes found within the network.")
    exit(1)

# Find shortest paths from each valid centroid node back to the central node, prioritizing shared paths
log("Calculating shortest paths from centroids to the central node...")
shortest_paths = {}
shared_paths = []

for centroid_node in valid_centroid_nodes:
    try:
        path = nx.shortest_path(G, source=centroid_node, target=central_node)
        shortest_paths[centroid_node] = path
        
        # Adding shared paths logic
        shared_paths.extend(path)
    except nx.NetworkXNoPath:
        log(f"No path found from centroid node {centroid_node} to central node {central_node}.")

# Create KML output
log("Creating KML output...")
k = kml.KML()
document = kml.Document()
k.append(document)

# Add central node to KML
placemark_central = kml.Placemark()
placemark_central.name = "Central Node"
placemark_central.geometry = Point(central_point.x, central_point.y)
document.append(placemark_central)

# Add connections to KML with boundary boxes
for idx, boundary in enumerate(boundary_boxes):
    placemark_boundary = kml.Placemark()
    placemark_boundary.name = boundary_box_labels[idx]
    placemark_boundary.geometry = Point(boundary.centroid.x, boundary.centroid.y)
    document.append(placemark_boundary)

# Add connections to KML
for centroid, path in shortest_paths.items():
    coords = [(G.nodes[node]['x'], G.nodes[node]['y']) for node in path]
    line = LineString(coords)
    placemark_line = kml.Placemark()
    placemark_line.name = f"Path from Node {centroid} to Central Node"
    placemark_line.geometry = line
    document.append(placemark_line)

# Add centroid points to KML
for idx, centroid in enumerate(centroids_gdf.geometry):
    placemark_centroid = kml.Placemark()
    placemark_centroid.name = f"Centroid {idx+1}"
    placemark_centroid.geometry = centroid
    document.append(placemark_centroid)

# Save KML
output_kml_path = "output.kml"
with open(output_kml_path, "w") as f:
    f.write(k.to_string(prettyprint=True))

log(f"KML output saved to '{output_kml_path}'")