diff --git a/.github/workflows/lint.yaml b/.github/workflows/lint.yaml
index f483af9..d2c348c 100644
--- a/.github/workflows/lint.yaml
+++ b/.github/workflows/lint.yaml
@@ -24,7 +24,7 @@ jobs:
         python -m pip install isort black flake8
     - name: isort
       run: python -m isort . --check --diff
-    # - name: black
-    #   run: python -m black --check --diff .
+    - name: black
+      run: python -m black --check --diff .
     # - name: flake8
     #   run: python -m flake8 .
diff --git a/pyproject.toml b/pyproject.toml
index e99e514..2e85df4 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -65,3 +65,7 @@ addopts = [
 
 [tool.isort]
 profile = "black"
+
+
+[tool.black]
+line-length = 99
diff --git a/src/openeo_gfmap/backend.py b/src/openeo_gfmap/backend.py
index c38c5a3..52a103f 100644
--- a/src/openeo_gfmap/backend.py
+++ b/src/openeo_gfmap/backend.py
@@ -52,9 +52,7 @@ def vito_connection(capfd: Optional = None) -> openeo.Connection:
 
 def cdse_connection(capfd: Optional = None) -> openeo.Connection:
     """Performs a connection to the CDSE backend using oidc authentication."""
-    connection = openeo.connect(
-        "https://openeo.dataspace.copernicus.eu/openeo/1.2"
-    )
+    connection = openeo.connect("https://openeo.dataspace.copernicus.eu/openeo/1.2")
     if capfd is not None:
         with capfd.disabled():
             # Temporarily disable output capturing, to make sure that OIDC device
diff --git a/src/openeo_gfmap/extractions/s2.py b/src/openeo_gfmap/extractions/s2.py
index bce191c..342f963 100644
--- a/src/openeo_gfmap/extractions/s2.py
+++ b/src/openeo_gfmap/extractions/s2.py
@@ -76,9 +76,7 @@ def s2_l2a_fetch_default(
             ), "CRS not defined within GeoJSON collection."
             spatial_extent = dict(spatial_extent)
 
-        cube = connection.load_collection(
-            collection_name, spatial_extent, temporal_extent, bands
-        )
+        cube = connection.load_collection(collection_name, spatial_extent, temporal_extent, bands)
 
         # Apply if the collection is a GeoJSON Feature collection
         if isinstance(spatial_extent, GeoJSON):
diff --git a/src/openeo_gfmap/fetching/__init__.py b/src/openeo_gfmap/fetching/__init__.py
index f1bbd9f..fccda37 100644
--- a/src/openeo_gfmap/fetching/__init__.py
+++ b/src/openeo_gfmap/fetching/__init__.py
@@ -10,6 +10,8 @@
 from .s2 import build_sentinel2_l2a_extractor
 
 __all__ = [
-    "build_sentinel2_l2a_extractor", "CollectionFetcher", "FetchType",
-    "build_sentinel1_grd_extractor"
+    "build_sentinel2_l2a_extractor",
+    "CollectionFetcher",
+    "FetchType",
+    "build_sentinel1_grd_extractor",
 ]
diff --git a/src/openeo_gfmap/fetching/commons.py b/src/openeo_gfmap/fetching/commons.py
index a589dc8..898c86d 100644
--- a/src/openeo_gfmap/fetching/commons.py
+++ b/src/openeo_gfmap/fetching/commons.py
@@ -128,8 +128,7 @@ def load_collection(
     pre_mask = params.get("pre_mask", None)
     if pre_mask is not None:
         assert isinstance(pre_mask, openeo.DataCube), (
-            f"The 'pre_mask' parameter must be an openeo datacube, "
-            f"got {pre_mask}."
+            f"The 'pre_mask' parameter must be an openeo datacube, " f"got {pre_mask}."
         )
         cube = cube.mask(pre_mask.resample_cube_spatial(cube))
 
diff --git a/src/openeo_gfmap/fetching/s1.py b/src/openeo_gfmap/fetching/s1.py
index 4eac573..1dd2fee 100644
--- a/src/openeo_gfmap/fetching/s1.py
+++ b/src/openeo_gfmap/fetching/s1.py
@@ -84,9 +84,7 @@ def s1_grd_fetch_default(
     return s1_grd_fetch_default
 
 
-def get_s1_grd_default_processor(
-    collection_name: str, fetch_type: FetchType
-) -> Callable:
+def get_s1_grd_default_processor(collection_name: str, fetch_type: FetchType) -> Callable:
     """Builds the preprocessing function from the collection name as it is stored
     in the target backend.
     """
@@ -94,7 +92,7 @@ def get_s1_grd_default_processor(
     def s1_grd_default_processor(cube: openeo.DataCube, **params):
         """Default collection preprocessing method.
         This method performs:
-        
+
         * Compute the backscatter of all the S1 products. By default, the
         "sigma0-ellipsoid" method is used with "COPERNICUS_30" DEM, but those
         can be changed by specifying "coefficient" and "elevation_model" in
@@ -113,9 +111,7 @@ def s1_grd_default_processor(cube: openeo.DataCube, **params):
         )
 
         cube = resample_reproject(
-            cube,
-            params.get("target_resolution", 10.0),
-            params.get("target_crs", None) 
+            cube, params.get("target_resolution", 10.0), params.get("target_crs", None)
         )
 
         cube = rename_bands(cube, BASE_SENTINEL1_GRD_MAPPING)
@@ -127,28 +123,20 @@ def s1_grd_default_processor(cube: openeo.DataCube, **params):
 
 SENTINEL1_GRD_BACKEND_MAP = {
     Backend.TERRASCOPE: {
-        "default": partial(
-            get_s1_grd_default_fetcher, collection_name="SENTINEL1_GRD"
-        ),
-        "preprocessor": partial(
-            get_s1_grd_default_processor, collection_name="SENTINEL1_GRD"
-        )
+        "default": partial(get_s1_grd_default_fetcher, collection_name="SENTINEL1_GRD"),
+        "preprocessor": partial(get_s1_grd_default_processor, collection_name="SENTINEL1_GRD"),
     },
     Backend.CDSE: {
-        "default": partial(
-            get_s1_grd_default_fetcher, collection_name="SENTINEL1_GRD"
-        ),
-        "preprocessor": partial(
-            get_s1_grd_default_processor, collection_name="SENTINEL1_GRD"
-        )
-    }
+        "default": partial(get_s1_grd_default_fetcher, collection_name="SENTINEL1_GRD"),
+        "preprocessor": partial(get_s1_grd_default_processor, collection_name="SENTINEL1_GRD"),
+    },
 }
 
 
 def build_sentinel1_grd_extractor(
     backend_context: BackendContext, bands: list, fetch_type: FetchType, **params
 ) -> CollectionFetcher:
-    """ Creates a S1 GRD collection extractor for the given backend."""
+    """Creates a S1 GRD collection extractor for the given backend."""
     backend_functions = SENTINEL1_GRD_BACKEND_MAP.get(backend_context.backend)
 
     fetcher, preprocessor = (
@@ -156,6 +144,4 @@ def build_sentinel1_grd_extractor(
         backend_functions["preprocessor"](fetch_type=fetch_type),
     )
 
-    return CollectionFetcher(
-        backend_context, bands, fetcher, preprocessor, **params
-    )
+    return CollectionFetcher(backend_context, bands, fetcher, preprocessor, **params)
diff --git a/src/openeo_gfmap/fetching/s2.py b/src/openeo_gfmap/fetching/s2.py
index de7e201..0f7b3b3 100644
--- a/src/openeo_gfmap/fetching/s2.py
+++ b/src/openeo_gfmap/fetching/s2.py
@@ -114,9 +114,7 @@ def s2_l2a_fetch_default(
     return s2_l2a_fetch_default
 
 
-def get_s2_l2a_element84_fetcher(
-    collection_name: str, fetch_type: FetchType
-) -> Callable:
+def get_s2_l2a_element84_fetcher(collection_name: str, fetch_type: FetchType) -> Callable:
     """Fetches the collections from the Sentinel-2 Cloud-Optimized GeoTIFFs
     bucket provided by Amazon and managed by Element84.
     """
@@ -157,9 +155,7 @@ def s2_l2a_element84_fetcher(
     return s2_l2a_element84_fetcher
 
 
-def get_s2_l2a_default_processor(
-    collection_name: str, fetch_type: FetchType
-) -> Callable:
+def get_s2_l2a_default_processor(collection_name: str, fetch_type: FetchType) -> Callable:
     """Builds the preprocessing function from the collection name as it stored
     in the target backend.
     """
@@ -188,15 +184,11 @@ def s2_l2a_default_processor(cube: openeo.DataCube, **params):
 SENTINEL2_L2A_BACKEND_MAP = {
     Backend.TERRASCOPE: {
         "fetch": partial(get_s2_l2a_default_fetcher, collection_name="SENTINEL2_L2A"),
-        "preprocessor": partial(
-            get_s2_l2a_default_processor, collection_name="SENTINEL2_L2A"
-        ),
+        "preprocessor": partial(get_s2_l2a_default_processor, collection_name="SENTINEL2_L2A"),
     },
     Backend.CDSE: {
         "fetch": partial(get_s2_l2a_default_fetcher, collection_name="SENTINEL2_L2A"),
-        "preprocessor": partial(
-            get_s2_l2a_default_processor, collection_name="SENTINEL2_L2A"
-        ),
+        "preprocessor": partial(get_s2_l2a_default_processor, collection_name="SENTINEL2_L2A"),
     },
 }
 
diff --git a/src/openeo_gfmap/preprocessing/__init__.py b/src/openeo_gfmap/preprocessing/__init__.py
index 9f83f02..0f25ca8 100644
--- a/src/openeo_gfmap/preprocessing/__init__.py
+++ b/src/openeo_gfmap/preprocessing/__init__.py
@@ -19,4 +19,4 @@
     "get_bap_score",
     "get_bap_mask",
     "bap_masking",
-]
\ No newline at end of file
+]
diff --git a/src/openeo_gfmap/preprocessing/cloudmasking.py b/src/openeo_gfmap/preprocessing/cloudmasking.py
index 487b705..e605689 100644
--- a/src/openeo_gfmap/preprocessing/cloudmasking.py
+++ b/src/openeo_gfmap/preprocessing/cloudmasking.py
@@ -9,15 +9,16 @@
 SCL_HARMONIZED_NAME: str = "S2-SCL"
 BAPSCORE_HARMONIZED_NAME: str = "S2-BAPSCORE"
 
+
 def mask_scl_dilation(cube: openeo.DataCube, **params: dict) -> openeo.DataCube:
     """Creates a mask from the SCL, dilates it and applies the mask to the optical
     bands of the datacube. The other bands such as DEM, SAR and METEO will not
     be affected by the mask.
     """
     # Asserts if the SCL layer exists
-    assert SCL_HARMONIZED_NAME in cube.metadata.band_names, (
-        f"The SCL band ({SCL_HARMONIZED_NAME}) is not present in the datacube."
-    )
+    assert (
+        SCL_HARMONIZED_NAME in cube.metadata.band_names
+    ), f"The SCL band ({SCL_HARMONIZED_NAME}) is not present in the datacube."
 
     kernel1_size = params.get("kernel1_size", 17)
     kernel2_size = params.get("kernel2_size", 3)
@@ -28,15 +29,11 @@ def mask_scl_dilation(cube: openeo.DataCube, **params: dict) -> openeo.DataCube:
 
     # Only applies the filtering to the optical part of the cube
     optical_cube = cube.filter_bands(
-        bands=list(
-            filter(lambda band: band.startswith("S2"), cube.metadata.band_names)
-        )
+        bands=list(filter(lambda band: band.startswith("S2"), cube.metadata.band_names))
     )
 
     nonoptical_cube = cube.filter_bands(
-        bands=list(
-            filter(lambda band: not band.startswith("S2"), cube.metadata.band_names)
-        )
+        bands=list(filter(lambda band: not band.startswith("S2"), cube.metadata.band_names))
     )
 
     optical_cube = optical_cube.process(
@@ -47,7 +44,7 @@ def mask_scl_dilation(cube: openeo.DataCube, **params: dict) -> openeo.DataCube:
         kernel2_size=kernel2_size,
         mask1_values=[2, 4, 5, 6, 7],
         mask2_values=[3, 8, 9, 10, 11],
-        erosion_kernel_size=erosion_kernel_size
+        erosion_kernel_size=erosion_kernel_size,
     )
 
     if len(nonoptical_cube.metadata.band_names) == 0:
@@ -55,6 +52,7 @@ def mask_scl_dilation(cube: openeo.DataCube, **params: dict) -> openeo.DataCube:
 
     return optical_cube.merge_cubes(nonoptical_cube)
 
+
 def get_bap_score(cube: openeo.DataCube, **params: dict) -> openeo.DataCube:
     """Calculates the Best Available Pixel (BAP) score for the given datacube,
     computed from the SCL layer.
@@ -70,7 +68,7 @@ def get_bap_score(cube: openeo.DataCube, **params: dict) -> openeo.DataCube:
     * Coverage Score: Per date, the percentage of all pixels that are classified
       as a cloud over the entire spatial extent is calculated. The Coverage
       Score is then equal to 1 - the cloud percentage.
-    * Date Score: In order to favor pixels that are observed in the middle of a 
+    * Date Score: In order to favor pixels that are observed in the middle of a
       month, a date score is calculated, which follows a Gaussian shape. I.e.
       the largest scores are given for days in the middle of the month, the
       lowest scores are given for days at the beginning and end of the month.
@@ -113,7 +111,7 @@ def get_bap_score(cube: openeo.DataCube, **params: dict) -> openeo.DataCube:
             kernel2_size=kernel2_size,
             mask1_values=[2, 4, 5, 6, 7],
             mask2_values=[3, 8, 9, 10, 11],
-            erosion_kernel_size=erosion_kernel_size
+            erosion_kernel_size=erosion_kernel_size,
         )
 
     # Replace NaN to 0 to avoid issues in the UDF
@@ -121,15 +119,22 @@ def get_bap_score(cube: openeo.DataCube, **params: dict) -> openeo.DataCube:
 
     score = scl_cube.apply_neighborhood(
         process=openeo.UDF.from_file(str(udf_path)),
-        size=[{'dimension': 'x', 'unit': 'px', 'value': 256}, {'dimension': 'y', 'unit': 'px', 'value': 256}],
-        overlap=[{'dimension': 'x', 'unit': 'px', 'value': 16}, {'dimension': 'y', 'unit': 'px', 'value': 16}],
+        size=[
+            {"dimension": "x", "unit": "px", "value": 256},
+            {"dimension": "y", "unit": "px", "value": 256},
+        ],
+        overlap=[
+            {"dimension": "x", "unit": "px", "value": 16},
+            {"dimension": "y", "unit": "px", "value": 16},
+        ],
     )
 
-    score = score.rename_labels('bands', [BAPSCORE_HARMONIZED_NAME])
+    score = score.rename_labels("bands", [BAPSCORE_HARMONIZED_NAME])
 
     # Merge the score to the scl cube
     return score
 
+
 def get_bap_mask(cube: openeo.DataCube, period: Union[str, list], **params: dict):
     """Computes the bap score and masks the optical bands of the datacube using
     the best scores for each pixel on a given time period. This method both
@@ -155,13 +160,14 @@ def get_bap_mask(cube: openeo.DataCube, period: Union[str, list], **params: dict
         The datacube with the BAP mask applied.
     """
     # Checks if the S2-SCL band is present in the datacube
-    assert SCL_HARMONIZED_NAME in cube.metadata.band_names, (
-        f"The {SCL_HARMONIZED_NAME} band is not present in the datacube."
-    )
+    assert (
+        SCL_HARMONIZED_NAME in cube.metadata.band_names
+    ), f"The {SCL_HARMONIZED_NAME} band is not present in the datacube."
 
     bap_score = get_bap_score(cube, **params)
 
     if isinstance(period, str):
+
         def max_score_selection(score):
             max_score = score.max()
             return score.array_apply(lambda x: x != max_score)
@@ -171,27 +177,26 @@ def max_score_selection(score):
             size=[
                 {"dimension": "x", "unit": "px", "value": 1},
                 {"dimension": "y", "unit": "px", "value": 1},
-                {"dimension": "t", "value": period}
+                {"dimension": "t", "value": period},
             ],
-            overlap=[]
+            overlap=[],
         )
     elif isinstance(period, list):
         udf_path = Path(__file__).parent / "udf_rank.py"
         rank_mask = bap_score.apply_neighborhood(
-            process=openeo.UDF.from_file(
-                str(udf_path),
-                context={"intervals": period}
-            ),
+            process=openeo.UDF.from_file(str(udf_path), context={"intervals": period}),
             size=[
-                {'dimension': 'x', 'unit': 'px', 'value': 256},
-                {'dimension': 'y', 'unit': 'px', 'value': 256}
+                {"dimension": "x", "unit": "px", "value": 256},
+                {"dimension": "y", "unit": "px", "value": 256},
             ],
             overlap=[],
         )
     else:
-        raise ValueError(f"'period' must be a string or a list of dates (in YYYY-mm-dd format), got {period}.")
+        raise ValueError(
+            f"'period' must be a string or a list of dates (in YYYY-mm-dd format), got {period}."
+        )
 
-    return rank_mask.rename_labels('bands', ['S2-BAPMASK'])
+    return rank_mask.rename_labels("bands", ["S2-BAPMASK"])
 
 
 def bap_masking(cube: openeo.DataCube, period: Union[str, list], **params: dict):
@@ -213,22 +218,16 @@ def bap_masking(cube: openeo.DataCube, period: Union[str, list], **params: dict)
         The datacube with the BAP mask applied.
     """
     optical_cube = cube.filter_bands(
-        bands=list(
-            filter(lambda band: band.startswith("S2"), cube.metadata.band_names)
-        )
+        bands=list(filter(lambda band: band.startswith("S2"), cube.metadata.band_names))
     )
 
     nonoptical_cube = cube.filter_bands(
-        bands=list(
-            filter(lambda band: not band.startswith("S2"), cube.metadata.band_names)
-        )
+        bands=list(filter(lambda band: not band.startswith("S2"), cube.metadata.band_names))
     )
 
     rank_mask = get_bap_mask(optical_cube, period, **params)
 
-    optical_cube = optical_cube.mask(
-        rank_mask.resample_cube_spatial(cube)
-    )
+    optical_cube = optical_cube.mask(rank_mask.resample_cube_spatial(cube))
 
     # Do not merge if bands are empty!
     if len(nonoptical_cube.metadata.band_names) == 0:
diff --git a/src/openeo_gfmap/preprocessing/compositing.py b/src/openeo_gfmap/preprocessing/compositing.py
index 5f7eaf4..a0b12f2 100644
--- a/src/openeo_gfmap/preprocessing/compositing.py
+++ b/src/openeo_gfmap/preprocessing/compositing.py
@@ -14,6 +14,7 @@ def median_compositing(cube: openeo.DataCube, period: Union[str, list]) -> opene
     elif isinstance(period, list):
         return cube.aggregate_temporal(intervals=period, reducer="median", dimension="t")
 
+
 def mean_compositing(cube: openeo.DataCube, period: str) -> openeo.DataCube:
     """Perfrom mean compositing on the given datacube."""
     if isinstance(period, str):
diff --git a/src/openeo_gfmap/preprocessing/interpolation.py b/src/openeo_gfmap/preprocessing/interpolation.py
index 79fe256..c404403 100644
--- a/src/openeo_gfmap/preprocessing/interpolation.py
+++ b/src/openeo_gfmap/preprocessing/interpolation.py
@@ -4,8 +4,8 @@
 import openeo
 
 
-def linear_interpolation(cube: openeo.DataCube,) -> openeo.DataCube:
+def linear_interpolation(
+    cube: openeo.DataCube,
+) -> openeo.DataCube:
     """Perform linear interpolation on the given datacube."""
-    return cube.apply_dimension(
-        dimension="t", process="array_interpolate_linear"
-    )
\ No newline at end of file
+    return cube.apply_dimension(dimension="t", process="array_interpolate_linear")
diff --git a/src/openeo_gfmap/preprocessing/udf_rank.py b/src/openeo_gfmap/preprocessing/udf_rank.py
index fa8e679..c3c48b0 100644
--- a/src/openeo_gfmap/preprocessing/udf_rank.py
+++ b/src/openeo_gfmap/preprocessing/udf_rank.py
@@ -16,7 +16,7 @@ def apply_datacube(cube: XarrayDataCube, context: dict) -> XarrayDataCube:
     """
     # First check if the period is defined in the context
     intervals = context.get("intervals", None)
-    array = cube.get_array().transpose('t', 'bands', 'y', 'x')
+    array = cube.get_array().transpose("t", "bands", "y", "x")
 
     bap_score = array.sel(bands="S2-BAPSCORE")
 
@@ -24,7 +24,6 @@ def select_maximum(score: xr.DataArray):
         max_score = score.max(dim="t")
         return score == max_score
 
-
     if isinstance(intervals, str):
         raise NotImplementedError(
             "Period as string is not implemented yet, please provide a list of interval tuples."
@@ -32,10 +31,8 @@ def select_maximum(score: xr.DataArray):
     elif isinstance(intervals, list):
         # Convert YYYY-mm-dd to datetime64 objects
         time_bins = [np.datetime64(interval[0]) for interval in intervals]
-        
-        rank_mask = bap_score.groupby_bins('t', bins=time_bins).map(
-            select_maximum
-        )
+
+        rank_mask = bap_score.groupby_bins("t", bins=time_bins).map(select_maximum)
     else:
         raise ValueError("Period is not defined in the UDF. Cannot run it.")
 
diff --git a/src/openeo_gfmap/preprocessing/udf_score.py b/src/openeo_gfmap/preprocessing/udf_score.py
index 056e2c2..64f20f4 100644
--- a/src/openeo_gfmap/preprocessing/udf_score.py
+++ b/src/openeo_gfmap/preprocessing/udf_score.py
@@ -11,11 +11,12 @@
 
 
 def apply_datacube(cube: XarrayDataCube, context: dict) -> XarrayDataCube:
-
     cube_array: xr.DataArray = cube.get_array()
-    cube_array = cube_array.transpose('t', 'bands', 'y', 'x')
+    cube_array = cube_array.transpose("t", "bands", "y", "x")
 
-    clouds = np.logical_or(np.logical_and(cube_array < 11, cube_array >= 8), cube_array == 3).isel(bands=0)
+    clouds = np.logical_or(np.logical_and(cube_array < 11, cube_array >= 8), cube_array == 3).isel(
+        bands=0
+    )
 
     weights = [1, 0.8, 0.5]
 
@@ -24,13 +25,16 @@ def apply_datacube(cube: XarrayDataCube, context: dict) -> XarrayDataCube:
     sigma = 5
     mu = 15
     score_doy = 1 / (sigma * math.sqrt(2 * math.pi)) * np.exp(-0.5 * ((times - mu) / sigma) ** 2)
-    score_doy = np.broadcast_to(score_doy[:, np.newaxis, np.newaxis],
-                                [cube_array.sizes['t'], cube_array.sizes['y'], cube_array.sizes['x']])
+    score_doy = np.broadcast_to(
+        score_doy[:, np.newaxis, np.newaxis],
+        [cube_array.sizes["t"], cube_array.sizes["y"], cube_array.sizes["x"]],
+    )
 
     # Calculate the Distance To Cloud score
     # Erode
     # Source: https://github.com/dzanaga/satio-pc/blob/e5fc46c0c14bba77e01dca409cf431e7ef22c077/src/satio_pc/preprocessing/clouds.py#L127
     e = footprints.disk(3)
+
     # Define a function to apply binary erosion
     def erode(image, selem):
         return ~binary_erosion(image, selem)
@@ -39,12 +43,12 @@ def erode(image, selem):
     eroded = xr.apply_ufunc(
         erode,  # function to apply
         clouds,  # input DataArray
-        input_core_dims=[['y', 'x']],  # dimensions over which to apply function
-        output_core_dims=[['y', 'x']],  # dimensions of the output
+        input_core_dims=[["y", "x"]],  # dimensions over which to apply function
+        output_core_dims=[["y", "x"]],  # dimensions of the output
         vectorize=True,  # vectorize the function over non-core dimensions
         dask="parallelized",  # enable dask parallelization
         output_dtypes=[np.int32],  # data type of the output
-        kwargs={'selem': e}  # additional keyword arguments to pass to erode
+        kwargs={"selem": e},  # additional keyword arguments to pass to erode
     )
 
     # Distance to cloud = dilation
@@ -53,33 +57,38 @@ def erode(image, selem):
     d = xr.apply_ufunc(
         distance_transform_cdt,
         eroded,
-        input_core_dims=[['y', 'x']],
-        output_core_dims=[['y', 'x']],
+        input_core_dims=[["y", "x"]],
+        output_core_dims=[["y", "x"]],
         vectorize=True,
         dask="parallelized",
-        output_dtypes=[np.int32]
+        output_dtypes=[np.int32],
     )
     d = xr.where(d == -1, d_req, d)
     score_clouds = 1 / (1 + np.exp(-0.2 * (np.minimum(d, d_req) - (d_req - d_min) / 2)))
 
     # Calculate the Coverage score
-    score_cov = 1 - clouds.sum(dim='x').sum(dim='y') / (
-            cube_array.sizes['x'] * cube_array.sizes['y'])
-    score_cov = np.broadcast_to(score_cov.values[:, np.newaxis, np.newaxis],
-                                [cube_array.sizes['t'], cube_array.sizes['y'], cube_array.sizes['x']])
+    score_cov = 1 - clouds.sum(dim="x").sum(dim="y") / (
+        cube_array.sizes["x"] * cube_array.sizes["y"]
+    )
+    score_cov = np.broadcast_to(
+        score_cov.values[:, np.newaxis, np.newaxis],
+        [cube_array.sizes["t"], cube_array.sizes["y"], cube_array.sizes["x"]],
+    )
 
     # Final score is weighted average
-    score = (weights[0] * score_clouds + weights[1] * score_doy + weights[2] * score_cov) / sum(weights)
-    score = np.where(cube_array.values[:,0,:,:]==0, 0, score)
+    score = (weights[0] * score_clouds + weights[1] * score_doy + weights[2] * score_cov) / sum(
+        weights
+    )
+    score = np.where(cube_array.values[:, 0, :, :] == 0, 0, score)
 
     score_da = xr.DataArray(
         score,
         coords={
-            't': cube_array.coords['t'],
-            'y': cube_array.coords['y'],
-            'x': cube_array.coords['x'],
+            "t": cube_array.coords["t"],
+            "y": cube_array.coords["y"],
+            "x": cube_array.coords["x"],
         },
-        dims=['t', 'y', 'x']
+        dims=["t", "y", "x"],
     )
 
     score_da = score_da.expand_dims(
@@ -88,6 +97,6 @@ def erode(image, selem):
         },
     )
 
-    score_da = score_da.transpose('t', 'bands', 'y', 'x')
+    score_da = score_da.transpose("t", "bands", "y", "x")
 
     return XarrayDataCube(score_da)
diff --git a/src/openeo_gfmap/spatial.py b/src/openeo_gfmap/spatial.py
index 454ef91..2bc0caf 100644
--- a/src/openeo_gfmap/spatial.py
+++ b/src/openeo_gfmap/spatial.py
@@ -29,13 +29,15 @@ def __dict__(self):
         }
 
     def __iter__(self):
-        return iter([
-            ("west", self.west),
-            ("south", self.south),
-            ("east", self.east),
-            ("north", self.north),
-            ("crs", self.epsg)
-        ])
+        return iter(
+            [
+                ("west", self.west),
+                ("south", self.south),
+                ("east", self.east),
+                ("north", self.north),
+                ("crs", self.epsg),
+            ]
+        )
 
 
 SpatialContext = Union[GeoJSON, BoundingBoxExtent]
diff --git a/src/openeo_gfmap/utils/__init__.py b/src/openeo_gfmap/utils/__init__.py
index b9f6608..9b9efd3 100644
--- a/src/openeo_gfmap/utils/__init__.py
+++ b/src/openeo_gfmap/utils/__init__.py
@@ -11,6 +11,11 @@
 )
 
 __all__ = [
-    "load_json", "normalize_array", "select_optical_bands", "array_bounds",
-    "select_sar_bands", "arrays_cosine_similarity", "quintad_intervals"
+    "load_json",
+    "normalize_array",
+    "select_optical_bands",
+    "array_bounds",
+    "select_sar_bands",
+    "arrays_cosine_similarity",
+    "quintad_intervals",
 ]
