chatgpt_geograph.py

gistfile1.py

from shapely.geometry import LineString, Point
from shapely.ops import unary_union
from geopy.distance import geodesic

def snap_to_nearest_node(point, nodes, tolerance):
    for node in nodes:
        if geodesic((point.y, point.x), (node.y, node.x)).meters <= tolerance:
            return node
    return point

def create_nodes_edges(line_strings, tolerance=1.0):
    # Step 1: Merge all line strings into a single MultiLineString
    merged_lines = unary_union(line_strings)

    # Step 2: Extract nodes and edges with snapping using geodesic distances
    nodes = []
    edges = []
    
    for line in merged_lines:
        if isinstance(line, LineString):
            coords = list(line.coords)
            snapped_start = snap_to_nearest_node(Point(coords[0]), nodes, tolerance)
            snapped_end = snap_to_nearest_node(Point(coords[-1]), nodes, tolerance)

            if snapped_start not in nodes:
                nodes.append(snapped_start)
            if snapped_end not in nodes:
                nodes.append(snapped_end)

            edges.append((snapped_start, snapped_end))
        else:
            for segment in line:
                coords = list(segment.coords)
                snapped_start = snap_to_nearest_node(Point(coords[0]), nodes, tolerance)
                snapped_end = snap_to_nearest_node(Point(coords[-1]), nodes, tolerance)

                if snapped_start not in nodes:
                    nodes.append(snapped_start)
                if snapped_end not in nodes:
                    nodes.append(snapped_end)

                edges.append((snapped_start, snapped_end))

    return nodes, edges

# Example usage with lat/lon data
line1 = LineString([(-73.9857, 40.7484), (-73.9851, 40.7485)])
line2 = LineString([(-73.9851, 40.7485), (-73.9845, 40.7486)])
line3 = LineString([(-73.9845, 40.7486), (-73.9839, 40.7487)])

# The tolerance is set in meters
nodes, edges = create_nodes_edges([line1, line2, line3], tolerance=5.0)

print("Nodes:")
for node in nodes:
    print(node)

print("\nEdges:")
for edge in edges:
    print(edge)