Strip the code from this script that generates the lines from the dissolved polygon edges so the script only creates these: 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)