diff --git a/src/openeo_gfmap/utils/catalogue.py b/src/openeo_gfmap/utils/catalogue.py
index b552da4..c4b1d8d 100644
--- a/src/openeo_gfmap/utils/catalogue.py
+++ b/src/openeo_gfmap/utils/catalogue.py
@@ -10,7 +10,7 @@ def _check_cdse_catalogue(
     collection: str,
     spatial_extent: SpatialContext,
     temporal_extent: TemporalContext,
-    **additional_parameters: dict
+    **additional_parameters: dict,
 ) -> bool:
     """Checks if there is at least one product available in the
     given spatio-temporal context for a collection in the CDSE catalogue,
@@ -37,16 +37,23 @@ def _check_cdse_catalogue(
         # Transform geojson into shapely geometry and compute bounds
         bounds = shape(spatial_extent).bounds
     elif isinstance(spatial_extent, SpatialContext):
-        bounds = [spatial_extent.west, spatial_extent.south, spatial_extent.east, spatial_extent.north]
+        bounds = [
+            spatial_extent.west,
+            spatial_extent.south,
+            spatial_extent.east,
+            spatial_extent.north,
+        ]
     else:
-        raise ValueError('Provided spatial extent is not a valid GeoJSON or SpatialContext object.')
+        raise ValueError(
+            "Provided spatial extent is not a valid GeoJSON or SpatialContext object."
+        )
 
     minx, miny, maxx, maxy = bounds
 
     # The date format should be YYYY-MM-DD
-    start_date = f'{temporal_extent.start_date}T00:00:00Z'
-    end_date = f'{temporal_extent.end_date}T00:00:00Z'
-    
+    start_date = f"{temporal_extent.start_date}T00:00:00Z"
+    end_date = f"{temporal_extent.end_date}T00:00:00Z"
+
     url = (
         f"https://catalogue.dataspace.copernicus.eu/resto/api/collections/"
         f"{collection}/search.json?box={minx},{miny},{maxx},{maxy}"
@@ -66,7 +73,9 @@ def _check_cdse_catalogue(
 
     body = response.json()
     grd_tiles = list(
-        filter(lambda feature: feature["properties"]["productType"].contains("GRD"), body["features"])
+        filter(
+            lambda feature: feature["properties"]["productType"].contains("GRD"), body["features"]
+        )
     )
 
-    return len(grd_tiles) > 0 
+    return len(grd_tiles) > 0
diff --git a/src/openeo_gfmap/utils/intervals.py b/src/openeo_gfmap/utils/intervals.py
index f53bc38..6c9f3ba 100644
--- a/src/openeo_gfmap/utils/intervals.py
+++ b/src/openeo_gfmap/utils/intervals.py
@@ -1,4 +1,4 @@
-"""Utilitary function for intervals, useful for temporal aggregation 
+"""Utilitary function for intervals, useful for temporal aggregation
 methods.
 """
 
