diff --git a/bar_plot.py b/bar_plot.py
index a93c897dbecc5fa5e1829384a2341b5c512af4e2..15807dce746a362e8bdabf654104081de8f4c488 100644
--- a/bar_plot.py
+++ b/bar_plot.py
@@ -5,41 +5,39 @@
 Reads data of oil spill arrival time from shipping routes to different sites and computes its monthly mean and standard deviation
 
 Author: Thomas Dobbelaere, Earth and Life Institute, UCLouvain, Belgium
-Last modified: 19 July 2022
+Last modified: 9 November 2022
 """
+import pandas as pd
+import calendar
 import matplotlib.pyplot as plt
 import numpy as np
-from matplotlib import pyplot
-import calendar
-import matplotlib as mpl
 
 months = [calendar.month_abbr[i] for i in range(1,13)]
-
-fd = np.load("arrival_by_route.npy", allow_pickle=True).item()
-routes = fd["routes"]
-meanv = np.nanmean(fd["data"], axis=2)
-stdv = np.nanstd(fd["data"], axis=2)
-route_names = {
-    "Barhain" : "Barhain",
-    "Iran" : "Qatar-Iran",
-    "Ras_laffan_field": 'Ras Laffan to North Field',
-    "Ras_laffan": "Ras Laffan LNG export",
-    "Doha" : 'Doha-Mesaieed'
-}
-x = np.arange(len(months))
+x = np.arange(12)
+routes = ["Bahrain", "Qatar-Iran", 'Ras Laffan to North Field',"Ras Laffan LNG export",'Doha-Mesaieed']
 width = 0.18
 offset = np.arange(-2,3)*width
-clrs = mpl.cm.Reds_r(0.1 + 0.8*np.linspace(0,1,5))
+cmap = plt.get_cmap("Reds_r")
+clrs = cmap(0.1 + 0.8*np.linspace(0,1,5))
+
+df = pd.read_csv("arrival_time_route.csv")
+df['month'] = df.month.apply(lambda x: months.index(x))
+data_by_route = df.groupby(by=['route'])
+means = df.groupby(by=['month','route']).mean()
 
 fig, ax = plt.subplots(figsize=(14,7))
-ibar = 0
-for key,label in route_names.items():
-    row = routes.index(key)
-    ax.bar( x+offset[ibar], meanv[row], width=width, yerr=stdv[row], \
-                color=clrs[ibar], ecolor="k", capsize=3, label=label)
-    ibar += 1
 
-ax.set_ylabel('Mean arrival time (days)', fontsize=14)
+for i,route in enumerate(routes):
+    meanv = [means.loc[(i,route),'arrival_time'] for i in range(12)]
+    ax.bar(
+        x+offset[i], meanv, width=width, color=clrs[i], ecolor='k', capsize=3, label=route
+    )
+    grp = data_by_route.get_group(route)
+    ax.scatter(x[grp.month]+offset[i], grp.arrival_time, s=10, facecolors='none', edgecolors="k")
+
+ax.set_ylabel('Arrival time (days)', fontsize=14)
+ax.set_yticks([0,1,2,3,4])
+ax.set_ylim(0,5)
 ax.set_yticklabels((0,1,2,3,4),fontsize=14)
 ax.set_xticks(x)
 ax.set_xticklabels(months, fontsize=14)
@@ -47,6 +45,5 @@ ax.spines['top'].set_visible(False)
 ax.spines['right'].set_visible(False)
 for i in range(1,5):
     ax.axhline(y=i, color='grey', linestyle=(0, (2, 10)))
-ax.legend(loc=1, fontsize=12)
-
-plt.savefig('mean_arrival_by_route.png', bbox_inches='tight', dpi=200)
+ax.legend(loc='upper center', fontsize=12,  bbox_to_anchor=(0.5, -0.075), ncol=5)
+plt.savefig("barplot_mean_dots.pdf", dpi=300, bbox_inches='tight')
\ No newline at end of file