Basic trip over a real-world graph

Basic trip over a real-world graph#

In this notebook, we set up a basic simulation where one vessel moves over a real-world graph. The vessel will sail from a location in the Maasvlakte, to a location in the Port of Moerdijk.

0. Import libraries#

import pathlib

# package(s) used for creating and geo-locating the graph
import networkx as nx
import shapely

# package(s) related to the simulation (creating the vessel, running the simulation)
import datetime
import simpy
import opentnsim
import opentnsim.fis as fis
from opentnsim.core.logutils import logbook2eventtable
from opentnsim.core.plotutils import generate_vessel_gantt_chart

# package(s) needed for inspecting the output
import pandas as pd
import geopandas as gpd
import movingpandas as mpd

# plot libraries
import folium
import pyarrow as pa
from lonboard import Map, PathLayer, viz
from lonboard.colormap import apply_categorical_cmap
from lonboard.experimental import TripsLayer

print("This notebook is executed with OpenTNSim version {}".format(opentnsim.__version__))

This notebook is executed with OpenTNSim version 1.3.7

1. Define object classes#

# make your preferred Vessel class out of available mix-ins.
Vessel = type(
    "Vessel", 
    (
        opentnsim.core.Identifiable, # allows to give the object a name and a random ID,
        opentnsim.core.Movable,      # allows the object to move, with a fixed speed, while logging this activity
    ), 
    {}
)

2. Create graph#

Next we create a network (a graph) along which the vessel can move. For this case we use the Fairway Information System graph, and make the vessel sail from a container terminal at the Maasvlakte to an container terminal at the Port of Moerdijk.

# load the processed version from the Fairway Information System graph provided by Rijkswaterstaat
FG = fis.load_network(version="0.3")

maasvlakte = shapely.Point(4.0566946, 51.9471624)
moerdijk = shapely.Point(4.5944738, 51.6829037)

nodes_gdf = gpd.GeoDataFrame(FG.nodes.values(), index=FG.nodes.keys())

distances, idx = nodes_gdf.sindex.nearest(maasvlakte)
maasvlakte_node = nodes_gdf.iloc[idx[0]]

distances, idx = nodes_gdf.sindex.nearest(moerdijk)
moerdijk_node = nodes_gdf.iloc[idx[0]]

route = nx.shortest_path(FG, maasvlakte_node.name, moerdijk_node.name, weight='length_m')

# Create a map centered between the two points
m = folium.Map(location=[51.83, 4.33], zoom_start = 10)

for edge in FG.edges(data = True):
    points_x = list(edge[2]["geometry"].coords.xy[0])
    points_y = list(edge[2]["geometry"].coords.xy[1])
    
    line = []
    for i, _ in enumerate(points_x):
        line.append((points_y[i], points_x[i]))

    if edge[0] in route and edge[1] in route:
        folium.PolyLine(line, color = "black", weight = 3, popup = edge[2]["Name"]).add_to(m)
    else:
        folium.PolyLine(line, color = "blue", weight = 1, popup = edge[2]["Name"]).add_to(m)


# Add round marker for Maasvlakte with popup
folium.CircleMarker(
    location=[maasvlakte_node.Y, maasvlakte_node.X],
    radius=8,
    color='green',
    fill=True,
    fill_color='green',
    fill_opacity=0.7,
    popup='Maasvlakte'
).add_to(m)

# Add round marker for Moerdijk with popup
folium.CircleMarker(
    location=[moerdijk_node.Y, moerdijk_node.X],
    radius=8,
    color='blue',
    fill=True,
    fill_color='blue',
    fill_opacity=0.7,
    popup='Moerdijk'
).add_to(m)

# Display the map
m
Make this Notebook Trusted to load map: File -> Trust Notebook
FG.nodes['8864054']

{'n': '8864054',
 'X': 3.56379034355024,
 'Y': 51.7453141508194,
 'geometry': <POINT (3.564 51.745)>,
 'Wkt': 'POINT (3.5637903435502398 51.7453141508193966)'}

FG.edges[('8861095', '8864054')]