@@ -8,14 +8,14 @@
 
 
 def quintad_intervals(temporal_extent: TemporalContext) -> list:
-    """ Returns a list of tuples (start_date, end_date) of quintad intervals
-        from the input temporal extent. Quintad intervals are intervals of
-        generally 5 days, that never overlap two months.
+    """Returns a list of tuples (start_date, end_date) of quintad intervals
+    from the input temporal extent. Quintad intervals are intervals of
+    generally 5 days, that never overlap two months.
 
-        All months are divided in 6 quintads, where the 6th quintad might
-        contain 6 days for months of 31 days.
-        For the month of February, the 6th quintad is only of three days, or
-        four days for the leap year.
+    All months are divided in 6 quintads, where the 6th quintad might
+    contain 6 days for months of 31 days.
+    For the month of February, the 6th quintad is only of three days, or
+    four days for the leap year.
     """
     start_date, end_date = temporal_extent.to_datetime()
     quintads = []
@@ -27,7 +27,7 @@ def quintad_intervals(temporal_extent: TemporalContext) -> list:
         offset = (start_date - timedelta(days=1)).day % 5
         current_date = current_date - timedelta(days=offset)
     else:
-        offset = 0    
+        offset = 0
 
     while current_date <= end_date:
         # Get the last day of the current month
diff --git a/src/openeo_gfmap/utils/tile_processing.py b/src/openeo_gfmap/utils/tile_processing.py
index c7add8a..fa1efab 100644
--- a/src/openeo_gfmap/utils/tile_processing.py
+++ b/src/openeo_gfmap/utils/tile_processing.py
@@ -4,63 +4,49 @@
 import xarray as xr
 
 
-def normalize_array(
-    inarr: xr.DataArray,
-    percentile: float = 0.99
-) -> xr.DataArray:
-    """ Performs normalization between 0.0 and 1.0 using the given
+def normalize_array(inarr: xr.DataArray, percentile: float = 0.99) -> xr.DataArray:
+    """Performs normalization between 0.0 and 1.0 using the given
     percentile.
     """
     quantile_value = inarr.quantile(percentile, dim=["x", "y", "t"])
     minimum = inarr.min(dim=["x", "y", "t"])
 
     inarr = (inarr - minimum) / (quantile_value - minimum)
