diff --git a/cats/CI_api_interface.py b/cats/CI_api_interface.py
index dc3e3fb..8942a17 100644
--- a/cats/CI_api_interface.py
+++ b/cats/CI_api_interface.py
@@ -1,6 +1,7 @@
 import sys
 from collections import namedtuple
 from datetime import datetime
+from zoneinfo import ZoneInfo
 
 from .forecast import CarbonIntensityPointEstimate
 
@@ -63,9 +64,13 @@ def invalid_code(r: dict):
         raise InvalidLocationError
 
     datefmt = "%Y-%m-%dT%H:%MZ"
+    # The "Z" at the end of the format string indicates UTC,
+    # however, strptime does not know how to parse this, so we 
+    # need to add tzinfo data.
+    utc = ZoneInfo('UTC')
     return [
         CarbonIntensityPointEstimate(
-            datetime=datetime.strptime(d["from"], datefmt),
+            datetime=datetime.strptime(d["from"], datefmt).replace(tzinfo=utc),
             value=d["intensity"]["forecast"],
         )
         for d in response["data"]["data"]
diff --git a/cats/optimise_starttime.py b/cats/optimise_starttime.py
index 0c8a8c4..b497479 100644
--- a/cats/optimise_starttime.py
+++ b/cats/optimise_starttime.py
@@ -16,5 +16,10 @@ def get_avg_estimates(data, duration=None):
     #     raise ValueError(
     #         "Windowed method timespan cannot be greater than the cached timespan"
     #     )
-    wf = WindowedForecast(data, duration, start=datetime.now())
+    # NB: datetime.now().astimezone() gives us a timezone aware datetime object
+    # in the current system timezone. The resulting start time from the forecast
+    # works in terms of this timezone (even if the data is in another timezone)
+    # so we end up with something in system time if data is in utc and system
+    # time is in bst (for example)
+    wf = WindowedForecast(data, duration, start=datetime.now().astimezone())
     return wf[0], min(wf)
diff --git a/tests/test_api_query.py b/tests/test_api_query.py
index 8fb745d..719647b 100644
--- a/tests/test_api_query.py
+++ b/tests/test_api_query.py
@@ -6,7 +6,10 @@
 
 def test_api_call():
     """
-    This just checks the API call runs and returns a list of tuples
+    This just checks the API call runs and returns a list of point estimates
+
+    Also confirms that datetime objects are timezone aware, as per
+    https://docs.python.org/3/library/datetime.html#determining-if-an-object-is-aware-or-naive
     """
 
     api_interface = API_interfaces["carbonintensity.org.uk"]
@@ -15,6 +18,8 @@ def test_api_call():
     assert isinstance(response, list)
     for item in response:
         assert isinstance(item, CarbonIntensityPointEstimate)
+        assert ((item.datetime.tzinfo is not None) and
+                (item.datetime.tzinfo.utcoffset(item.datetime) is not None))
 
 def test_bad_postcode():
     api_interface = API_interfaces["carbonintensity.org.uk"]
@@ -24,3 +29,4 @@ def test_bad_postcode():
 
     with pytest.raises(InvalidLocationError):
         response = cats.CI_api_query.get_CI_forecast('A', api_interface)
+
diff --git a/tests/test_windowed_forecast.py b/tests/test_windowed_forecast.py
index 8d59aad..faaafb6 100644
--- a/tests/test_windowed_forecast.py
+++ b/tests/test_windowed_forecast.py
@@ -1,5 +1,6 @@
 import csv
-from datetime import datetime, timedelta
+from datetime import datetime, timedelta, timezone
+from zoneinfo import ZoneInfo
 import math
 from pathlib import Path
 from numpy.testing import assert_allclose
@@ -61,7 +62,7 @@ def test_minimise_average():
         next(csvfile)  # Skip header line
         data = [
             CarbonIntensityPointEstimate(
-                datetime=datetime.fromisoformat(datestr[:-1]),
+                datetime=datetime.fromisoformat(datestr[:-1]+'+00:00'),
                 value=float(intensity_value),
             )
             for datestr, _, _, intensity_value in csvfile
@@ -75,8 +76,8 @@ def test_minimise_average():
         # Intensity point estimates over best runtime period
         v = [10, 8, 7, 7, 5, 8, 8]
         expected = CarbonIntensityAverageEstimate(
-            start=datetime.fromisoformat("2023-05-05T12:00"),
-            end=datetime.fromisoformat("2023-05-05T15:00"),
+            start=datetime.fromisoformat("2023-05-05T12:00+00:00"),
+            end=datetime.fromisoformat("2023-05-05T15:00+00:00"),
             value=sum(
                 [0.5 * (a + b) for a, b in zip(v[:-1], v[1:])]
             ) / window_size
@@ -84,13 +85,49 @@ def test_minimise_average():
         assert (result == expected)
 
 
+def test_minimise_average_bst():
+    # We should get a start time in BST if we provide the starting time
+    # in that timezone, even if the intensity estimate is in UTC. This
+    # is needed as the `at` command works in local system time (and that's
+    # what we put in)
+    with open(TEST_DATA, "r") as f:
+        csvfile = csv.reader(f, delimiter=",")
+        next(csvfile)  # Skip header line
+        data = [
+            CarbonIntensityPointEstimate(
+                datetime=datetime.fromisoformat(datestr[:-1]+'+00:00'),
+                value=float(intensity_value),
+            )
+            for datestr, _, _, intensity_value in csvfile
+        ]
+
+        window_size = 6
+        # Data points separated by 30 minutes intervals
+        duration = window_size * 30
+        start_time_bst = data[0].datetime.replace(tzinfo=timezone(timedelta(seconds=-3600)))
+        result = min(WindowedForecast(data, duration, start=start_time_bst))
+
+        # Intensity point estimates over best runtime period
+        v = [10, 8, 7, 7, 5, 8, 8]
+        expected = CarbonIntensityAverageEstimate(
+            start=datetime.fromisoformat("2023-05-05T11:00-01:00"),
+            end=datetime.fromisoformat("2023-05-05T14:00-01:00"),
+            value=sum(
+                [0.5 * (a + b) for a, b in zip(v[:-1], v[1:])]
+            ) / window_size
+        )
+        assert (result == expected)
+        assert (result.start.tzinfo == expected.start.tzinfo)
+        assert (result.end.tzinfo == expected.end.tzinfo)
+
+
 def test_average_intensity_now():
     with open(TEST_DATA, "r") as f:
         csvfile = csv.reader(f, delimiter=",")
         next(csvfile)  # Skip header line
         data = [
             CarbonIntensityPointEstimate(
-                datetime=datetime.fromisoformat(datestr[:-1]),
+                datetime=datetime.fromisoformat(datestr[:-1]+'+00:00'),
                 value=float(intensity_value),
             )
             for datestr, _, _, intensity_value in csvfile
@@ -118,24 +155,25 @@ def test_average_intensity_with_offset():
     # carbon intensity data points. In this case cats interpolate the
     # intensity value at beginning and end of each potential job
     # duration window.
+    utc = ZoneInfo('UTC')
     CI_forecast = [
-        CarbonIntensityPointEstimate(26, datetime(2023,1,1,8,30)),
-        CarbonIntensityPointEstimate(40, datetime(2023,1,1,9,0)),
-        CarbonIntensityPointEstimate(50, datetime(2023,1,1,9,30)),
-        CarbonIntensityPointEstimate(60, datetime(2023,1,1,10,0)),
-        CarbonIntensityPointEstimate(25, datetime(2023,1,1,10,30)),
+        CarbonIntensityPointEstimate(26, datetime(2023,1,1,8,30,tzinfo=utc)),
+        CarbonIntensityPointEstimate(40, datetime(2023,1,1,9,0,tzinfo=utc)),
+        CarbonIntensityPointEstimate(50, datetime(2023,1,1,9,30,tzinfo=utc)),
+        CarbonIntensityPointEstimate(60, datetime(2023,1,1,10,0,tzinfo=utc)),
+        CarbonIntensityPointEstimate(25, datetime(2023,1,1,10,30,tzinfo=utc)),
     ]
     duration = 70  # in minutes
     # First available data point is for 08:00 but the job
     # starts 15 minutes later.
-    job_start = datetime.fromisoformat("2023-01-01T08:45")
+    job_start = datetime.fromisoformat("2023-01-01T08:45+00:00")
     result = WindowedForecast(CI_forecast, duration, start=job_start)[1]
 
     interp1 = 40 + 15 * (50 - 40) / 30
     interp2 = 60 + 25 * (25 - 60) / 30
     expected = CarbonIntensityAverageEstimate(
-        start=datetime(2023,1,1,9,15),
-        end=datetime(2023,1,1,10,25),
+        start=datetime(2023,1,1,9,15,tzinfo=utc),
+        end=datetime(2023,1,1,10,25,tzinfo=utc),
         value=(
             0.5 * (interp1 + 50) * 15 +
             0.5 * (50 + 60) * 30 +
@@ -151,7 +189,7 @@ def test_average_intensity_with_offset():
 
     # When start at 09:15 the WindowedForecast's 'data' attribute
     # should discard the first data point at 08:30.
-    job_start = datetime.fromisoformat("2023-01-01T09:15")
+    job_start = datetime.fromisoformat("2023-01-01T09:15+00:00")
     result = WindowedForecast(CI_forecast, duration, start=job_start)[0]
     assert result == expected
 
@@ -165,7 +203,7 @@ def test_average_intensity_with_offset_long_job():
         next(csvfile)  # Skip header line
         data = [
             CarbonIntensityPointEstimate(
-                datetime=datetime.fromisoformat(datestr[:-1]),
+                datetime=datetime.fromisoformat(datestr[:-1]+'+00:00'),
                 value=float(intensity_value),
             )
             for datestr, _, _, intensity_value in csvfile
@@ -174,7 +212,8 @@ def test_average_intensity_with_offset_long_job():
         duration = 194  # in minutes
         # First available data point is for 12:30 but the job
         # starts 18 minutes later.
-        job_start = datetime.fromisoformat("2023-05-04T12:48")
+        # Start time in BST
+        job_start = datetime.fromisoformat("2023-05-04T13:48+01:00")
         result = WindowedForecast(data, duration, start=job_start)[2]
 
         # First and last element in v are interpolated intensity value.
@@ -191,6 +230,8 @@ def test_average_intensity_with_offset_long_job():
             ) / duration
         )
         assert (result == expected)
+        assert (result.start.tzinfo == expected.start.tzinfo)
+        assert (result.end.tzinfo == expected.end.tzinfo)
 
 def test_average_intensity_with_offset_short_job():
     # Case where job is short: start and end time fall between two
@@ -200,7 +241,7 @@ def test_average_intensity_with_offset_short_job():
         next(csvfile)  # Skip header line
         data = [
             CarbonIntensityPointEstimate(
-                datetime=datetime.fromisoformat(datestr[:-1]),
+                datetime=datetime.fromisoformat(datestr[:-1]+"+00:00"),
                 value=float(intensity_value),
             )
             for datestr, _, _, intensity_value in csvfile
@@ -209,7 +250,7 @@ def test_average_intensity_with_offset_short_job():
         duration = 6  # in minutes
         # First available data point is for 12:30 but the job
         # starts 6 minutes later.
-        job_start = datetime.fromisoformat("2023-05-04T12:48")
+        job_start = datetime.fromisoformat("2023-05-04T12:48+00:00")
         result = WindowedForecast(data, duration, start=job_start)[2]
 
         # Job starts at 12:48 and ends at 12:54. For each candidate