the boundary boxes and buffered area output is still wrong. There is no buffer for this and the Boundary boxes are POINTS not BOXES. Make sure that instead of drawing boxes and buffering useless data, the script is placing BBs 1 for every 2 ervens on the lines of the buffered polygon edges. Place the points preferably as close as possible to where there are two ervens touching sides. import geopandas as gpd
import osmnx as ox
import numpy as np
from shapely.geometry import Point, LineString
from shapely.ops import unary_union, nearest_points
import networkx as nx
from fastkml import kml
import warnings
import os
import sys
import matplotlib.pyplot as plt

# Suppress specific warnings
warnings.filterwarnings("ignore", category=FutureWarning)

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

# Function to visualize geometries
def visualize_geometries(geometries, title):
    plt.figure(figsize=(10, 10))
    for geom in geometries:
        if geom.type == 'Polygon':
            x, y = geom.exterior.xy
            plt.fill(x, y, alpha=0.5, fc='blue', ec='black')
        elif geom.type == 'MultiPolygon':
            for poly in geom.geoms:
                x, y = poly.exterior.xy
                plt.fill(x, y, alpha=0.5, fc='blue', ec='black')
        elif geom.type == 'LineString':
            x, y = geom.xy
            plt.plot(x, y, color='red', linewidth=2)
    plt.title(title)
    plt.xlabel("Longitude")
    plt.ylabel("Latitude")
    plt.axis("equal")
    plt.grid()
    plt.show()

# Parameters
SIDE_PREFERENCE = 'left'  # Options: 'left', 'right'
OFFSET_DISTANCE = 2        # Meters; adjust as needed based on CRS
BUFFER_DISTANCE = 5        # Meters; adjust as needed based on CRS
TARGET_UTM_EPSG = 32734  # EPSG code for WGS 84 / UTM zone 34S

# Define shapefile path
input_layer_path = r"D:\Spatial\Piketberg Erven.shp"

# Check if all necessary shapefile components are present
required_extensions = ['.shp', '.shx', '.dbf']
basename, _ = os.path.splitext(input_layer_path)
missing_files = []
for ext in required_extensions:
    file_path = basename + ext
    if not os.path.exists(file_path):
        missing_files.append(file_path)

if missing_files:
    log("Error: The following essential shapefile components are missing:")
    for file in missing_files:
        log(f" - {file}")
    log("Please ensure all necessary shapefile components are present.")
    sys.exit(1)
else:
    log("All essential shapefile components are present.")

# Attempt to read the shapefile in latitude/longitude (EPSG:4326)
try:
    log("Loading shapefile in EPSG:4326...")
    erven_gdf_latlon = gpd.read_file(input_layer_path)  # Keep original in lat/lon
    log("Shapefile loaded successfully.")
except Exception as e:
    log(f"Error loading shapefile: {e}")
    sys.exit(1)

# Check if GeoDataFrame is empty
if erven_gdf_latlon.empty:
    log("Error: The input shapefile is empty. Please check the file contents.")
    sys.exit(1)

# Validate geometries
log("Validating geometries...")
erven_gdf_latlon['geometry'] = erven_gdf_latlon['geometry'].buffer(0)  # Attempt to fix invalid geometries
invalid_geoms = erven_gdf_latlon[~erven_gdf_latlon.is_valid]
if not invalid_geoms.empty:
    log(f"Found {len(invalid_geoms)} invalid geometries. Attempting to fix...")
    erven_gdf_latlon.loc[~erven_gdf_latlon.is_valid, 'geometry'] = erven_gdf_latlon.loc[~erven_gdf_latlon.is_valid, 'geometry'].buffer(0)
    # Recheck
    invalid_geoms = erven_gdf_latlon[~erven_gdf_latlon.is_valid]
    if not invalid_geoms.empty:
        log("Some geometries could not be fixed. Please review the input data.")
        sys.exit(1)
    else:
        log("All geometries are now valid.")
else:
    log("All geometries are valid.")

