python 3.x - OSMNX + NetworkX graph_from_gdfs incorrect nodes returned using nx.ego_graph

Question

Attempting to load a geopackage edge and node layers as NetworkX Graph but failing to get the correct nodes returned when using 'nx.ego_graph'. Using version OSMNX 1.0, NetworkX 2.5.

In a previous question G Boeing assisted in using the Geopackage to Graph function.

I now have a Graph object but the analysis doesn't appear to work as it does when using the True Graph object e.g. G = ox.graph_from_point((53.956748, -1.081676)).

I want to create isochrones for different travel times and reuse the clever work Mr Boeing has already done in this example. I would like to use the 'buffer' option predominantly- rather than the convex hull.

The analysis from the example above currently works OK when using G = ox.graph_from_point((53.956748, -1.081676),dist=2000, dist_type='network') i.e. dynamically downloading the osm data but I have a specific workflow where I export data to a geopackage and edit, add some nodes, and edges to represent possible new roads.

Even without editing the Geopackage data, I cannot get the analysis to work, the nodes returned for different travel times appear wrong. The isochrone polygon for the 10 min trip time is the same as the 5 min trip time. And I don't believe the 5 min trip time polygon is correct as it's a very small area.

If you run the script below, you can already see when setting the node colours there is an issue as trip time 5 and 10 return the same exact nodes (subgraph = nx.ego_graph(G=G, n=origin_node, radius=trip_time, distance='time')).

Does anyone know if this indicates that perhaps the Graph object hasn't been generated properly when using 'ox.graph_from_gdfs'? Or is there a change that is required?

I tried using 'undirected=True' in nx.ego_graph and that returned...possibly the correct nodes in the colour setting section but when generating the polygons, I got this error:

edge_lookup = G.get_edge_data(n_fr, n_to)[0].get('geometry',  LineString([f,t]))
KeyError: 0

In the testing I did I found:

• using the buffer isochrone method along with nx.ego_graph: undirected=True - gives the error above
• using the buffer isochrone method along with nx.ego_graph: undirected=False - gives the most bizarre pear shaped polygon which covers half of Northern Europe (should be a small site in the UK < 2000m wide)- same for all trips
• using the convex hull isochrone method along with nx.ego_graph: undirected=False- gives identical polygons for each trip, and neither polygon is the correct shape - far too small
• using the convex hull isochrone method along with nx.ego_graph: undirected=True- gives the correct polygons for each trip, 10 mins bigger than 5 mins, etc

See example again for the convex hull method.

Hopefully someone has an idea how to correct the error.

Code:

geopPath = "Transport_Analysis.gpkg"
geopPath2 = "Transport_Analysis2.gpkg"

yG = 53.9582559586775  # EPSG:4326 Geographic
xG = -1.07662649048548 # EPSG:4326 Geographic
yP = 451746.474534685  # EPSG:27700 Projected
xP = 460685.875721401  # EPSG:27700 Projected

distance = 2000
network_type = 'walk'
trip_times = [5, 10, 15, 20, 25, 30]
travel_speed = 4.5

trip_times = trip_times[0:2]  # limit to just first 2 for testing

useWGS84 = True # True False Choose to use Geographic or Projected Coordinates

if useWGS84 == True:
    epsg = 4326
    orig = (yG,xG) # EPSG:4326
    print("Using Geographic coords: {},{}".format(xG,yG))
    graph_attrs = {'crs': 'epsg:4326', 'simplified': True}
else:
    epsg = 27700 # for my project site in the UK
    orig = (yP,xP) # EPSG:27700
    print("Using Projected coords: {},{}".format(xP,yP))
    graph_attrs = {'crs': 'epsg:{}'.format(epsg), 'simplified': True}

download = 0
if download == 1:
    print("Downloading OSMNX")
    G = ox.graph_from_point(orig, dist=distance, dist_type='network')
    origin_node = ox.get_nearest_node(G, orig)
    print("Downloaded OSMNX")                       
    ox.io.save_graph_geopackage(G, filepath=geopPath2, encoding='utf-8')
    geopPath = geopPath2
