diff --git a/src/worldcereal/seasons.py b/src/worldcereal/seasons.py
index 1a343845..a0df9382 100644
--- a/src/worldcereal/seasons.py
+++ b/src/worldcereal/seasons.py
@@ -1,5 +1,6 @@
 import datetime
 import importlib.resources as pkg_resources
+import math
 
 import numpy as np
 import pandas as pd
@@ -16,14 +17,50 @@ class NoSeasonError(Exception):
     pass
 
 
-def doy_from_tiff(season: str, kind: str, bounds: tuple, epsg: str, resolution: int):
+def doy_to_angle(day_of_year, total_days=365):
+    return 2 * math.pi * (day_of_year / total_days)
+
+
+def angle_to_doy(angle, total_days=365):
+    return (angle / (2 * math.pi)) * total_days
+
+
+def circular_median_day_of_year(doy_array, total_days=365):
+    """This function computes the median doy from a given array
+    taking into account its circular nature. Still has to be used with caution!
+    """
+    angles = np.array([doy_to_angle(day, total_days) for day in doy_array])
+
+    def circular_distance(a, b):
+        return np.angle(np.exp(1j * (a - b)))
+
+    median_angle = None
+    min_distance = float("inf")
+
+    for angle in angles:
+        total_distance = np.sum(
+            [abs(circular_distance(angle, other)) for other in angles]
+        )
+        if total_distance < min_distance:
+            min_distance = total_distance
+            median_angle = angle
+
+    median_day = round(angle_to_doy(median_angle, total_days))
+    median_day = median_day % total_days
+    if median_day < 0:
+        median_day += total_days
+
+    return median_day
+
+
+def doy_from_tiff(season: str, kind: str, bonds: tuple, epsg: int, resolution: int):
     """Function to read SOS/EOS DOY from TIFF
 
     Args:
         season (str): crop season to process
         kind (str): which DOY to read (SOS/EOS)
         bounds (tuple): the bounds to read
-        epsg (str): epsg in which bounds are expressed
+        epsg (int): epsg in which bounds are expressed
         resolution (int): resolution in meters of resulting array
 
     Raises:
@@ -35,6 +72,11 @@ def doy_from_tiff(season: str, kind: str, bounds: tuple, epsg: str, resolution:
         np.ndarray: resulting DOY array
     """
 
+    if epsg == 4326:
+        raise ValueError(
+            "EPSG 4326 not supported for DOY data. Use a projected CRS instead."
+        )
+
     if season not in SUPPORTED_SEASONS:
         raise ValueError(f"Season `{season}` not supported.")
     else:
@@ -53,7 +95,9 @@ def doy_from_tiff(season: str, kind: str, bounds: tuple, epsg: str, resolution:
     logger.info(f"Loading DOY data from: {doy_file}")
 
     with pkg_resources.open_binary(cropcalendars, doy_file) as doy_file:
-        doy_data = load_reproject(doy_file, bounds, epsg, resolution)
+        doy_data = load_reproject(
+            doy_file, bounds, epsg, resolution, nodata_value=0, fill_value=np.nan
+        )
 
     return doy_data
 
@@ -189,14 +233,18 @@ def get_processing_dates_for_extent(
     if season not in SUPPORTED_SEASONS:
         raise ValueError(f"Season `{season}` not supported!")
 
-    bounds = (extent.east, extent.south, extent.west, extent.north)
+    bounds = (extent.west, extent.south, extent.east, extent.north)
     epsg = extent.epsg
 
     # Get DOY of EOS for the target season
-    eos_doy = doy_from_tiff(season, "EOS", bounds, epsg, 10000)
+    eos_doy = doy_from_tiff(season, "EOS", bounds, epsg, 10000).flatten()
+
+    # Check if we have seasonality
+    if not np.isfinite(eos_doy).any():
+        raise NoSeasonError(f"No valid EOS DOY found for season `{season}`")
 
-    # Get the median EOS
-    eos_doy_median = np.median(eos_doy)
+    # Compute median DOY
+    eos_doy_median = circular_median_day_of_year(eos_doy)
 
     # We can derive the end date from year and EOS
     end_date = datetime.datetime(year, 1, 1) + pd.Timedelta(days=eos_doy_median)
diff --git a/src/worldcereal/utils/geoloader.py b/src/worldcereal/utils/geoloader.py
index 53365269..9af73e39 100644
--- a/src/worldcereal/utils/geoloader.py
+++ b/src/worldcereal/utils/geoloader.py
@@ -115,5 +115,3 @@ def load_reproject(
         dst = dst[border_buff:-border_buff, border_buff:-border_buff]
 
     return dst
-
-    return dst
diff --git a/tests/worldcerealtests/test_seasons.py b/tests/worldcerealtests/test_seasons.py
index 31e744f6..71be2a0e 100644
--- a/tests/worldcerealtests/test_seasons.py
+++ b/tests/worldcerealtests/test_seasons.py
@@ -1,9 +1,11 @@
 import datetime
 
+import numpy as np
 import pandas as pd
 
 from worldcereal import BoundingBoxExtent
 from worldcereal.seasons import (
+    circular_median_day_of_year,
     doy_from_tiff,
     doy_to_date_after,
     get_processing_dates_for_extent,
@@ -43,7 +45,7 @@ def test_doy_to_date_after():
 def test_get_processing_dates_for_extent():
     # Test to check if we can infer processing dates for default season
     # tc-annual
-    bounds = (574680, 5621800, 575320, 5622440)
+    bounds = (167286, 553423, 943774, 997257)
     epsg = 32631
     year = 2021
     extent = BoundingBoxExtent(*bounds, epsg)
@@ -54,3 +56,15 @@ def test_get_processing_dates_for_extent():
     assert pd.to_datetime(end_date) - pd.to_datetime(start_date) == pd.Timedelta(
         days=364
     )
+
+
+def test_compute_median_doy():
+    # Test to check if we can compute median day of year
+    # if array crosses the calendar year
+    doy_array = np.array([360, 362, 365, 1, 3, 5, 7])
+    assert circular_median_day_of_year(doy_array) == 1
+
+    # Other tests
+    assert circular_median_day_of_year([1, 2, 3, 4, 5]) == 3
+    assert circular_median_day_of_year([320, 330, 340, 360]) == 330
+    assert circular_median_day_of_year([320, 330, 340, 360, 10]) == 340