Visualizing the sailed path#
Imports#
Import the required libraries
# package(s) related to time, space and id
import datetime
import math
import os
import pathlib
import platform
import random
import time
import warnings
import geopandas as gpd
import matplotlib.animation as animation
import matplotlib.pyplot as plt
# Used for making the graph to visualize our problem
import networkx as nx
import numpy as np
# OpenTNSIM
import opentnsim
import opentnsim.core as core
import opentnsim.utils
# package(s) for data handling
import pandas as pd
# spatial libraries
import pyproj
import shapely.geometry
# you need these dependencies (you can get these from anaconda)
# package(s) related to the simulation
import simpy
from shapely.errors import ShapelyDeprecationWarning
from simplekml import Kml, Style
warnings.filterwarnings("ignore", category=ShapelyDeprecationWarning)
print("This notebook has been tested with OpenTNSim version {}".format(opentnsim.__version__))
This notebook has been tested with OpenTNSim version 1.3.7
# Graph location
src_dir = pathlib.Path(opentnsim.__file__).parent.parent
# Graph location
location_graph = src_dir / "notebooks"
name_graph = location_graph / "Shape-Files" / "Rotterdam-Antwerpen-corridor" / "edges_2.shp"
Create graph#
Important:
If you use windows and get the following error “ImportError: read_shp requires OGR: http://www.gdal.org/”, you probably have this issue. Solving it is possible by running the following commands in your terminal (as explained here):
#Create a new virtual environment
conda create -n testgdal -c conda-forge gdal vs2015_runtime=14
#Activate virtual environment
activate testgdal
#Open Jupyter notebook
jupyer notebook
def read_shp(name_graph):
gdf = gpd.read_file(name_graph)
def source_target(geom):
source, *middle, target = geom.coords
# source = shapely.Point(source).wkt
# target = shapely.Point(target).wkt
return pd.Series({"source": source, "target": target})
nodes = gdf["geometry"].apply(source_target)
gdf["source"] = nodes["source"]
gdf["target"] = nodes["target"]
FG = nx.from_pandas_edgelist(gdf)
return FG
# Read the shape-file
FG = read_shp(name_graph)
# Draw the shape-file to get a first impression
plt.figure(figsize=(18, 18))
nx.draw(FG)
# Show the drawing
plt.show()
# calculate distance between two points
def calculate_distance(orig, dest):
wgs84 = pyproj.Geod(ellps="WGS84")
distance = wgs84.inv(
shapely.geometry.shape(orig).x,
shapely.geometry.shape(orig).y,
shapely.geometry.shape(dest).x,
shapely.geometry.shape(dest).y,
)[2]
return distance
H_G = nx.Graph()
lat_lon_to_index = {}
edge_id_counter = 0
for i, node in enumerate(FG.nodes(data=True)):
H_G.add_node(i, pos=node[0], name="Node {}".format(i), geometry=shapely.geometry.Point(node[0][0], node[0][1]))
lat_lon_to_index[node[0]] = i
for edge in FG.edges(data=True):
H_G.add_edge(
lat_lon_to_index[edge[0]],
lat_lon_to_index[edge[1]],
dis=calculate_distance(
nx.get_node_attributes(H_G, "geometry")[lat_lon_to_index[edge[1]]],
nx.get_node_attributes(H_G, "geometry")[lat_lon_to_index[edge[0]]],
),
)
FG = H_G.to_directed()
plt.figure(figsize=(18, 18))
nx.draw(FG, nx.get_node_attributes(FG, "pos"), with_labels=True, node_size=0.5, font_size=2, width=0.2, arrowsize=3)
plt.show()
Create vessels#
Vessel without graph, but with shortest path.
# Make a class out of mix-ins
TransportResource = type(
"TransportResource", (core.Identifiable, core.Movable, core.HasResource, core.Routable, core.ExtraMetadata), {}
)
# For testing purposes we only want v to be 1
def compute_v_provider(v_empty, v_full):
return lambda x: 1
data_vessel = {
"env": None,
"name": "Vessel number 1",
"route": None,
"geometry": shapely.geometry.Point(0, 0), # lon, lat
"capacity": 1_000,
"v": 1,
}
# create the transport processing resource
vessel = TransportResource(**data_vessel)
Define paths#
# First simulation is from random node 1 to random node 2
source = 80
target = 287
path = nx.dijkstra_path(FG, source, target)
Run simulation#
def start(env, vessel):
while True:
vessel.log_entry_v0("Start sailing", env.now, "", vessel.geometry)
yield from vessel.move()
vessel.log_entry_v0("Stop sailing", env.now, "", vessel.geometry)
if vessel.geometry == nx.get_node_attributes(FG, "geometry")[vessel.route[-1]]:
break
# Start simpy environment
simulation_start = datetime.datetime.now()
env = simpy.Environment(initial_time=time.mktime(simulation_start.timetuple()))
env.epoch = time.mktime(simulation_start.timetuple())
# Add graph to environment
env.FG = FG
# Add environment and path to the vessel
vessel.env = env
vessel.route = path
vessel.geometry = nx.get_node_attributes(FG, "geometry")[path[0]]
# Start the simulation
env.process(start(env, vessel))
env.run()
print("Simulation of path {} took {} seconds".format(path, int(env.now)))
Simulation of path [80, 81, 94, 95, 96, 303, 74, 304, 315, 395, 383, 253, 254, 255, 88, 89, 90, 109, 287] took 1717070818 seconds
Obtain vessel log information#
The cel below uses the vessel log. The core function log_entry is used, which takes four arguments:
Log. A text to describe what is logged.
t. The timestamp.
Value. The value for the log (for sailing this is the distance).
Geometry The location of the vessel while loggin.
vessel_log = pd.DataFrame.from_dict(vessel.logbook)
vessel_log.head()
| Message | Timestamp | Value | Geometry | |
|---|---|---|---|---|
| 0 | Start sailing | 2024-05-30 09:35:02.000000 | POINT (4.2898595 51.3462538) | |
| 1 | Sailing from node 80 to node 81 start | 2024-05-30 09:35:02.000000 | 0 | POINT (4.2898595 51.3462538) |
| 2 | Sailing from node 80 to node 81 stop | 2024-05-30 09:35:30.779383 | 0 | POINT (4.2898218 51.3459962) |
| 3 | Sailing from node 81 to node 94 start | 2024-05-30 09:35:30.779383 | 0 | POINT (4.2898218 51.3459962) |
| 4 | Sailing from node 81 to node 94 stop | 2024-05-30 09:44:03.014840 | 0 | POINT (4.2825009 51.3464185) |
Visualization of path#
If you get an error regarding ffmpeg use this answer. You have to install ffmpeg in your Conda environment. It can be done using the following command.
#Install ffmpeg using Conda
conda install -c conda-forge ffmpeg
%%time
# Get the coordinates of every step
coordinates = []
for i in vessel_log["Geometry"]:
coordinates.append((i.x, i.y))
# Get the time of every step
timesteps = []
for i in vessel_log["Timestamp"]:
timesteps.append(i.timestamp())
# Make the animation
fig = plt.figure(figsize=[12, 12])
nx.draw(FG, nx.get_node_attributes(FG, "pos"), with_labels=True, node_size=0.5, font_size=2, width=0.2, arrowsize=3)
(location,) = plt.plot([], [], "kx", markersize=15)
(path,) = plt.plot([], [], "go", markersize=8)
time_text = plt.title("")
def init():
location.set_data([], [])
path.set_data([], [])
time_text.set_text("Time is 0")
def animate(i):
this_x = [coordinates[i][0]]
this_y = [coordinates[i][1]]
location.set_data(this_x, this_y)
time_text.set_text("Time is {}".format(int(timesteps[i])))
if 0 < i:
past_x = [coordinate[0] for coordinate in coordinates[:i]]
past_y = [coordinate[1] for coordinate in coordinates[:i]]
path.set_data(past_x, past_y)
return location, path, time_text
else:
return location, time_text
ani = animation.FuncAnimation(fig, animate, np.arange(0, len(timesteps)), init_func=init)
ani.save("Example 06 - route.mp4", fps=5)
MovieWriter ffmpeg unavailable; using Pillow instead.
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/animation.py:243, in AbstractMovieWriter.saving(self, fig, outfile, dpi, *args, **kwargs)
242 try:
--> 243 yield self
244 finally:
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/animation.py:1105, in Animation.save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs, progress_callback)
1103 for anim, d in zip(all_anim, data):
1104 # TODO: See if turning off blit is really necessary
-> 1105 anim._draw_next_frame(d, blit=False)
1106 if progress_callback is not None:
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/animation.py:1141, in Animation._draw_next_frame(self, framedata, blit)
1140 self._draw_frame(framedata)
-> 1141 self._post_draw(framedata, blit)
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/animation.py:1166, in Animation._post_draw(self, framedata, blit)
1165 else:
-> 1166 self._fig.canvas.draw_idle()
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/backend_bases.py:1919, in FigureCanvasBase.draw_idle(self, *args, **kwargs)
1918 with self._idle_draw_cntx():
-> 1919 self.draw(*args, **kwargs)
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/backends/backend_agg.py:387, in FigureCanvasAgg.draw(self)
385 with (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
386 else nullcontext()):
--> 387 self.figure.draw(self.renderer)
388 # A GUI class may be need to update a window using this draw, so
389 # don't forget to call the superclass.
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/artist.py:95, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
93 @wraps(draw)
94 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 95 result = draw(artist, renderer, *args, **kwargs)
96 if renderer._rasterizing:
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
70 renderer.start_filter()
---> 72 return draw(artist, renderer)
73 finally:
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/figure.py:3155, in Figure.draw(self, renderer)
3154 self.patch.draw(renderer)
-> 3155 mimage._draw_list_compositing_images(
3156 renderer, self, artists, self.suppressComposite)
3158 renderer.close_group('figure')
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/image.py:132, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
131 for a in artists:
--> 132 a.draw(renderer)
133 else:
134 # Composite any adjacent images together
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
70 renderer.start_filter()
---> 72 return draw(artist, renderer)
73 finally:
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/axes/_base.py:3109, in _AxesBase.draw(self, renderer)
3107 _draw_rasterized(self.figure, artists_rasterized, renderer)
-> 3109 mimage._draw_list_compositing_images(
3110 renderer, self, artists, self.figure.suppressComposite)
3112 renderer.close_group('axes')
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/image.py:132, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
131 for a in artists:
--> 132 a.draw(renderer)
133 else:
134 # Composite any adjacent images together
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/artist.py:39, in _prevent_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
37 renderer._rasterizing = False
---> 39 return draw(artist, renderer, *args, **kwargs)
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/patches.py:4446, in FancyArrowPatch.draw(self, renderer)
4445 self._dpi_cor = renderer.points_to_pixels(1.)
-> 4446 path, fillable = self._get_path_in_displaycoord()
4448 if not np.iterable(fillable):
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/patches.py:4421, in FancyArrowPatch._get_path_in_displaycoord(self)
4420 (posA, posB) = self.get_transform().transform((posA, posB))
-> 4421 _path = self.get_connectionstyle()(posA, posB,
4422 patchA=self.patchA,
4423 patchB=self.patchB,
4424 shrinkA=self.shrinkA * dpi_cor,
4425 shrinkB=self.shrinkB * dpi_cor
4426 )
4427 else:
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/networkx/drawing/nx_pylab.py:813, in draw_networkx_edges.<locals>._draw_networkx_edges_fancy_arrow_patch.<locals>._connectionstyle(posA, posB, *args, **kwargs)
811 # if not, fall back to the user specified behavior
812 else:
--> 813 ret = base_connection_style(posA, posB, *args, **kwargs)
815 return ret
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/patches.py:2822, in ConnectionStyle._Base.__call__(self, posA, posB, shrinkA, shrinkB, patchA, patchB)
2817 path = self._clip(
2818 path,
2819 self._in_patch(patchA) if patchA else None,
2820 self._in_patch(patchB) if patchB else None,
2821 )
-> 2822 path = self._clip(
2823 path,
2824 inside_circle(*path.vertices[0], shrinkA) if shrinkA else None,
2825 inside_circle(*path.vertices[-1], shrinkB) if shrinkB else None
2826 )
2827 return path
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/patches.py:2800, in ConnectionStyle._Base._clip(self, path, in_start, in_stop)
2799 try:
-> 2800 _, path = split_path_inout(path, in_start)
2801 except ValueError:
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/bezier.py:378, in split_path_inout(path, inside, tolerance, reorder_inout)
377 bp = bezier_path.reshape((-1, 2))
--> 378 left, right = split_bezier_intersecting_with_closedpath(
379 bp, inside, tolerance)
380 if len(left) == 2:
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/bezier.py:341, in split_bezier_intersecting_with_closedpath(bezier, inside_closedpath, tolerance)
339 bezier_point_at_t = bz.point_at_t
--> 341 t0, t1 = find_bezier_t_intersecting_with_closedpath(
342 bezier_point_at_t, inside_closedpath, tolerance=tolerance)
344 _left, _right = split_de_casteljau(bezier, (t0 + t1) / 2.)
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/bezier.py:169, in find_bezier_t_intersecting_with_closedpath(bezier_point_at_t, inside_closedpath, t0, t1, tolerance)
168 middle_t = 0.5 * (t0 + t1)
--> 169 middle = bezier_point_at_t(middle_t)
170 middle_inside = inside_closedpath(middle)
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/bezier.py:230, in BezierSegment.point_at_t(self, t)
227 """
228 Evaluate the curve at a single point, returning a tuple of *d* floats.
229 """
--> 230 return tuple(self(t))
KeyboardInterrupt:
During handling of the above exception, another exception occurred:
KeyError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/PIL/Image.py:2436, in Image.save(self, fp, format, **params)
2435 try:
-> 2436 format = EXTENSION[ext]
2437 except KeyError as e:
KeyError: '.mp4'
The above exception was the direct cause of the following exception:
ValueError Traceback (most recent call last)
File <timed exec>:47
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/animation.py:1109, in Animation.save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs, progress_callback)
1107 progress_callback(frame_number, total_frames)
1108 frame_number += 1
-> 1109 writer.grab_frame(**savefig_kwargs)
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/contextlib.py:137, in _GeneratorContextManager.__exit__(self, typ, value, traceback)
135 value = typ()
136 try:
--> 137 self.gen.throw(typ, value, traceback)
138 except StopIteration as exc:
139 # Suppress StopIteration *unless* it's the same exception that
140 # was passed to throw(). This prevents a StopIteration
141 # raised inside the "with" statement from being suppressed.
142 return exc is not value
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/animation.py:245, in AbstractMovieWriter.saving(self, fig, outfile, dpi, *args, **kwargs)
243 yield self
244 finally:
--> 245 self.finish()
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib/animation.py:515, in PillowWriter.finish(self)
514 def finish(self):
--> 515 self._frames[0].save(
516 self.outfile, save_all=True, append_images=self._frames[1:],
517 duration=int(1000 / self.fps), loop=0)
File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/PIL/Image.py:2439, in Image.save(self, fp, format, **params)
2437 except KeyError as e:
2438 msg = f"unknown file extension: {ext}"
-> 2439 raise ValueError(msg) from e
2441 if format.upper() not in SAVE:
2442 init()
ValueError: unknown file extension: .mp4
Visualisation improved#
geom_x = []
geom_y = []
for geom in vessel_log["Geometry"]:
geom_x.append(geom.x)
geom_y.append(geom.y)
vessel_log["Geometry - x"] = geom_x
vessel_log["Geometry - y"] = geom_y
vessel_log.head()
| Message | Timestamp | Value | Geometry | Geometry - x | Geometry - y | |
|---|---|---|---|---|---|---|
| 0 | Start sailing | 2024-05-30 10:31:38.000000 | POINT (4.2898595 51.3462538) | 4.289860 | 51.346254 | |
| 1 | Sailing from node 80 to node 81 start | 2024-05-30 10:31:38.000000 | 0 | POINT (4.2898595 51.3462538) | 4.289860 | 51.346254 |
| 2 | Sailing from node 80 to node 81 stop | 2024-05-30 10:32:06.779383 | 0 | POINT (4.2898218 51.3459962) | 4.289822 | 51.345996 |
| 3 | Sailing from node 81 to node 94 start | 2024-05-30 10:32:06.779383 | 0 | POINT (4.2898218 51.3459962) | 4.289822 | 51.345996 |
| 4 | Sailing from node 81 to node 94 stop | 2024-05-30 10:40:39.014840 | 0 | POINT (4.2825009 51.3464185) | 4.282501 | 51.346418 |
geom_x = []
geom_y = []
for geom in vessel_log["Geometry"]:
geom_x.append(geom.x)
geom_y.append(geom.y)
vessel_log["Geometry - x"] = geom_x
vessel_log["Geometry - y"] = geom_y
time_stamp_min = min(vessel_log["Timestamp"]).timestamp()
time_stamp_max = max(vessel_log["Timestamp"]).timestamp()
steps = int(np.floor((time_stamp_max - time_stamp_min) / 60))
steps = vessel_log.shape[0]
timestamps_t = np.linspace(time_stamp_min, time_stamp_max, steps)
times = []
for t in vessel_log["Timestamp"]:
times.append(t.timestamp())
vessel_log["timestamps_t"] = timestamps_t
vessel_log["timestamps_x"] = np.interp(timestamps_t, times, vessel_log["Geometry - x"])
vessel_log["timestamps_y"] = np.interp(timestamps_t, times, vessel_log["Geometry - y"])
timestamps_t = vessel_log["timestamps_t"]
timestamps_x = vessel_log["timestamps_x"]
timestamps_y = vessel_log["timestamps_y"]
%%time
# Make the animation
fig = plt.figure(figsize=[12, 12])
nx.draw(FG, nx.get_node_attributes(FG, "pos"), with_labels=True, node_size=0.5, font_size=2, width=0.2, arrowsize=3)
(location,) = plt.plot([], [], "ko", markersize=15)
(path,) = plt.plot([], [], "g")
time_text = plt.title("")
def init():
location.set_data([], [])
path.set_data([], [])
time_text.set_text("Time is 0")
def animate(i):
this_x = [timestamps_x[i]]
this_y = [timestamps_x[i]]
location.set_data(this_x, this_y)
time_text.set_text("Time is {}".format(int(timestamps_t[i])))
if 0 < i:
past_x = [x for x in timestamps_x[:i]]
past_y = [y for y in timestamps_y[:i]]
path.set_data(past_x, past_y)
return location, path, time_text
else:
return location, time_text
ani = animation.FuncAnimation(fig, animate, np.arange(0, len(timestamps_t)), init_func=init)
ani.save("Example 06 - route - improved.mp4", fps=5, writer='ffmpeg')
CPU times: user 22.5 s, sys: 2 s, total: 24.5 s
Wall time: 21.8 s