else:
    print("Using existing geopackage: {}".format(geopPath))

# geopPath edges and nodes in native osm EPSG:4326
# load GeoPackage as node/edge GeoDataFrames indexed as described in OSMnx docs
gdf_nodes = gpd.read_file(geopPath, layer='nodes').set_index('osmid')
gdf_edges = gpd.read_file(geopPath, layer='edges').set_index(['u', 'v', 'key'])                
assert gdf_nodes.index.is_unique and gdf_edges.index.is_unique

print("### Loading Existing geopackage road edges and nodes as Graph")
# convert the node/edge GeoDataFrames to a MultiDiGraph
G = ox.graph_from_gdfs(gdf_nodes, gdf_edges, graph_attrs)

if useWGS84 != True:  # choose to use projected coordinates for origin/graph instead of geographic
    crsObj = CRS.from_user_input(epsg)
    # Finally converting graph to projected CRS
    G = ox.project_graph(G,to_crs=crsObj)  

origin_node = ox.get_nearest_node(G, orig)
print("closest node to coordinate: {} (coord Y,X: {})".format(origin_node,orig))

# add an edge attribute for time in minutes required to traverse each edge
meters_per_minute = travel_speed * 1000 / 60 #km per hour to m per minute
for u, v, k, data in G.edges(data=True, keys=True):
    #print(data['time'])
    data['time'] = data['length'] / meters_per_minute
iso_colors = ox.plot.get_colors(n=len(trip_times), cmap='plasma', start=0, return_hex=True)

# color the nodes according to isochrone then plot the street network
node_colors = {}
for trip_time, color in zip(sorted(trip_times, reverse=False), iso_colors):
    subgraph = nx.ego_graph(G=G, n=origin_node, radius=trip_time,undirected=True, distance='time') # 
    print("Trip time: {}. Subgraph.nodes: {}. Subnodes: {}".format(trip_time,len(subgraph.nodes()),subgraph.nodes()))
    for node in subgraph.nodes():
        node_colors[node] = color

tripReverse = False

edge_buff=25
node_buff=0
infill=False
isochrone_polys = []

#trip_times = list(trip_times)

for trip_time in sorted(trip_times, reverse=tripReverse):  # ,undirected=True
    subgraph = nx.ego_graph(G=G, n=origin_node, radius=trip_time,undirected=True, distance='time')

    node_points = [Point((data['x'], data['y'])) for node, data in subgraph.nodes(data=True)]
    nodes_gdf = gpd.GeoDataFrame({'id': list(subgraph.nodes())}, geometry=node_points)
    nodes_gdf = nodes_gdf.set_index('id')
    
    print("Trip time: {}. Nodes: {}".format(trip_time,len(node_points)))
    #print("node_points: {}".format(node_points))
    edge_lines = []
    for n_fr, n_to in subgraph.edges():
        f = nodes_gdf.loc[n_fr].geometry
        t = nodes_gdf.loc[n_to].geometry
        edge_lookup = G.get_edge_data(n_fr, n_to)[0].get('geometry',  LineString([f,t]))
        edge_lines.append(edge_lookup)

    n = nodes_gdf.buffer(node_buff).geometry
    e = gpd.GeoSeries(edge_lines).buffer(edge_buff).geometry
    all_gs = list(n) + list(e)
    new_iso = gpd.GeoSeries(all_gs).unary_union
    
    # try to fill in surrounded areas so shapes will appear solid and blocks without white space inside them
    if infill:
        new_iso = Polygon(new_iso.exterior)
    isochrone_polys.append(new_iso)
    #print("new_iso: {}".format(new_iso))

        

iso_df = pd.DataFrame({'trip_times_mins':trip_times[::-1],'color':colors,'geom':polys[::-1]})
iso_df = gpd.GeoDataFrame(iso_df,crs='epsg:{}'.format(epsg),geometry='geom')
# export isochrones to geopackage 
iso_df.to_file(geopPath, driver='GPKG',layer='isochrone_{}'.format(network_type))
print("Written isochrone")

Answer 1

等待大神答复