# Ensure CRS is set
if erven_gdf_latlon.crs is None:
    log("Error: CRS is not set for the shapefile. Please assign the correct CRS before proceeding.")
    sys.exit(1)
else:
    log(f"Shapefile CRS: {erven_gdf_latlon.crs}")

# Create a projected version for buffering and spatial operations
log(f"Projecting GeoDataFrame to UTM Zone 34S (EPSG:{TARGET_UTM_EPSG}) for accurate distance operations...")
try:
    erven_gdf_utm = erven_gdf_latlon.to_crs(epsg=TARGET_UTM_EPSG)
    log("Projection successful.")
except Exception as e:
    log(f"Projection failed: {e}")
    sys.exit(1)

# Generate centroids for all the polygons in UTM
log("Calculating centroids...")
erven_gdf_utm['centroid'] = erven_gdf_utm.geometry.centroid
centroids_gdf_utm = gpd.GeoDataFrame(geometry=erven_gdf_utm['centroid'], crs=erven_gdf_utm.crs)
log(f"Calculated {len(centroids_gdf_utm)} centroids.")

# Create a buffered area around the polygons in UTM
log("Creating buffered area...")
buffered_gdf_utm = erven_gdf_utm.buffer(BUFFER_DISTANCE)
dissolved_buffer_utm = unary_union(buffered_gdf_utm)
log("Buffered area created and dissolved.")

# Visualize the buffered areas
visualize_geometries([dissolved_buffer_utm, *erven_gdf_utm.geometry], title="Buffered Area and Polygons")

# Reproject the dissolved buffer back to lat/lon for OSMnx bounding box
log("Reprojecting dissolved buffer back to EPSG:4326 for OSMnx...")
dissolved_buffer_latlon = gpd.GeoSeries([dissolved_buffer_utm], crs=erven_gdf_utm.crs).to_crs(epsg=4326)
dissolved_buffer_latlon = dissolved_buffer_latlon.iloc[0]

# Get the bounds in latitude and longitude for OSM data download
bounds = dissolved_buffer_latlon.bounds  # (minx, miny, maxx, maxy) in lat/lon
minx, miny, maxx, maxy = bounds
log(f"Buffered area bounds in lat/lon: minx={minx}, miny={miny}, maxx={maxx}, maxy={maxy}")

# Download road network data within the buffered area in lat/lon
try:
    log("Downloading road network data from OSM...")
    roads_G = ox.graph_from_bbox(maxy, miny, maxx, minx, network_type='drive', simplify=True)
    log("OSM road data successfully loaded.")
except ox._errors.InsufficientResponseError:
    log("Error: No road data found for the specified area. Please check the bounding box or network type.")
    sys.exit(1)
except Exception as e:
    log(f"Error downloading road network data: {e}")
    sys.exit(1)

# Convert roads graph to GeoDataFrame
try:
    roads_node_gdf_latlon, roads_edge_gdf_latlon = ox.graph_to_gdfs(roads_G)
    log("Converted road network graph to GeoDataFrames.")
except Exception as e:
    log(f"Error converting road network to GeoDataFrames: {e}")
    sys.exit(1)

# Validate road geometries
log("Validating road geometries...")
roads_edge_gdf_latlon['geometry'] = roads_edge_gdf_latlon['geometry'].buffer(0)  # Attempt to fix invalid geometries
invalid_roads = roads_edge_gdf_latlon[~roads_edge_gdf_latlon.is_valid]
if not invalid_roads.empty:
    log(f"Found {len(invalid_roads)} invalid road geometries. Attempting to fix...")
    roads_edge_gdf_latlon.loc[~roads_edge_gdf_latlon.is_valid, 'geometry'] = roads_edge_gdf_latlon.loc[~roads_edge_gdf_latlon.is_valid, 'geometry'].buffer(0)
    # Recheck
    invalid_roads = roads_edge_gdf_latlon[~roads_edge_gdf_latlon.is_valid]
    if not invalid_roads.empty:
        log("Some road geometries could not be fixed. Please review the OSM data.")
        sys.exit(1)
    else:
        log("All road geometries are now valid.")