-    
+
     # Perform clipping on values that are higher than the computed quantile
     return inarr.where(inarr < 1.0, 1.0)
 
-def select_optical_bands(
-    inarr: xr.DataArray
-) -> xr.DataArray:
+
+def select_optical_bands(inarr: xr.DataArray) -> xr.DataArray:
     """Filters and keep only the optical bands for a given array."""
     return inarr.sel(
-        bands=[
-            band
-            for band in inarr.coords["bands"].to_numpy()
-            if band.startswith("S2-B")
-        ]
+        bands=[band for band in inarr.coords["bands"].to_numpy() if band.startswith("S2-B")]
     )
 
-def select_sar_bands(
-    inarr: xr.DataArray
-) -> xr.DataArray:
+
+def select_sar_bands(inarr: xr.DataArray) -> xr.DataArray:
     """Filters and keep only the SAR bands for a given array."""
     return inarr.sel(
         bands=[
-            band
-            for band in inarr.coords["bands"].to_numpy()
-            if band in ["VV", "VH", "HH", "HV"]
+            band for band in inarr.coords["bands"].to_numpy() if band in ["VV", "VH", "HH", "HV"]
         ]
     )
 
-def array_bounds(
-    inarr: xr.DataArray
-) -> tuple:
+
+def array_bounds(inarr: xr.DataArray) -> tuple:
     """Returns the 4 bounds values for the x and y coordinates of the tile"""
     return (
         inarr.coords["x"].min().item(),
         inarr.coords["y"].min().item(),
         inarr.coords["x"].max().item(),
-        inarr.coords["y"].max().item()
+        inarr.coords["y"].max().item(),
     )
 
