fix, revise and improve: 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

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

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

# Parameters
SIDE_PREFERENCE = 'left'   # Options: 'left', 'right'
OFFSET_DISTANCE = 2        # Meters; adjust as needed based on CRS
TARGET_UTM_EPSG = 32734    # EPSG code for WGS 84 / UTM zone 34S
SAMPLING_DISTANCE = 100    # Meters; distance between boundary box points

# 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
# Note: User requested to remove buffering; hence, we skip buffering steps related to aggregation boxes

# Find shared boundaries between erven polygons
log("Identifying shared boundaries between erven polygons...")
# Extract boundaries of each polygon
boundaries = erven_gdf_utm.boundary

# Explode boundaries into individual LineStrings
log("Exploding boundaries into individual LineStrings...")
boundaries_exploded = boundaries.explode().reset_index(drop=True)

# Remove any non-LineString geometries
log("Filtering LineString geometries...")
boundaries_exploded = boundaries_exploded[boundaries_exploded.type == 'LineString']

# Create a GeoDataFrame of boundaries
log("Creating GeoDataFrame for boundaries...")
boundaries_gdf = gpd.GeoDataFrame(geometry=boundaries_exploded, crs=erven_gdf_utm.crs)

# Canonicalize LineStrings to have a consistent order for comparison
log("Canonicalizing LineStrings for comparison...")
boundaries_gdf['geometry'] = boundaries_gdf['geometry'].apply(lambda x: LineString(sorted(x.coords)))

# Count occurrences of each boundary LineString
log("Counting occurrences of each boundary LineString...")
boundary_counts = boundaries_gdf['geometry'].value_counts()

# Select boundaries that appear exactly twice (shared by two polygons)
shared_boundaries = boundary_counts[boundary_counts == 2].index.tolist()
log(f"Number of shared boundaries identified: {len(shared_boundaries)}")

# Create a GeoDataFrame for shared boundaries
shared_boundaries_gdf = gpd.GeoDataFrame(geometry=shared_boundaries, crs=erven_gdf_utm.crs)

# Create BB Points along shared boundaries
log("Creating BB Points along shared boundaries...")
bb_points = []
for boundary in shared_boundaries_gdf.geometry:
    # Place BB Point at the midpoint of the shared boundary
    midpoint = boundary.interpolate(0.5, normalized=True)
    # Calculate direction from central centroid to midpoint
    central_centroid = erven_gdf_utm.geometry.unary_union.centroid
    vector_x = midpoint.x - central_centroid.x
    vector_y = midpoint.y - central_centroid.y
    length = np.hypot(vector_x, vector_y)
    if length == 0:
        continue  # Skip if midpoint coincides with centroid
    unit_x = vector_x / length
    unit_y = vector_y / length
    # Offset the point outward by OFFSET_DISTANCE
    offset_x = midpoint.x + unit_x * OFFSET_DISTANCE
    offset_y = midpoint.y + unit_y * OFFSET_DISTANCE
    bb_point = Point(offset_x, offset_y)
    bb_points.append(bb_point)

log(f"Created {len(bb_points)} BB Points.")

# Create GeoDataFrame for BB Points
boundary_boxes_gdf_utm = gpd.GeoDataFrame(geometry=bb_points, crs=erven_gdf_utm.crs)

# Add BB Points to trench network
log("Integrating BB Points into trench network...")
for idx, bb_point in boundary_boxes_gdf_utm.iterrows():
    # Find the nearest node in the trench network graph
    try:
        nearest_node = ox.distance.nearest_nodes(
            roads_G, 
            X=bb_point.geometry.x, 
            Y=bb_point.geometry.y
        )
        # Convert nearest_node key back to coordinates in UTM
        if nearest_node in roads_node_gdf.index:
            nearest_node_coords = (roads_node_gdf.loc[nearest_node, 'x'], roads_node_gdf.loc[nearest_node, 'y'])
        else:
            log(f"Nearest node {nearest_node} not found in roads GeoDataFrame.")
            continue
        # Add the BB Point as a node
        bb_node = (bb_point.geometry.x, bb_point.geometry.y)
        trench_G.add_node(bb_node)
        # Connect BB Point to the nearest trench node
        trench_G.add_edge(
            bb_node,
            nearest_node_coords,
            geometry=LineString([bb_node, nearest_node_coords]),
            length=Point(bb_node).distance(Point(nearest_node_coords))
        )
        log(f"BB Point {idx + 1} connected to trench network at node {nearest_node_coords}.")
    except Exception as e:
        log(f"Error connecting BB Point {idx + 1}: {e}")
        continue

# Snapping centroids to the nearest trench network 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 = (erven_gdf_utm.geometry.unary_union.centroid.x, erven_gdf_utm.geometry.unary_union.centroid.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)
log(f"Trench lines GeoDataFrame created with {len(trench_gdf_utm)} lines.")

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

# Reproject BB Points to lat/lon for KML
boundary_boxes_latlon = boundary_boxes_gdf_utm.to_crs(epsg=4326)

# Reproject central point to lat/lon for KML
central_point_latlon = gpd.GeoSeries([Point(central_node)], 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_boxes_latlon.iterrows():
    placemark_boundary = kml.Placemark(ns='http://www.opengis.net/kml/2.2')
    placemark_boundary.name = f'Boundary Box {idx + 1}'
    placemark_boundary.geometry = centroid.geometry
    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)