else:
    log("All road geometries are valid.")

# Reproject roads to UTM for spatial operations
log(f"Projecting road network to UTM Zone 34S (EPSG:{TARGET_UTM_EPSG})...")
try:
    roads_node_gdf = roads_node_gdf_latlon.to_crs(epsg=TARGET_UTM_EPSG)
    roads_edge_gdf = roads_edge_gdf_latlon.to_crs(epsg=TARGET_UTM_EPSG)
    log("Projection of road network successful.")
except Exception as e:
    log(f"Projection of road network failed: {e}")
    sys.exit(1)

# Get dissolved polygon edges as network lines in UTM
log("Extracting dissolved polygon edges...")
dissolved_polygons_utm = unary_union(erven_gdf_utm.geometry)

# Debugging Information
log(f"Type of dissolved_polygons_utm: {type(dissolved_polygons_utm)}")
log(f"Geometry type: {dissolved_polygons_utm.type}")

if dissolved_polygons_utm.type == 'Polygon':
    dissolved_lines = [LineString(dissolved_polygons_utm.exterior)]
    log("Processed as a single Polygon.")
elif dissolved_polygons_utm.type == 'MultiPolygon':
    dissolved_lines = [LineString(poly.exterior) for poly in dissolved_polygons_utm.geoms]
    log(f"Processed MultiPolygon with {len(dissolved_polygons_utm.geoms)} individual Polygons.")
else:
    log(f"Error: Encountered unsupported geometry type: {dissolved_polygons_utm.type}")
    sys.exit(1)

# Create a GeoDataFrame from the dissolved lines
dissolved_lines_gdf = gpd.GeoDataFrame(geometry=dissolved_lines, crs=erven_gdf_utm.crs)
log(f"Extracted {len(dissolved_lines_gdf)} dissolved polygon edge(s).")

# Visualize the dissolved lines
visualize_geometries(dissolved_lines, title="Dissolved Polygon Edges")

# Create a NetworkX graph from dissolved polygon edges
log("Creating trench network graph from dissolved polygon edges...")
trench_G = nx.Graph()

for idx, row in dissolved_lines_gdf.iterrows():
    line = row.geometry
    if line is None or line.is_empty:
        continue
    coords = list(line.coords)
    for u, v in zip(coords[:-1], coords[1:]):
        trench_G.add_edge(u, v, geometry=LineString([u, v]), length=Point(u).distance(Point(v)))
log(f"Trench network graph created with {trench_G.number_of_nodes()} nodes and {trench_G.number_of_edges()} edges.")

# Create aggregation nodes (boundary boxes) in UTM
log("Creating boundary boxes for aggregation nodes...")
boundary_boxes = []
for idx in range(0, len(erven_gdf_utm), 2):  # One box for every two erven polygons
    subset = erven_gdf_utm.geometry.iloc[idx:idx + 2]
    union = unary_union(subset)
    # Buffer slightly outward to avoid overlapping with erven
    box_geometry = union.buffer(2).envelope  # Buffer of 2 meters; adjust as needed
    boundary_boxes.append(box_geometry)

boundary_boxes_gdf_utm = gpd.GeoDataFrame(geometry=boundary_boxes, crs=erven_gdf_utm.crs)
log(f"Created {len(boundary_boxes_gdf_utm)} boundary box(es).")

# Visualize boundary boxes
visualize_geometries([dissolved_buffer_utm, *erven_gdf_utm.geometry, *boundary_boxes], title="Boundary Boxes and Buffer Area")

# Place a central node in UTM
log("Calculating central point...")
central_point_utm = boundary_boxes_gdf_utm.geometry.unary_union.centroid
log(f"Central point at: ({central_point_utm.x}, {central_point_utm.y})")

# Add central node to trench graph
trench_G.add_node((central_point_utm.x, central_point_utm.y))
log("Central node added to trench network graph.")