-def arrays_cosine_similarity(
-    first_array: xr.DataArray,
-    second_array: xr.DataArray
-) -> float:
+
+def arrays_cosine_similarity(first_array: xr.DataArray, second_array: xr.DataArray) -> float:
     """Returns a similarity score based on normalized cosine distance. The
     input arrays must have similar ranges to obtain a valid score.
-    1.0 represents the best score (same tiles), while 0.0 is the worst score. 
+    1.0 represents the best score (same tiles), while 0.0 is the worst score.
     """
     dot_product = np.sum(first_array * second_array)
     first_norm = np.linalg.norm(first_array)
@@ -68,4 +54,3 @@ def arrays_cosine_similarity(
     similarity = (dot_product / (first_norm * second_norm)).item()
 
     return similarity
-
diff --git a/tests/test_openeo_gfmap/test_cloud_masking.py b/tests/test_openeo_gfmap/test_cloud_masking.py
index 6963319..649314d 100644
--- a/tests/test_openeo_gfmap/test_cloud_masking.py
+++ b/tests/test_openeo_gfmap/test_cloud_masking.py
@@ -23,13 +23,12 @@
     south=51.215806593713,
     east=5.060320484557499,
     north=51.22149744530769,
-    epsg=4326
+    epsg=4326,
 )
 
 # November 2022 to February 2023
-temporal_extent = TemporalContext(
-    start_date="2022-11-01", end_date="2023-02-28"
-)
+temporal_extent = TemporalContext(start_date="2022-11-01", end_date="2023-02-28")
+
 
 @pytest.mark.parametrize("backend", backends)
 def test_bap_score(backend: Backend):
@@ -37,33 +36,25 @@ def test_bap_score(backend: Backend):
     backend_context = BackendContext(backend=backend)
 
     # Additional parameters
-    fetching_parameters = {
-        "fetching_resolution": 10.0
-    }
+    fetching_parameters = {"fetching_resolution": 10.0}
 
-    preprocessing_parameters = {
-        "apply_scl_dilation": True
-    }
+    preprocessing_parameters = {"apply_scl_dilation": True}
 
     # Fetch the datacube
     s2_extractor = build_sentinel2_l2a_extractor(
         backend_context=backend_context,
         bands=["S2-B04", "S2-B08", "S2-SCL"],
         fetch_type=FetchType.TILE,
-        **fetching_parameters
+        **fetching_parameters,
     )
 
-    cube = s2_extractor.get_cube(
-        connection, spatial_extent, temporal_extent
-    )
+    cube = s2_extractor.get_cube(connection, spatial_extent, temporal_extent)
 
     # Compute the BAP score
     bap_score = get_bap_score(cube, **preprocessing_parameters)
     ndvi = cube.ndvi(nir="S2-B08", red="S2-B04")
-    
-    cube = bap_score.merge_cubes(ndvi).rename_labels(
-        'bands', ['S2-BAPSCORE', 'S2-NDVI']
-    )
+
+    cube = bap_score.merge_cubes(ndvi).rename_labels("bands", ["S2-BAPSCORE", "S2-NDVI"])
 
     job = cube.create_job(
         title="BAP score unittest",
@@ -71,12 +62,11 @@ def test_bap_score(backend: Backend):
     )
 
     job.start_and_wait()
-    
+
     for asset in job.get_results().get_assets():
         if asset.metadata["type"].startswith("application/x-netcdf"):
-            asset.download(
-                Path(__file__).parent / f"results/bap_score_{backend.value}.nc"
-            )
+            asset.download(Path(__file__).parent / f"results/bap_score_{backend.value}.nc")
+
 
 @pytest.mark.parametrize("backend", backends)
 def test_bap_masking(backend: Backend):