{'GeoType': 'section',
 'Name': 'Vaarwegvak van 0 tot 2 - H',
 'Length': 2.346,
 'GeneralDepth': nan,
 'GeneralLength': nan,
 'GeneralWidth': nan,
 'SeaFairingDepth': nan,
 'PushedLength': nan,
 'PushedWidth': nan,
 'GeneralHeight': nan,
 'SeaFairingLength': nan,
 'SeaFairingWidth': nan,
 'CoupledLength': nan,
 'CoupledWidth': nan,
 'PushedDepth': nan,
 'WidePushedDepth': nan,
 'CoupledDepth': nan,
 'WidePushedLength': nan,
 'WidePushedWidth': nan,
 'SeaFairingHeight': nan,
 'Id_navigability': nan,
 'Classification': nan,
 'Code': nan,
 'Description': nan,
 'length_deg': nan,
 'length': 0.0255813258548747,
 'Wkt': 'LINESTRING (3.54535894046351 51.727661900382, 3.5602008295784 51.742791566037, 3.56379034355024 51.7453141508194)',
 'StartJunctionId': '8861095',
 'EndJunctionId': '8864054',
 'subgraph': 0,
 'length_m': 2835.225880186762,
 'geometry': <LINESTRING (3.545 51.728, 3.56 51.743, 3.564 51.745)>}

3. Run simulation#

def mission(env, vessel):
    """
    Method that defines the mission of the vessel. 
    
    In this case: 
        keep moving along the path until its end point is reached
    """
    while True:
        yield from vessel.move()
        
        if vessel.geometry == nx.get_node_attributes(env.graph, "geometry")[vessel.route[-1]]:
            break

# start simpy environment
simulation_start = datetime.datetime(2024, 1, 1, 0, 0, 0)
env = simpy.Environment(initial_time=simulation_start.timestamp())
env.epoch = simulation_start

# add graph to environment
env.graph = FG

# create vessel from a dict 
path = route
data_vessel = {
    "env": env,                                       # needed for simpy simulation
    "name": "Vessel",                                 # required by Identifiable
    "geometry": env.graph.nodes[path[0]]['geometry'], # required by Locatable
    "route": path,                                    # required by Routeable
    "v": 1,                                           # required by Movable, 1 m/s to check if the distance is covered in the expected time
}  # 

# create an instance of the Vessel class using the input dict data_vessel
vessel = Vessel(**data_vessel)

# start the simulation
env.process(mission(env, vessel))
env.run()

4. Inspect output#

We can now inspect the simulation output by inspecting the vessel.logbook. Note that the Log mix-in was included when we added Movable. The vessel.logbook keeps track of the moving activities of the vessel. For each discrete event OpenTNSim logs an event message, the start/stop time and the location. The vessel.logbook is of type dict. For convenient inspection it can be loaded into a Pandas dataframe.

# load the logbook data into a dataframe
df = pd.DataFrame.from_dict(vessel.logbook)

print("'{}' logbook data:".format(vessel.name))  
print('')

display(df)

'Vessel' logbook data:
Message Timestamp Value Geometry
0 Sailing from node 8865973 to node 8867633 start 2024-01-01 00:00:00.000000 0.000000 POINT (4.03967999329435 51.9453492703995)
1 Sailing from node 8865973 to node 8867633 stop 2024-01-01 00:44:24.379610 2664.379609 POINT (4.07753550656458 51.9504653911648)
2 Sailing from node 8867633 to node 8862102 start 2024-01-01 00:44:24.379610 2664.379609 POINT (4.07753550656458 51.9504653911648)
3 Sailing from node 8867633 to node 8862102 stop 2024-01-01 00:49:28.777612 2968.777612 POINT (4.08193446085462 51.95077489336861)
4 Sailing from node 8862102 to node 8861217 start 2024-01-01 00:49:28.777612 2968.777612 POINT (4.08193446085462 51.95077489336861)
... ... ... ... ...
83 Sailing from node 8860742 to node 30986654 stop 2024-01-01 16:27:29.435584 59249.435585 POINT (4.58885224011778 51.6948284745577)
84 Sailing from node 30986654 to node 8865570 start 2024-01-01 16:27:29.435584 59249.435585 POINT (4.58885224011778 51.6948284745577)
85 Sailing from node 30986654 to node 8865570 stop 2024-01-01 16:37:38.270048 59858.270049 POINT (4.59124447578966 51.6895622155408)
86 Sailing from node 8865570 to node 8864978 start 2024-01-01 16:37:38.270048 59858.270049 POINT (4.59124447578966 51.6895622155408)
87 Sailing from node 8865570 to node 8864978 stop 2024-01-01 16:57:15.383138 61035.383139 POINT (4.5959831688951 51.6794008186706)