# Connecting boundary boxes to the central node in trench graph
log("Connecting boundary boxes to the central node...")
for boundary in boundary_boxes_gdf_utm.geometry:
    # Find the nearest point on trench network to the boundary centroid
    nearest_geom = nearest_points(boundary.centroid, dissolved_lines_gdf.unary_union)[1]
    if nearest_geom.is_empty:
        log("Error: Nearest point on trench network is empty.")
        continue
    # Connect to central node
    trench_G.add_edge(
        (nearest_geom.x, nearest_geom.y),
        (central_point_utm.x, central_point_utm.y),
        geometry=LineString([nearest_geom, (central_point_utm.x, central_point_utm.y)]),
        length=Point(nearest_geom).distance(Point(central_point_utm))
    )
log("Boundary boxes connected to the central node.")

# Snapping centroids to the nearest trench graph nodes in UTM
log("Snapping centroids to the nearest trench network nodes...")
def snap_to_nearest_trench_node(point, graph):
    nearest_node = None
    min_dist = float('inf')
    for node in graph.nodes:
        dist = point.distance(Point(node))
        if dist < min_dist:
            min_dist = dist
            nearest_node = node
    return nearest_node

centroid_nodes = []
for idx, centroid in centroids_gdf_utm.iterrows():
    snapped_node = snap_to_nearest_trench_node(centroid.geometry, trench_G)
    centroid_nodes.append(snapped_node)
    log(f"Centroid {idx + 1} snapped to node at {snapped_node}")

# Connect centroids to the trench network if not already connected
log("Connecting centroids to the trench network...")
for idx, (centroid, node) in enumerate(zip(centroids_gdf_utm.geometry, centroid_nodes), start=1):
    if node not in trench_G.nodes:
        trench_G.add_node(node)
        trench_G.add_edge(node, node, geometry=LineString([node, node]), length=0)
        log(f"Connected Centroid {idx} to trench network at node {node}")
    else:
        log(f"Centroid {idx} is already connected to trench network at node {node}")

# Define the central node in UTM
central_node = (central_point_utm.x, central_point_utm.y)
log(f"Central node defined as: {central_node}")

# Ensure central node is in graph
if central_node not in trench_G.nodes:
    trench_G.add_node(central_node)
    log("Central node added to trench network graph.")
else:
    log("Central node is already in trench network graph.")

# Find shortest paths from each centroid node to the central node
log("Calculating shortest paths from centroids to the central node...")
shortest_paths = {}
for idx, centroid_node in enumerate(centroid_nodes, start=1):
    if centroid_node == central_node:
        shortest_paths[centroid_node] = [central_node]
        log(f"Centroid {idx} is the central node.")
        continue
    try:
        path = nx.shortest_path(trench_G, source=centroid_node, target=central_node, weight='length')
        shortest_paths[centroid_node] = path
        log(f"Shortest path found for Centroid {idx}: {len(path)} nodes.")
    except nx.NetworkXNoPath:
        log(f"No path found from Centroid {idx} at node {centroid_node} to central node.")
    except Exception as e:
        log(f"Error finding path for Centroid {idx}: {e}")

# Define function to offset LineString to a preferred side
def offset_line(line, side, distance):
    try:
        if side not in ['left', 'right']:
            raise ValueError("SIDE_PREFERENCE must be 'left' or 'right'")
        offset = line.parallel_offset(distance, side=side, join_style=2, resolution=16)
        # Handle MultiLineString by merging
        if offset.type == 'MultiLineString':
            offset = unary_union(offset)
        return offset
    except Exception as e:
        log(f"Offsetting failed: {e}")
        return line  # Return original line if offset fails

