diff --git a/eis_toolkit/cli.py b/eis_toolkit/cli.py
index dc5c5c60..c709aff5 100644
--- a/eis_toolkit/cli.py
+++ b/eis_toolkit/cli.py
@@ -1275,7 +1275,9 @@ def distance_to_anomaly_cli(
 
     Uses only the first band of the raster.
     """
-    from eis_toolkit.raster_processing.distance_to_anomaly import distance_to_anomaly
+    from sys import platform
+
+    from eis_toolkit.raster_processing.distance_to_anomaly import distance_to_anomaly, distance_to_anomaly_gdal
 
     typer.echo("Progress: 10%")
 
@@ -1286,13 +1288,23 @@ def distance_to_anomaly_cli(
 
     with rasterio.open(input_raster) as raster:
         typer.echo("Progress: 25%")
-        out_image, out_meta = distance_to_anomaly(
-            anomaly_raster_profile=raster.profile,
-            anomaly_raster_data=raster.read(1),
-            threshold_criteria_value=threshold_criteria_value,
-            threshold_criteria=get_enum_values(threshold_criteria),
-            max_distance=max_distance,
-        )
+        # Use optimized version if Windows
+        if platform == "win32":
+            out_image, out_meta = distance_to_anomaly_gdal(
+                anomaly_raster_profile=raster.profile,
+                anomaly_raster_data=raster.read(1),
+                threshold_criteria_value=threshold_criteria_value,
+                threshold_criteria=get_enum_values(threshold_criteria),
+                max_distance=max_distance,
+            )
+        else:
+            out_image, out_meta = distance_to_anomaly(
+                anomaly_raster_profile=raster.profile,
+                anomaly_raster_data=raster.read(1),
+                threshold_criteria_value=threshold_criteria_value,
+                threshold_criteria=get_enum_values(threshold_criteria),
+                max_distance=max_distance,
+            )
 
     typer.echo("Progress: 75%")
 
diff --git a/eis_toolkit/raster_processing/distance_to_anomaly.py b/eis_toolkit/raster_processing/distance_to_anomaly.py
index 598f95d2..c082d594 100644
--- a/eis_toolkit/raster_processing/distance_to_anomaly.py
+++ b/eis_toolkit/raster_processing/distance_to_anomaly.py
@@ -8,7 +8,7 @@
 from osgeo import gdal
 from rasterio import profiles
 
-from eis_toolkit.exceptions import EmptyDataException, InvalidParameterValueException
+from eis_toolkit.exceptions import EmptyDataException, InvalidParameterValueException, NumericValueSignException
 from eis_toolkit.utilities.checks.raster import check_raster_profile
 from eis_toolkit.utilities.miscellaneous import row_points, toggle_gdal_exceptions
 from eis_toolkit.vector_processing.distance_computation import distance_computation
@@ -84,26 +84,24 @@ def distance_to_anomaly(
 
 
 @beartype
-def distance_to_anomaly_gdal_compute_proximity(
+def distance_to_anomaly_gdal(
     anomaly_raster_profile: Union[profiles.Profile, dict],
     anomaly_raster_data: np.ndarray,
     threshold_criteria_value: Union[Tuple[Number, Number], Number],
     threshold_criteria: Literal["lower", "higher", "in_between", "outside"],
+    max_distance: Optional[Number] = None,
 ) -> Tuple[np.ndarray, profiles.Profile]:
     """Calculate distance from raster cell to nearest anomaly.
 
+    This tool is much faster than `distance_to_anomaly` but only available on
+    Windows.
+
     The criteria for what is anomalous can be defined as a single number and
     criteria text of "higher" or "lower". Alternatively, the definition can be
     a range where values inside (criteria text of "within") or outside are
     marked as anomalous (criteria text of "outside"). If anomaly_raster_profile does
     contain "nodata" key, np.nan is assumed to correspond to nodata values.
 
-    This function demonstrates superior performance compared to the distance_to_anomaly
-    and distance_to_anomaly_gdal functions, as it uses a low-level, C++-based API
-    within the GDAL library. By directly computing the proximity map from the
-    source dataset, it benefits from the core-level optimizations inherent to GDAL,
-    ensuring enhanced efficiency and speed.
-
     Args:
         anomaly_raster_profile: The raster profile in which the distances
             to the nearest anomalous value are determined.
@@ -114,6 +112,7 @@ def distance_to_anomaly_gdal_compute_proximity(
             the first value should be the minimum and the second
             the maximum value.
         threshold_criteria: Method to define anomalous.
+        max_distance: The maximum distance in the output array.
 
     Returns:
         A 2D numpy array with the distances to anomalies computed
@@ -124,12 +123,15 @@ def distance_to_anomaly_gdal_compute_proximity(
     _check_threshold_criteria_and_value(
         threshold_criteria=threshold_criteria, threshold_criteria_value=threshold_criteria_value
     )
+    if max_distance is not None and max_distance <= 0:
+        raise NumericValueSignException("Expected max distance to be a positive number.")
 
-    out_array, out_meta = _distance_to_anomaly_gdal_compute_proximity(
+    out_array, out_meta = _distance_to_anomaly_gdal(
         anomaly_raster_profile=anomaly_raster_profile,
         anomaly_raster_data=anomaly_raster_data,
         threshold_criteria=threshold_criteria,
         threshold_criteria_value=threshold_criteria_value,
+        max_distance=max_distance,
     )
 
     return out_array, out_meta
@@ -216,11 +218,12 @@ def _distance_to_anomaly(
     return distance_array
 
 
-def _distance_to_anomaly_gdal_compute_proximity(
+def _distance_to_anomaly_gdal(
     anomaly_raster_profile: Union[profiles.Profile, dict],
     anomaly_raster_data: np.ndarray,
     threshold_criteria_value: Union[Tuple[Number, Number], Number],
     threshold_criteria: Literal["lower", "higher", "in_between", "outside"],
+    max_distance: Optional[Number],
 ) -> Tuple[np.ndarray, profiles.Profile]:
 
     data_fits_criteria = _validate_threshold_criteria(
@@ -259,6 +262,8 @@ def _distance_to_anomaly_gdal_compute_proximity(
 
     # Create outputs
     out_array = out_band.ReadAsArray()
+    if max_distance is not None:
+        out_array[out_array > max_distance] = max_distance
     out_meta = anomaly_raster_profile.copy()
 
     # Update metadata
diff --git a/tests/raster_processing/test_distance_to_anomaly.py b/tests/raster_processing/test_distance_to_anomaly.py
index 0203a650..19585988 100644
--- a/tests/raster_processing/test_distance_to_anomaly.py
+++ b/tests/raster_processing/test_distance_to_anomaly.py
@@ -276,7 +276,7 @@ def test_distance_to_anomaly_nodata_handling(
     EXPECTED_PYTESTPARAMS,
 )
 @pytest.mark.xfail(sys.platform != "win32", reason="gdal_array available only on Windows.", raises=ModuleNotFoundError)
-def test_distance_to_anomaly_gdal_compute_proximity_expected(
+def test_distance_to_anomaly_gdal_expected(
     anomaly_raster_profile,
     anomaly_raster_data,
     threshold_criteria_value,
@@ -284,10 +284,10 @@ def test_distance_to_anomaly_gdal_compute_proximity_expected(
     expected_shape,
     expected_mean,
 ):
-    """Test distance_to_anomaly_gdal_compute_proximity_expected with expected result."""
+    """Test distance_to_anomaly_gdal with expected result."""
 
     assert not np.any(np.isnan(anomaly_raster_data))
-    out_image, out_profile = distance_to_anomaly.distance_to_anomaly_gdal_compute_proximity(
+    out_image, out_profile = distance_to_anomaly.distance_to_anomaly_gdal(
         anomaly_raster_profile=anomaly_raster_profile,
         anomaly_raster_data=anomaly_raster_data,
         threshold_criteria_value=threshold_criteria_value,
@@ -330,7 +330,7 @@ def test_distance_to_anomaly_gdal_compute_proximity_expected(
     ],
 )
 @pytest.mark.xfail(sys.platform != "win32", reason="gdal_array available only on Windows.", raises=ModuleNotFoundError)
-def test_distance_to_anomaly_gdal_compute_proximity_nodata_handling(
+def test_distance_to_anomaly_gdal_nodata_handling(
     anomaly_raster_profile,
     anomaly_raster_data,
     threshold_criteria_value,
@@ -339,12 +339,12 @@ def test_distance_to_anomaly_gdal_compute_proximity_nodata_handling(
     expected_mean_without_nodata,
     nodata_mask_value,
 ):
-    """Test distance_to_anomaly_gdal_compute_proximity with expected result."""
+    """Test distance_to_anomaly_gdal with expected result."""
 
     anomaly_raster_data_with_nodata = np.where(anomaly_raster_data > nodata_mask_value, np.nan, anomaly_raster_data)
     assert np.any(np.isnan(anomaly_raster_data_with_nodata))
 
-    out_image, out_profile = distance_to_anomaly.distance_to_anomaly_gdal_compute_proximity(
+    out_image, out_profile = distance_to_anomaly.distance_to_anomaly_gdal(
         anomaly_raster_profile=anomaly_raster_profile,
         anomaly_raster_data=anomaly_raster_data_with_nodata,
         threshold_criteria_value=threshold_criteria_value,
@@ -437,7 +437,7 @@ def test_distance_to_anomaly_gdal_compute_proximity_nodata_handling(
     ],
 )
 @pytest.mark.xfail(sys.platform != "win32", reason="gdal_array available only on Windows.", raises=ModuleNotFoundError)
-def test_distance_to_anomaly_gdal_compute_proximity_expected_check(
+def test_distance_to_anomaly_gdal_expected_check(
     anomaly_raster_profile,
     anomaly_raster_data,
     threshold_criteria_value,
@@ -445,12 +445,12 @@ def test_distance_to_anomaly_gdal_compute_proximity_expected_check(
     profile_additions,
     raises,
 ):
-    """Test distance_to_anomaly_gdal_compute_proximity checks."""
+    """Test distance_to_anomaly_gdal checks."""
 
     anomaly_raster_profile.update(profile_additions())
     anomaly_raster_profile_with_additions = anomaly_raster_profile
     with raises() as exc_info:
-        out_image, out_profile = distance_to_anomaly.distance_to_anomaly_gdal_compute_proximity(
+        out_image, out_profile = distance_to_anomaly.distance_to_anomaly_gdal(
             anomaly_raster_profile=anomaly_raster_profile_with_additions,
             anomaly_raster_data=anomaly_raster_data,
             threshold_criteria_value=threshold_criteria_value,