diff --git a/README.md b/README.md
index d955102..97437a5 100644
--- a/README.md
+++ b/README.md
@@ -148,7 +148,7 @@ plot_reservoir(
 Plot a particular output variable (as defined in `config.yaml`) over time:
 ```python
 from mosartwmpy import Model
-from mosartwmpy.plotting.plot import plot_reservoir
+from mosartwmpy.plotting.plot import plot_variable
 config_file = 'config.yaml'
 mosart_wm = Model()
 mosart_wm.initialize(config_file)
diff --git a/mosartwmpy/plotting/plot.py b/mosartwmpy/plotting/plot.py
index e2e5e31..9386ef3 100644
--- a/mosartwmpy/plotting/plot.py
+++ b/mosartwmpy/plotting/plot.py
@@ -1,5 +1,4 @@
 import contextily as ctx
-from epiweeks import Week
 import geopandas as gpd
 import matplotlib.pyplot as plt
 import numpy as np
@@ -12,6 +11,7 @@
 import warnings
 
 from mosartwmpy import Model
+from mosartwmpy.utilities.epiweek import get_epiweek_from_datetime
 
 warnings.filterwarnings("ignore", category=ShapelyDeprecationWarning)
 
@@ -186,7 +186,7 @@ def plot_reservoir(
 
     params = xr.open_dataset(model.config.get('water_management.reservoirs.parameters.path')).to_dataframe()
     dam = params[params[model.config.get('water_management.reservoirs.parameters.variables.reservoir_id')] == grand_id]
-    weeks = df.index.to_series().apply(lambda time: Week.fromdate(pd.to_datetime(time)).week)
+    weeks = df.index.to_series().apply(lambda time: get_epiweek_from_datetime(pd.to_datetime(time)))
     weeks[weeks > 52] = 52
 
     meanflow = pd.read_parquet(model.config.get('water_management.reservoirs.streamflow.path'))
diff --git a/mosartwmpy/reservoirs/istarf.py b/mosartwmpy/reservoirs/istarf.py
index 96acb60..60788d8 100644
--- a/mosartwmpy/reservoirs/istarf.py
+++ b/mosartwmpy/reservoirs/istarf.py
@@ -2,19 +2,19 @@
 import numpy as np
 
 from datetime import datetime
-from epiweeks import Week
 from benedict.dicts import benedict as Benedict
 
 from mosartwmpy.config.parameters import Parameters
 from mosartwmpy.state.state import State
 from mosartwmpy.grid.grid import Grid
+from mosartwmpy.utilities.epiweek import get_epiweek_from_datetime
 
 
 def istarf_release(state: State, grid: Grid, config: Benedict, parameters: Parameters, current_time: datetime):
     # estimate reservoir release using ISTARF which is based on harmonic functions
 
     # restrict epiweek to [1, 52]
-    epiweek = np.minimum(float(Week.fromdate(current_time).week), 52.0)
+    epiweek = np.minimum(float(get_epiweek_from_datetime(current_time)), 52.0)
 
     daily_release = np.zeros(len(grid.reservoir_id))
 
diff --git a/mosartwmpy/tests/test_epiweek.py b/mosartwmpy/tests/test_epiweek.py
new file mode 100644
index 0000000..5055020
--- /dev/null
+++ b/mosartwmpy/tests/test_epiweek.py
@@ -0,0 +1,18 @@
+from datetime import datetime
+import unittest
+
+from mosartwmpy.utilities.epiweek import get_epiweek_from_datetime
+
+
+class EpiweekTest(unittest.TestCase):
+
+    def test_epiweek(self):
+        self.assertEqual(get_epiweek_from_datetime(datetime(2017, 1, 1)), 1)
+        self.assertEqual(get_epiweek_from_datetime(datetime(2015, 1, 1)), 53)
+        self.assertEqual(get_epiweek_from_datetime(datetime(2016, 1, 1)), 52)
+        self.assertEqual(get_epiweek_from_datetime(datetime(2017, 3, 15)), 11)
+        self.assertEqual(get_epiweek_from_datetime(datetime(2017, 3, 31)), 13)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/mosartwmpy/tests/test_istarf.py b/mosartwmpy/tests/test_istarf.py
index 6892938..cf8e5c3 100644
--- a/mosartwmpy/tests/test_istarf.py
+++ b/mosartwmpy/tests/test_istarf.py
@@ -2,16 +2,16 @@
 import unittest
 
 from datetime import date
-from epiweeks import Week
 
 from mosartwmpy.reservoirs.istarf import compute_istarf_release
+from mosartwmpy.utilities.epiweek import get_epiweek_from_datetime
 
 
 class IstarfTest(unittest.TestCase):
 
     def test_is_correct(self):
 
-        epiweek = np.minimum(float(Week.fromdate(date(1994, 10, 1)).week), 52.0)
+        epiweek = np.minimum(float(get_epiweek_from_datetime(date(1994, 10, 1))), 52.0)
         reservoir_id = np.array([0.0])
         upper_min = np.array([-np.Inf])
         upper_max = np.array([np.Inf])
diff --git a/mosartwmpy/utilities/epiweek.py b/mosartwmpy/utilities/epiweek.py
new file mode 100644
index 0000000..7acba7b
--- /dev/null
+++ b/mosartwmpy/utilities/epiweek.py
@@ -0,0 +1,37 @@
+from datetime import date, datetime
+
+
+def get_epiweek_from_datetime(dt: datetime) -> int:
+    """Calculate the epiweek number from a datetime
+
+    Args:
+        dt (datetime): the datetime object to calculate the epiweek for
+
+    Returns:
+        int: the epiweek number
+
+    """
+    year = dt.year
+    dt_ordinal = dt.toordinal()
+    year_start_ordinal = _get_year_start_ordinal(year)
+    week = (dt_ordinal - year_start_ordinal) // 7
+    if week < 0:
+        year = year - 1
+        year_start_ordinal = _get_year_start_ordinal(year)
+        week = (dt_ordinal - year_start_ordinal) // 7
+    elif week >= 52:
+        year_start_ordinal = _get_year_start_ordinal(year + 1)
+        if date_ordinal > year_start_ordinal:
+            year = year + 1
+            week = 0
+    week = week + 1
+    return week
+
+
+def _get_year_start_ordinal(y):
+    year_start = date(y, 1, 1)
+    year_start_weekday = year_start.weekday()
+    year_start_ordinal = year_start.toordinal() - year_start_weekday - 1
+    if year_start_weekday > 2:
+        return year_start_ordinal + 7
+    return year_start_ordinal
\ No newline at end of file
diff --git a/setup.py b/setup.py
index e26e024..e784665 100644
--- a/setup.py
+++ b/setup.py
@@ -34,7 +34,6 @@ def readme():
         'click==8.0.1',
         'contextily==1.2.0',
         'dask[complete]==2021.10.0',
-        'epiweeks==2.1.3',
         'geopandas==0.10.2',
         'h5netcdf==0.11.0',
         'hvplot==0.7.3',