Is there a way to approximate a polygon to a polygon in a road network

Question

I'm working with polygons in maps. I have a df, df_polygons, each row in the df refers to a different polygon. And df_polygons['polygon'] retrieves a specific polygon, defined by a set of ordered points, an example of 2 polygons that are in the dataframe:

polygon_1 = [[41.2169494797, -8.6778177263],
 [41.2169902336, -8.692729496],
 [41.2282499337, -8.6926767872],
 [41.2282887721, -8.7075911445],
 [41.2339186214, -8.7075660584],
 [41.2451783036, -8.7075158648],
 [41.2451394423, -8.6925976679],
 [41.2507692675, -8.6925712798],
 [41.250749108, -8.6851115523],
 [41.2563789237, -8.6850845162],
 [41.2563582769, -8.6776241559],
 [41.2394688256, -8.6777071617],
 [41.2393814564, -8.6478734838],
 [41.2281218289, -8.6479339004],
 [41.2282091637, -8.6777624597],
 [41.2169494797, -8.6778177263]]

polygon_2 = [[41.194470767, -8.6928348239],
 [41.2057305114, -8.6927821749],
 [41.2056897735, -8.6778729616],
 [41.2169494797, -8.6778177263],
 [41.2169902336, -8.692729496],
 [41.2282499337, -8.6926767872],
 [41.2282887721, -8.7075911445],
 [41.2508081364, -8.7074907573],
 [41.2508800748, -8.737329796200001],
 [41.1945813563, -8.7375549721],
 [41.194470767, -8.6928348239]]

Plotting the two polygons, I get:

But, I'm wondering if it is possible to get an approximated polygon, where every edge of the polygon is in the road network. Until now, I have come up with this:

m = folium.Map(location=(41.25387369401857, -8.781923733405996), zoom_start=12)

def interpolate_points(coord1, coord2, n):
    lat1, lon1 = coord1
    lat2, lon2 = coord2

    lats = np.linspace(lat1, lat2, n)
    lons = np.linspace(lon1, lon2, n)

    points = [(lat, lon) for lat, lon in zip(lats, lons)]
    return points

df_polygons = pd.read_json()
df_polygons = df_polygons.loc[0:2]

for p in range(len(df_polygons)):
    coords = df_polygons.loc[p, 'polygon']

    north = max([lat for lat, lon in coords])
    south = min([lat for lat, lon in coords])
    east = max([lon for lat, lon in coords])
    west = min([lon for lat, lon in coords])
    G = ox.graph.graph_from_bbox(north, south, east, west, network_type='drive')
    G_edges = ox.graph_to_gdfs(G,nodes=False)
        
    edges = [[coords[i], coords[i+1]] for i in range(len(coords)-1)]
        
    all_data = []
    for i in range(len(edges)):
        coord1 = edges[i][0]  # (lat1, lon1)
        coord2 = edges[i][1] # (lat2, lon2)

        n = 10  # Number of points to interpolate
        interpolated_points = interpolate_points(coord1, coord2, n)
        nearest_edges = ox.distance.nearest_edges(G, [x for y, x in interpolated_points], [y for y, x in interpolated_points], interpolate=None, return_dist=False)
            
        last_point = None
        for j in range(len(nearest_edges)):
            selected_row = G_edges.loc[(nearest_edges[j][0], nearest_edges[j][1])]
            road_data = list(selected_row.geometry.iloc[0].coords)
                
            all_data.append(road_data)
                
    road_data_1d = [item for sublist in all_data for item in sublist]

    folium.PolyLine([(road_data_1d[k][1], road_data_1d[k][0]) for k in range(len(road_data_1d))], color='blue', fill_color='blue').add_to(m)
    
m.save("test-figure.html")
m

First, I compute the graph G as the bounding box of the polygon. Then, basically, I interpolate every edge of the polygon, and I get n points belonging to the edge. Then, I compute the nearest edge (on G) to the point. The result of this (for the same 2 polygons) is the following:

There are several gaps between the two polygons, I was wondering if there is a better solution that mine, that could retrieve a more accurate approximation.

Thanks in advance!