# Function to create trench lines with side preference
def create_trench_lines(trench_graph, paths, side_pref, offset_dist, roads_gdf):
    trench_lines = []
    for path in paths.values():
        if len(path) < 2:
            continue
        path_coords = list(path)
        path_line = LineString(path_coords)
        # Offset the trench line to the preferred side
        offset_geom = offset_line(path_line, side_pref, offset_dist)
        # Check for intersections with roads
        intersects = roads_gdf.geometry.intersects(offset_geom)
        if intersects.any():
            log("Trench line intersects roads. Trying opposite side.")
            opposite_side = 'right' if side_pref == 'left' else 'left'
            offset_geom_opposite = offset_line(path_line, opposite_side, offset_dist)
            intersects_opposite = roads_gdf.geometry.intersects(offset_geom_opposite)
            if intersects_opposite.any():
                log("Trench line still intersects roads after opposite side offset. Using original path.")
                trench_lines.append(path_line)
            else:
                log("Trench line fixed by offsetting to the opposite side.")
                trench_lines.append(offset_geom_opposite)
        else:
            log("Trench line successfully offset to the preferred side.")
            trench_lines.append(offset_geom)
    return trench_lines

# Create trench lines with side preference
log("Creating trench lines with side preference...")
trench_lines = create_trench_lines(
    trench_G,
    shortest_paths,
    SIDE_PREFERENCE,
    OFFSET_DISTANCE,
    roads_edge_gdf
)
log(f"Created {len(trench_lines)} trench line(s).")

# Create a GeoDataFrame for trench lines in UTM
trench_gdf_utm = gpd.GeoDataFrame(geometry=trench_lines, crs=erven_gdf_utm.crs)

# Visualize the trench lines
visualize_geometries([*trench_lines, central_point_utm], title="Trench Lines with Central Point")

# Reproject trench lines back to lat/lon for KML output
trench_gdf_latlon = trench_gdf_utm.to_crs(epsg=4326)

# Create boundary boxes centroids for KML in UTM
boundary_centroids_utm = boundary_boxes_gdf_utm.geometry.apply(lambda box: box.centroid)
boundary_centroids_gdf_utm = gpd.GeoDataFrame(geometry=boundary_centroids_utm, crs=erven_gdf_utm.crs)

# Reproject boundary centroids to lat/lon for KML
boundary_centroids_latlon = boundary_centroids_gdf_utm.to_crs(epsg=4326)

# Reproject central point to lat/lon for KML
central_point_latlon = gpd.GeoSeries([central_point_utm], crs=erven_gdf_utm.crs).to_crs(epsg=4326).iloc[0]

# Reproject centroids to lat/lon for KML
centroids_gdf_latlon = centroids_gdf_utm.to_crs(epsg=4326)

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

document = kml.Document(ns='http://www.opengis.net/kml/2.2', name="Trench Network")
k.append(document)

# Add central node to KML
placemark_central = kml.Placemark(ns='http://www.opengis.net/kml/2.2')
placemark_central.name = "Central Node"
placemark_central.geometry = central_point_latlon
document.append(placemark_central)
log("Added Central Node to KML.")

# Add boundary boxes centroids to KML
for idx, centroid in boundary_centroids_latlon.iteritems():
    placemark_boundary = kml.Placemark(ns='http://www.opengis.net/kml/2.2')
    placemark_boundary.name = f'Boundary Box {idx + 1}'
    placemark_boundary.geometry = centroid
    document.append(placemark_boundary)
    log(f"Added Boundary Box {idx + 1} to KML.")

# Add trench lines to KML
for idx, trench in trench_gdf_latlon.iterrows():
    placemark_line = kml.Placemark(ns='http://www.opengis.net/kml/2.2')
    placemark_line.name = f"Trench Line {idx + 1}"
    placemark_line.geometry = trench.geometry
    document.append(placemark_line)
    log(f"Added Trench Line {idx + 1} to KML.")

# Add centroid points to KML
for idx, centroid in centroids_gdf_latlon.iterrows():
    placemark_centroid = kml.Placemark(ns='http://www.opengis.net/kml/2.2')
    placemark_centroid.name = f"Centroid {idx + 1}"
    placemark_centroid.geometry = centroid.geometry
    document.append(placemark_centroid)
    log(f"Added Centroid {idx + 1} to KML.")

# Save KML
output_kml_path = "output.kml"
try:
    with open(output_kml_path, "w") as f:
        f.write(k.to_string(prettyprint=True))
    log(f"KML output saved to '{output_kml_path}'")
except Exception as e:
    log(f"Error saving KML: {e}")
    sys.exit(1)