@@ -84,25 +74,21 @@ def test_bap_masking(backend: Backend):
     backend_context = BackendContext(backend=backend)
 
     # Additional parameters
-    fetching_parameters = {
-        "fetching_resolution": 10.0
-    }
+    fetching_parameters = {"fetching_resolution": 10.0}
 
     # Fetch the datacube
     s2_extractor = build_sentinel2_l2a_extractor(
         backend_context=backend_context,
         bands=["S2-B04", "S2-B03", "S2-B02", "S2-SCL"],
         fetch_type=FetchType.TILE,
-        **fetching_parameters
+        **fetching_parameters,
     )
 
-    cube = s2_extractor.get_cube(
-        connection, spatial_extent, temporal_extent
-    )
+    cube = s2_extractor.get_cube(connection, spatial_extent, temporal_extent)
 
     cube = cube.linear_scale_range(0, 65535, 0, 65535)
 
-    # Perform masking with BAP, masking optical bands 
+    # Perform masking with BAP, masking optical bands
     cube = bap_masking(cube, period="dekad")
 
     # Perform compositing, the cube should be only composed of optical bands
@@ -111,9 +97,7 @@ def test_bap_masking(backend: Backend):
     cube = cube.linear_scale_range(0, 65535, 0, 65535)
 
     # Remove SCL
