modify this script so that it interconnects the dissolved polygon edges by doing a nearest neighbour analysis and drawing a blue linestring between the seperate pieces. 
import geopandas as gpd
import osmnx as ox
import numpy as np
from shapely.geometry import LineString, Point
from shapely.ops import unary_union, nearest_points
from itertools import combinations
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 with separate colors for dissolved edges and connections
def visualize_geometries(dissolved_geometries, connection_geometries, title):
    plt.figure(figsize=(10, 10))
    
    # Plot dissolved polygon edges in red
    for geom in dissolved_geometries:
        if geom.type == 'LineString':
            x, y = geom.xy
            plt.plot(x, y, color='red', linewidth=2, label='Dissolved Edges' if 'Dissolved Edges' not in plt.gca().get_legend_handles_labels()[1] else "")
    
    # Plot connecting lines in blue
    for geom in connection_geometries:
        if geom.type == 'LineString':
            x, y = geom.xy
            plt.plot(x, y, color='blue', linestyle='--', linewidth=1, label='Connections' if 'Connections' not in plt.gca().get_legend_handles_labels()[1] else "")
    
    plt.title(title)
    plt.xlabel("Easting (m)")
    plt.ylabel("Northing (m)")
    plt.axis("equal")
    plt.grid()
    plt.legend()
    plt.show()

# Parameters
TARGET_UTM_EPSG = 32734  # EPSG code for WGS 84 / UTM zone 34S
input_layer_path = r"D:\Spatial\Piketberg Erven.shp"
CONNECTION_TOLERANCE = 10  # meters, adjust as needed

# 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)
    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 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)

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

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

# Function to extract endpoints from LineStrings
def get_endpoints(line):
    return [Point(line.coords[0]), Point(line.coords[-1])]

# Extract all endpoints
log("Extracting endpoints of dissolved lines...")
endpoints = []
for line in dissolved_lines:
    endpoints.extend(get_endpoints(line))

# Create a GeoDataFrame of endpoints
endpoints_gdf = gpd.GeoDataFrame(geometry=endpoints, crs=erven_gdf_utm.crs)

# Spatial index for efficient querying
endpoints_sindex = endpoints_gdf.sindex

# Initialize list to hold connection lines
connection_lines = []
connected = set()

log("Creating connections between nearby endpoints...")
# Iterate through each endpoint
for idx, point in endpoints_gdf.iterrows():
    if idx in connected:
        continue  # Skip already connected points

    possible_matches_index = list(endpoints_sindex.intersection(point.geometry.buffer(CONNECTION_TOLERANCE).bounds))
    possible_matches = endpoints_gdf.iloc[possible_matches_index]
    possible_matches = possible_matches[possible_matches.geometry.distance(point.geometry) <= CONNECTION_TOLERANCE]

    # Find the nearest point that is not the same point
    possible_matches = possible_matches[possible_matches.index != idx]
    if not possible_matches.empty:
        # Safely get the index of the nearest point
        nearest_idx = possible_matches.geometry.distance(point.geometry).idxmin()
        nearest_point = possible_matches.geometry.iloc[nearest_idx]  # Get the nearest point
        connection_line = LineString([point.geometry, nearest_point])
        connection_lines.append(connection_line)
        
        # Mark both points as connected
        connected.add(idx)
        connected.add(nearest_idx)  # Add nearest_idx

# Create a GeoDataFrame for connections
connections_gdf = gpd.GeoDataFrame(geometry=connection_lines, crs=erven_gdf_utm.crs)
log(f"Created {len(connections_gdf)} connection line(s).")