88 rows × 4 columns

The inspection of the logbook data shows that Vessel moved from its origin (Node 0) to its destination (Node 1). The print statements show that the length of the route from Node 0 to Node 1 is exactly 100 km. At a given speed of 1 m/s the trip duration should be exactly 100000 seconds, as is indeed shown to be the case.

df_eventtable = logbook2eventtable([vessel])
df_eventtable.head(3)
object id object name activity name start location stop location start time stop time distance (m) duration (s)
0 8202da4c-cc87-47a8-9019-45eb85f570a9 Vessel Sailing from node 8865973 to node 8867633 POINT (4.03967999329435 51.9453492703995) POINT (4.07753550656458 51.9504653911648) 2024-01-01 00:00:00.000000 2024-01-01 00:44:24.379610 2664.379609 2664.379610
1 8202da4c-cc87-47a8-9019-45eb85f570a9 Vessel Sailing from node 8867633 to node 8862102 POINT (4.07753550656458 51.9504653911648) POINT (4.08193446085462 51.95077489336861) 2024-01-01 00:44:24.379610 2024-01-01 00:49:28.777612 304.398002 304.398002
2 8202da4c-cc87-47a8-9019-45eb85f570a9 Vessel Sailing from node 8862102 to node 8861217 POINT (4.08193446085462 51.95077489336861) POINT (4.07844925662863 51.9415581666513) 2024-01-01 00:49:28.777612 2024-01-01 01:07:01.917087 1053.139475 1053.139475 
generate_vessel_gantt_chart(df_eventtable)

duration = df_eventtable['duration (s)'].sum()
distance = df_eventtable['distance (m)'].sum()

print('{} sailed {:.1f} m in {:.1f} s, which is {}'.format(vessel.name, distance, duration, str(datetime.timedelta(seconds=duration))))

Vessel sailed 61035.4 m in 61035.4 s, which is 16:57:15.383138

5. Visualize simulation#

Animation is only of points in the log, resulting in the ship sailing over land.

# Create a geopandas GeoDataFrame from the vessel_log
vessel_log_gdf = gpd.GeoDataFrame(df, geometry=df['Geometry'], crs='EPSG:4326')

# Drop the 'Geometry' column, since it is now a duplicate of the 'geometry' column created by geopandas
vessel_log_gdf = vessel_log_gdf.drop(columns=['Geometry'])

src_dir = pathlib.Path(opentnsim.__file__).parent.parent

# Graph location
location_graph = src_dir / "notebooks"

# Create a Trajectory using movingpandas with the id being the ship name.
traj = mpd.Trajectory(vessel_log_gdf, t='Timestamp', traj_id='ship_1')
# Convert Trajectory to GeoDataFrame to save to a .gpkg file
traj_gdf = traj.to_point_gdf()
traj_gdf.to_file(location_graph / '0006 - route.gpkg')

INFO:pyogrio._io:Created 45 records

# For visualization, the typing of the attributes needs to be consistent. Here we choose to convert them all to string.
traj_gdf['Value'] = traj_gdf['Value'].astype(str)
viz(traj_gdf)

# For the animation we need to create a TrajectoryCollection first. Here we construct it again from the point_gdf
traj_collection = mpd.TrajectoryCollection(traj_gdf, traj_id_col='traj_id', t='t')

# Now we only have one ship, but you can assign different colors to different ships
trajid_to_color = {
    "ship_1": "blue",
}

get_color = apply_categorical_cmap(pa.array(traj_gdf['traj_id'].unique()), trajid_to_color)

# Create the visualization from the TrajectoryCollection
trips_layer = TripsLayer.from_movingpandas(
    traj_collection,
    width_min_pixels=5,
    trail_length=1000,
    get_color=get_color
)

/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/lonboard/experimental/traits.py:151: UserWarning:

Reducing precision of input timestamp data to 's' to fit into available GPU precision.

# Show the visualization on the map
m = Map(trips_layer, _height=600)
m

# Animate the visualization
trips_layer.animate(step=datetime.timedelta(seconds=30), fps=30)
Contents