How can I find all the points that are a exactly certain distance, say 25km and 50km, from the start node?

Question

I was writing a Python script that asks input coordinates from the user and uses it to get the nearest node as well as finding all the points that are 25km and 50km away from the nearest node, along the network. The network infrastructure here is railway.

In brief, how I went about this is that from the input coordinates, I obtained the nearest node using this method ox.distance.nearest_nodes, then I used the coordinates of the nearest node as the source / origin node. While the target(s) / destination nodes were all the other nodes in the graph. I then used nx.shortest_path_length method to get distances from the source node to every other node on the graph. Then I used an if-condition to determine those points that are likely to be at exactly 25km and 50km respectively from the origin node. I created a range for this because it's almost unlikely to find a node at exactly 25km or 50km.

This may have been the best hack from my arsenal. However, I'm finding some issues. First is that I don't get the point at exactly 25km or 50km but rather nodes close to that, because of my filter range. So what happens is that if I'm using a range of let's say 24800m and 25200m, I'll find several number of nodes within this range. They could be along the same path or even adjacent to each other.

Here is my code:

# import the required libraries
import networkx as nx
import osmnx as ox
import requests
import matplotlib.cm as cm
import matplotlib.colors as colors
import folium
import warnings
%matplotlib inline
ox.config(log_console=True)
ox.__version__

# filter warnings
warnings.simplefilter(action="ignore", category=FutureWarning)
warnings.simplefilter(action="ignore", category=PendingDeprecationWarning)
warnings.simplefilter(action="ignore", category=UserWarning)

### Getting the railway infrastructure
# Let's ask input from user - coordinates
#latitude = float(input("Enter latitude: "))
#longitude = float(input("Enter longitude: "))

# graph from input coordinates
# default coordinates for testing are 51.5 (lat) and -0.13 (lon) - London
latitude = 51.5
longitude = -0.13
G = ox.graph_from_point(
    (latitude, longitude),
    retain_all=False,
    dist=70000, # I believe this also has an impact on speed. Since my ceiling target it 50km, can I put it at 50000?
    dist_type='network',
    truncate_by_edge=True,
    simplify=True,
    custom_filter='["railway"~"way"]',
)

# plotting the graph
fig, ax = ox.plot_graph(G, node_size=0, edge_color="w", edge_linewidth=0.2)

# finding the node that is closest to the input coordinates
gdf_nodes = ox.graph_to_gdfs(G, edges=False)
closest_node = ox.distance.nearest_nodes(G, longitude, latitude)

# getting the coordinates of the closest node
closest_node_lat = gdf_nodes.loc[[closest_node]].y.values[0]
closest_node_lon = gdf_nodes.loc[[closest_node]].x.values[0]

# calculating the length of shortest distances
start_node_coords = (latitude, longitude)
orig_node = ox.get_nearest_node(G, (latitude, longitude))
edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
z = 1
length_list = []

# let's create a color code for points at x distance
green = []
blue = []
red = []
for coordinate in gdf_nodes['geometry']:
    lat = coordinate.y
    lon = coordinate.x
    dest_node = ox.get_nearest_node(G, (lat, lon))
    try:
        length = nx.shortest_path_length(G, source=orig_node, target=dest_node, weight='length')
        print(length, coordinate, z)
        length_list.append(length)
        # let's create a filter that finds if the length is either
        # 10km, 25km or 50km. We're going to find these points from 
        # our origin node, along all the available paths
        if length > 9800 and length < 10200:
            ten_km_coordinates = (lat, lon)
            green.append(ten_km_coordinates)
        elif length > 24800 and length < 25200:
            twenty_five_km_coordinates = (lat, lon)
            blue.append(twenty_five_km_coordinates)
        elif length > 49900 and length < 50100:
            fifty_km_coordinates = (lat, lon)
            red.append(fifty_km_coordinates)
            
        z = z + 1
    except nx.exception.NetworkXNoPath:
        G.remove_nodes_from([dest_node])
        print(coordinate)

# creating a map element
map = folium.Map(location=[closest_node_lat, closest_node_lon], tiles="Stamen Terrain", zoom_start=12)
folium.Marker(location=[latitude, longitude], popup="My Input", tooltip="My Input").add_to(map)
folium.Marker(location=[closest_node_lat, closest_node_lon], popup="Closest Node", tooltip="Closest Node", icon=folium.Icon(color='green')).add_to(map)
for i in range(0, len(green)):
    folium.Marker([green[i][0], green[i][1]], popup="10KM Away", icon=folium.Icon(color='green')).add_to(map)
for i in range(0, len(blue)):
    folium.Marker([blue[i][0], blue[i][1]], popup="25KM Away", icon=folium.Icon(color='blue')).add_to(map)
for i in range(0, len(red)):
    folium.Marker([red[i][0], red[i][1]], popup="50KM Away", icon=folium.Icon(color='red')).add_to(map)

folium.Choropleth(edges.geometry, line_weight=2, line_color='yellow').add_to(map)

Is there a better way I could do this to achieve my objective or atleast improve this code to find the points at exactly 25km and 50km respectively and not the nodes, because for me I'm just finding the nodes? I want to find all the those points at every available path / route from the origin node, not just one path. More like saying, if I travel along this route, where will I be at 25km? What about along this other route? I also realized the speed of this depends on the number of nodes. The more the nodes, the more it will take the script to generate results. I will greatly appreciate any help. Thanks!