-    cube = cube.filter_bands(
-        [band for band in cube.metadata.band_names if band != "S2-SCL"]
-    )
+    cube = cube.filter_bands([band for band in cube.metadata.band_names if band != "S2-SCL"])
 
     job = cube.create_job(
         title="BAP compositing unittest",
@@ -124,9 +108,8 @@ def test_bap_masking(backend: Backend):
 
     for asset in job.get_results().get_assets():
         if asset.metadata["type"].startswith("application/x-netcdf"):
-            asset.download(
-                Path(__file__).parent / f"results/bap_composited_{backend.value}.nc"
-            )
+            asset.download(Path(__file__).parent / f"results/bap_composited_{backend.value}.nc")
+
 
 @pytest.mark.parametrize("backend", backends)
 def test_bap_quintad(backend: Backend):
@@ -134,54 +117,46 @@ def test_bap_quintad(backend: Backend):
     backend_context = BackendContext(backend=backend)
 
     # Additional parameters
-    fetching_parameters = {
-        "fetching_resolution": 10.0
-    }
-    preprocessing_parameters = {
-        "apply_scl_dilation": True 
-    }
+    fetching_parameters = {"fetching_resolution": 10.0}
+    preprocessing_parameters = {"apply_scl_dilation": True}
 
     # Fetch the datacube
     s2_extractor = build_sentinel2_l2a_extractor(
         backend_context=backend_context,
         bands=["S2-SCL"],
         fetch_type=FetchType.TILE,
-        **fetching_parameters
+        **fetching_parameters,
     )
 
-    cube = s2_extractor.get_cube(
-        connection, spatial_extent, temporal_extent
-    )
+    cube = s2_extractor.get_cube(connection, spatial_extent, temporal_extent)
 
-    compositing_intervals = quintad_intervals(
-        temporal_extent
-    )
+    compositing_intervals = quintad_intervals(temporal_extent)
 
     expected_intervals = [
-        ('2022-11-01', '2022-11-05'),
-        ('2022-11-06', '2022-11-10'),
-        ('2022-11-11', '2022-11-15'),
-        ('2022-11-16', '2022-11-20'),
-        ('2022-11-21', '2022-11-25'),
-        ('2022-11-26', '2022-11-30'),
-        ('2022-12-01', '2022-12-05'),
-        ('2022-12-06', '2022-12-10'),
-        ('2022-12-11', '2022-12-15'),
-        ('2022-12-16', '2022-12-20'),
-        ('2022-12-21', '2022-12-25'),
-        ('2022-12-26', '2022-12-31'),
-        ('2023-01-01', '2023-01-05'),
-        ('2023-01-06', '2023-01-10'),
-        ('2023-01-11', '2023-01-15'),
-        ('2023-01-16', '2023-01-20'),
-        ('2023-01-21', '2023-01-25'),
-        ('2023-01-26', '2023-01-31'),
-        ('2023-02-01', '2023-02-05'),
-        ('2023-02-06', '2023-02-10'),
-        ('2023-02-11', '2023-02-15'),
-        ('2023-02-16', '2023-02-20'),
-        ('2023-02-21', '2023-02-25'),
-        ('2023-02-26', '2023-02-28'),
+        ("2022-11-01", "2022-11-05"),
+        ("2022-11-06", "2022-11-10"),
+        ("2022-11-11", "2022-11-15"),
+        ("2022-11-16", "2022-11-20"),
+        ("2022-11-21", "2022-11-25"),
+        ("2022-11-26", "2022-11-30"),
+        ("2022-12-01", "2022-12-05"),
+        ("2022-12-06", "2022-12-10"),
+        ("2022-12-11", "2022-12-15"),
+        ("2022-12-16", "2022-12-20"),
+        ("2022-12-21", "2022-12-25"),
+        ("2022-12-26", "2022-12-31"),
+        ("2023-01-01", "2023-01-05"),
+        ("2023-01-06", "2023-01-10"),
+        ("2023-01-11", "2023-01-15"),
+        ("2023-01-16", "2023-01-20"),
+        ("2023-01-21", "2023-01-25"),
+        ("2023-01-26", "2023-01-31"),
+        ("2023-02-01", "2023-02-05"),
+        ("2023-02-06", "2023-02-10"),
+        ("2023-02-11", "2023-02-15"),
+        ("2023-02-16", "2023-02-20"),
+        ("2023-02-21", "2023-02-25"),
+        ("2023-02-26", "2023-02-28"),
     ]
 
     assert compositing_intervals == expected_intervals
@@ -192,20 +167,18 @@ def test_bap_quintad(backend: Backend):
     # Create a new extractor for the whole data now
     fetching_parameters = {
         "fetching_resolution": 10.0,
-        "pre_mask": bap_mask  # Use of the pre-computed bap mask to load inteligently the data
+        "pre_mask": bap_mask,  # Use of the pre-computed bap mask to load inteligently the data
     }
 
     s2_extractor = build_sentinel2_l2a_extractor(
         backend_context=backend_context,
         bands=["S2-B04", "S2-B03", "S2-B02", "S2-B08", "S2-SCL"],
         fetch_type=FetchType.TILE,
-        **fetching_parameters
+        **fetching_parameters,
     )
 
     # Performs quintal compositing
-    cube = s2_extractor.get_cube(
-        connection, spatial_extent, temporal_extent
-    )
+    cube = s2_extractor.get_cube(connection, spatial_extent, temporal_extent)
 
     cube = median_compositing(cube, period=compositing_intervals)
 
@@ -220,6 +193,4 @@ def test_bap_quintad(backend: Backend):
 
     for asset in job.get_results().get_assets():
         if asset.metadata["type"].startswith("application/x-netcdf"):
-            asset.download(
-                Path(__file__).parent / f"results/bap_quintad_{backend.value}.nc"
-            )
+            asset.download(Path(__file__).parent / f"results/bap_quintad_{backend.value}.nc")
diff --git a/tests/test_openeo_gfmap/test_intervals.py b/tests/test_openeo_gfmap/test_intervals.py
index 2ce833d..876d0d0 100644
--- a/tests/test_openeo_gfmap/test_intervals.py
+++ b/tests/test_openeo_gfmap/test_intervals.py
@@ -19,6 +19,7 @@ def test_quintad_january():
 
     assert quintad_intervals(temporal_extent) == expected
 
+
 def test_quintad_april():
     start_date = "2023-04-01"
     end_date = "2023-04-30"
@@ -36,6 +37,7 @@ def test_quintad_april():
 
     assert quintad_intervals(temporal_extent) == expected
 
+
 def test_quintad_february_nonleap():
     start_date = "2023-02-01"
     end_date = "2023-02-28"
@@ -53,6 +55,7 @@ def test_quintad_february_nonleap():
 
     assert quintad_intervals(temporal_extent) == expected
 
+
 def test_quitad_february_leapyear():
     start_date = "2024-02-01"
     end_date = "2024-02-29"
@@ -70,12 +73,13 @@ def test_quitad_february_leapyear():
 
     assert quintad_intervals(temporal_extent) == expected
 
+
 def test_quintad_four_months():
     start_date = "2023-01-01"
     end_date = "2023-04-30"
 
     temporal_extent = TemporalContext(start_date, end_date)
-    
+
     expected = [
         ("2023-01-01", "2023-01-05"),
         ("2023-01-06", "2023-01-10"),
@@ -105,6 +109,7 @@ def test_quintad_four_months():
 
     assert quintad_intervals(temporal_extent) == expected
 
+
 def test_quintad_july_august():
     start_date = "2023-07-01"
     end_date = "2023-08-31"
@@ -128,6 +133,7 @@ def test_quintad_july_august():
 
     assert quintad_intervals(temporal_extent) == expected
 
+
 def test_quintad_mid_month():
     start_date = "2023-01-02"
     end_date = "2023-01-31"
@@ -145,6 +151,7 @@ def test_quintad_mid_month():
 
     assert quintad_intervals(temporal_extent) == expected
 
+
 def test_quintad_full_year():
     # non-leap year
     start_date = "2023-01-01"
@@ -162,6 +169,7 @@ def test_quintad_full_year():
 
     assert len(quintad_intervals(temporal_extent)) == 72
 
+
 def test_quintad_mid_month_february():
     start_date = "2024-01-31"
     end_date = "2024-03-02"
@@ -181,6 +189,7 @@ def test_quintad_mid_month_february():
 
     assert quintad_intervals(temporal_extent) == expected
 
+
 def test_quintad_single_day():
     start_date = "2024-02-29"
     end_date = "2024-02-29"
@@ -193,6 +202,7 @@ def test_quintad_single_day():
 
     assert quintad_intervals(temporal_extent) == expected
 
+
 def test_quintad_end_month():
     start_date = "2024-02-14"
     end_date = "2024-03-01"
@@ -209,6 +219,7 @@ def test_quintad_end_month():
 
     assert quintad_intervals(temporal_extent) == expected
 
+
 def test_quintad_new_year():
     start_date = "2023-12-04"
     end_date = "2024-01-01"
@@ -227,4 +238,4 @@ def test_quintad_new_year():
 
     print(quintad_intervals(temporal_extent))
 
-    assert quintad_intervals(temporal_extent) == expected
\ No newline at end of file
+    assert quintad_intervals(temporal_extent) == expected
diff --git a/tests/test_openeo_gfmap/test_s1_fetchers.py b/tests/test_openeo_gfmap/test_s1_fetchers.py
index 9d81fb2..60e483c 100644
--- a/tests/test_openeo_gfmap/test_s1_fetchers.py
+++ b/tests/test_openeo_gfmap/test_s1_fetchers.py
@@ -42,7 +42,7 @@ def sentinel1_grd(
         spatial_extent: SpatialContext,
         temporal_extent: TemporalContext,
         backend: Backend,
-        connection=openeo.Connection
+        connection=openeo.Connection,
     ):
         context = BackendContext(backend)
         country = spatial_extent["country"]
@@ -55,7 +55,7 @@ def sentinel1_grd(
             "coefficient": "gamma0-ellipsoid",
             "load_collection": {
                 "polarization": lambda polar: (polar == "VV") or (polar == "VH"),
-            }
+            },
         }
 
         extractor: CollectionFetcher = build_sentinel1_grd_extractor(
@@ -74,18 +74,11 @@ def sentinel1_grd(
             start_date=temporal_extent[0], end_date=temporal_extent[1]
         )
 
-        cube = extractor.get_cube(
-            connection, spatial_extent, temporal_extent
-        )
+        cube = extractor.get_cube(connection, spatial_extent, temporal_extent)
 
-        output_file = (
-            Path(__file__).parent
-            / f"results/{country}_{backend.value}_sentinel1_grd.nc"
-        )
+        output_file = Path(__file__).parent / f"results/{country}_{backend.value}_sentinel1_grd.nc"
 
-        job = cube.create_job(
-            title="Sentinel1 GRD tile extraction", out_format="NetCDF"
-        )
+        job = cube.create_job(title="Sentinel1 GRD tile extraction", out_format="NetCDF")
 
         job.start_and_wait()
 
@@ -115,8 +108,7 @@ def compare_sentinel1_tiles():
             loaded_tiles = []
             for backend in backend_types:
                 tile_path = (
-                    Path(__file__).parent
-                    / f"results/{country}_{backend.value}_sentinel1_grd.nc"
+                    Path(__file__).parent / f"results/{country}_{backend.value}_sentinel1_grd.nc"
                 )
                 loaded_tiles.append(xr.open_dataset(tile_path, engine="h5netcdf"))
 
@@ -131,7 +123,7 @@ def compare_sentinel1_tiles():
                     dtype = array.dtype
                 else:
                     assert dtype == array.dtype
-            
+
             bounds = None
             for tile in loaded_tiles:
                 tile_bounds = array_bounds(tile)
@@ -147,16 +139,14 @@ def compare_sentinel1_tiles():
             first_tile = normalized_tiles[0]
             for tile_idx in range(1, len(normalized_tiles)):
                 tile_to_compare = normalized_tiles[tile_idx]
-                similarity_score = arrays_cosine_similarity(
-                    first_tile, tile_to_compare
-                )
+                similarity_score = arrays_cosine_similarity(first_tile, tile_to_compare)
                 assert similarity_score >= 0.95
 
     def sentinel1_grd_point_based(
         spatial_context: SpatialContext,
         temporal_context: TemporalContext,
         backend: Backend,
-        connection: openeo.Connection
+        connection: openeo.Connection,
     ):
         """Test the point based extraction from the spatial aggregation of the
         given polygons.
@@ -173,24 +163,17 @@ def sentinel1_grd_point_based(
             "coefficient": "gamma0-ellipsoid",
             "load_collection": {
                 "polarization": lambda polar: (polar == "VV") or (polar == "VH"),
-            }
+            },
         }
         extractor = build_sentinel1_grd_extractor(
-            backend_context=context,
-            bands=bands,
-            fetch_type=FetchType.POINT,
-            **fetching_parameters
+            backend_context=context, bands=bands, fetch_type=FetchType.POINT, **fetching_parameters
         )
 
-        cube = extractor.get_cube(
-            connection, spatial_context, temporal_context
-        )
+        cube = extractor.get_cube(connection, spatial_context, temporal_context)
 
         cube = cube.aggregate_spatial(spatial_context, reducer="mean")
 
-        output_file = (
-            Path(__file__).parent / f"results/points_{backend.value}_sentinel1_grd.nc"
-        )
+        output_file = Path(__file__).parent / f"results/points_{backend.value}_sentinel1_grd.nc"
 
         cube.download(output_file, format="JSON")
 
@@ -202,7 +185,7 @@ def sentinel1_grd_point_based(
                 if band in col:
                     exists = True
             assert exists, f"Couldn't find a single column for band {band}"
-        
+
         assert len(df.columns) % len(bands) == 0, (
             f"The number of columns ({len(df.columns)}) should be a multiple"
             f"of the number of bands ({len(bands)})"
@@ -214,7 +197,7 @@ def sentinel1_grd_polygon_based(
         spatial_context: SpatialContext,
         temporal_context: TemporalContext,
         backend: Backend,
-        connection: openeo.Connection
+        connection: openeo.Connection,
     ):
         context = BackendContext(backend)
         bands = ["S1-VV", "S1-VH"]
@@ -226,14 +209,14 @@ def sentinel1_grd_polygon_based(
             "coefficient": "gamma0-ellipsoid",
             "load_collection": {
                 "polarization": lambda polar: (polar == "VV") or (polar == "VH"),
-            }
+            },
         }
 
         extractor = build_sentinel1_grd_extractor(
             backend_context=context,
             bands=bands,
             fetch_type=FetchType.POLYGON,
-            **fetching_parameters
+            **fetching_parameters,
         )
 
         cube = extractor.get_cube(connection, spatial_context, temporal_context)
@@ -251,30 +234,24 @@ def sentinel1_grd_polygon_based(
         results.download_files(output_folder)
 
         # List all the files available in the folder
-        extracted_files = list(
-            filter(lambda file: file.suffix == ".nc", output_folder.iterdir())
-        )
+        extracted_files = list(filter(lambda file: file.suffix == ".nc", output_folder.iterdir()))
         # Check if there is one file for each polygon
         assert len(extracted_files) == len(spatial_context["features"])
 
 
-
-@pytest.mark.parametrize(
-    "spatial_context, temporal_context, backend", test_configurations
-)
+@pytest.mark.parametrize("spatial_context, temporal_context, backend", test_configurations)
 def test_sentinel1_grd(
-    spatial_context: SpatialContext, temporal_context: TemporalContext,
-    backend: Backend
+    spatial_context: SpatialContext, temporal_context: TemporalContext, backend: Backend
 ):
     connection = BACKEND_CONNECTIONS[backend]()
-    TestS1Extractors.sentinel1_grd(
-        spatial_context, temporal_context, backend, connection
-    )
+    TestS1Extractors.sentinel1_grd(spatial_context, temporal_context, backend, connection)
+
 
 @pytest.mark.depends(on=["test_sentinel1_grd"])
 def test_compare_sentinel1_tiles():
     TestS1Extractors.compare_sentinel1_tiles()
 
+
 @pytest.mark.parametrize("backend", test_backends)
 def test_sentinel1_grd_point_based(backend: Backend):
     connection = BACKEND_CONNECTIONS[backend]()
diff --git a/tests/test_openeo_gfmap/test_s2_fetchers.py b/tests/test_openeo_gfmap/test_s2_fetchers.py
index 722c05b..2acaa06 100644
--- a/tests/test_openeo_gfmap/test_s2_fetchers.py
+++ b/tests/test_openeo_gfmap/test_s2_fetchers.py
@@ -52,21 +52,17 @@
 TEMPORAL_EXTENT_2 = ["2023-01-01", "2023-02-01"]
 
 # Dataset of polygons for POINT based extraction
-POINT_EXTRACTION_DF = (
-    Path(__file__).parent / "resources/malawi_extraction_polygons.gpkg"
-)
+POINT_EXTRACTION_DF = Path(__file__).parent / "resources/malawi_extraction_polygons.gpkg"
 
 # Datase of polygons for Polygon based extraction
-POLYGON_EXTRACTION_DF = (
-    Path(__file__).parent / "resources/puglia_extraction_polygons.gpkg"
-)
+POLYGON_EXTRACTION_DF = Path(__file__).parent / "resources/puglia_extraction_polygons.gpkg"
 
-#test_backends = [Backend.TERRASCOPE, Backend.CDSE]
+# test_backends = [Backend.TERRASCOPE, Backend.CDSE]
 test_backends = [Backend.CDSE]
 
 test_spatio_temporal_extends = [
     (SPATIAL_EXTENT_1, TEMPORAL_EXTENT_1),
-#    (SPATIAL_EXTENT_2, TEMPORAL_EXTENT_2),
+    #    (SPATIAL_EXTENT_2, TEMPORAL_EXTENT_2),
 ]
 
 test_configurations = [
@@ -101,10 +97,9 @@ def sentinel2_l2a(
             "S2-SCL",
             "S2-AOT",
         ]
-        fetching_parameters = {
-            "target_resolution": 10.0,
-            "target_crs": 3035
-        } if country == "Belgium" else {}
+        fetching_parameters = (
+            {"target_resolution": 10.0, "target_crs": 3035} if country == "Belgium" else {}
+        )
         extractor: CollectionFetcher = build_sentinel2_l2a_extractor(
             context=context,
             bands=bands,
@@ -117,7 +112,7 @@ def sentinel2_l2a(
             south=spatial_extent["south"],
             east=spatial_extent["east"],
             north=spatial_extent["north"],
-            epsg=spatial_extent["crs"]
+            epsg=spatial_extent["crs"],
         )
 
         temporal_extent = TemporalContext(
@@ -126,10 +121,7 @@ def sentinel2_l2a(
 
         cube = extractor.get_cube(connection, spatial_extent, temporal_extent)
 
-        output_file = (
-            Path(__file__).parent
-            / f"results/{country}_{backend.value}_sentinel2_l2a.nc"
-        )
+        output_file = Path(__file__).parent / f"results/{country}_{backend.value}_sentinel2_l2a.nc"
 
         cube.download(output_file, format="NetCDF")
 
@@ -154,8 +146,7 @@ def compare_sentinel2_tiles():
                 if backend == Backend.EODC:  # TODO fix EDOC backend first
                     continue
                 tile_path = (
-                    Path(__file__).parent
-                    / f"results/{country}_{backend.value}_sentinel2_l2a.nc"
+                    Path(__file__).parent / f"results/{country}_{backend.value}_sentinel2_l2a.nc"
                 )
                 loaded_tiles.append(xr.open_dataset(tile_path, engine="h5netcdf"))
 
@@ -188,9 +179,7 @@ def compare_sentinel2_tiles():
             first_tile = normalized_tiles[0]
             for tile_idx in range(1, len(normalized_tiles)):
                 tile_to_compare = normalized_tiles[tile_idx]
-                similarity_score = arrays_cosine_similarity(
-                    first_tile, tile_to_compare
-                )
+                similarity_score = arrays_cosine_similarity(first_tile, tile_to_compare)
                 assert similarity_score >= 0.95
 
     def sentinel2_l2a_point_based(
@@ -219,9 +208,7 @@ def sentinel2_l2a_point_based(
 
         cube = cube.aggregate_spatial(spatial_context, reducer="mean")
 
-        output_file = (
-            Path(__file__).parent / f"results/points_{backend.value}_sentinel2_l2a.json"
-        )
+        output_file = Path(__file__).parent / f"results/points_{backend.value}_sentinel2_l2a.json"
 
         cube.download(output_file, format="JSON")
 
@@ -274,29 +261,24 @@ def sentinel2_l2a_polygon_based(
         results.download_files(output_folder)
 
         # List all the files available in the folder
-        extracted_files = list(
-            filter(lambda file: file.suffix == ".nc", output_folder.iterdir())
-        )
+        extracted_files = list(filter(lambda file: file.suffix == ".nc", output_folder.iterdir()))
         # Check if there is one file for each polygon
         assert len(extracted_files) == len(spatial_context["features"])
 
 
-@pytest.mark.parametrize(
-    "spatial_context, temporal_context, backend", test_configurations
-)
+@pytest.mark.parametrize("spatial_context, temporal_context, backend", test_configurations)
 def test_sentinel2_l2a(
     spatial_context: SpatialContext, temporal_context: TemporalContext, backend: Backend
 ):
     connection = BACKEND_CONNECTIONS[backend]()
-    TestS2Extractors.sentinel2_l2a(
-        spatial_context, temporal_context, backend, connection
-    )
+    TestS2Extractors.sentinel2_l2a(spatial_context, temporal_context, backend, connection)
 
 
 @pytest.mark.depends(on=["test_sentinel2_l2a"])
 def test_compare_sentinel2_tiles():
     TestS2Extractors.compare_sentinel2_tiles()
 
+
 @pytest.mark.parametrize("backend", test_backends)
 def test_sentinel2_l2a_point_based(backend: Backend):
     connection = BACKEND_CONNECTIONS[backend]()