create_isochrone_boundary.py

import numpy as np
import geopandas as gpd

# keep a geodataframe of all points creating the edges
all_points = gpd.GeoSeries(
    clipped_edges_gdf.get_coordinates(ignore_index=True).apply(
        lambda row: Point(row["x"], row["y"]), axis=1
    ),
    crs=4326,
)

# get only farthest points of clipped edges
end_points = clipped_edges_gdf.geometry.apply(
    lambda ls: Point(ls.coords[-1])
)

def get_angle(point1: Point, point2: Point):
    rads = np.arctan2(point2.y - point1.y, point2.x - point1.x)
    return np.degrees(rads)

# sort the points by the angle from the center
sorted_end_points = gpd.GeoSeries(
    sorted(end_points, key=lambda pt: get_angle(center_node_point, pt)), crs=4326
).drop_duplicates()
# create the first boundary using edges and sorted end points
isochrone_boundary = unary_union(
    polygonize(edges_union.union(LinearRing(sorted_end_points)).to_list()).geoms
)

# check if any point is outside the boundary
points_outside_boundary = all_points[~all_points.intersects(isochrone_boundary)]
if len(points_outside_boundary) > 0:
    # create new boundary from the endpoints and points outside
    end_points = gpd.pd.concat([end_points, points_outside_boundary])
    sorted_end_points = gpd.GeoSeries(
        sorted(end_points, key=lambda pt: get_angle(center_node_point, pt)),
        crs=4326,
    ).drop_duplicates()
    isochrone_boundary = unary_union(
        polygonize(edges_union.union(LinearRing(sorted_end_points)).to_list()).geoms
    )

# if there is a previous boundary - union it to avoid edges intersections
if previous_boundary is not None:
    isochrone_boundary = isochrone_boundary.union(previous_boundary)

# if the boundary is a multipolygon - keep only the biggest component
if isochrone_boundary.geom_type == "MultiPolygon":
    isochrone_boundary = sorted(
        isochrone_boundary.geoms, key=lambda g: g.area, reverse=True
    )[0]