diff --git a/envs/conda-forge.yml b/envs/conda-forge.yml
index 1b6c024384..deced92d43 100644
--- a/envs/conda-forge.yml
+++ b/envs/conda-forge.yml
@@ -12,6 +12,7 @@ dependencies:
   # Included in improver-feedstock requirements
   - cartopy
   - clize
+  - geopandas
   - iris=3.12
   - numpy=2.2
   - scipy=1.15
diff --git a/envs/environment_a.yml b/envs/environment_a.yml
index 4ab3ad7f55..4fb9dcaf0d 100644
--- a/envs/environment_a.yml
+++ b/envs/environment_a.yml
@@ -16,6 +16,7 @@ dependencies:
   # Optional
   - fastparquet
   - statsmodels
+  - geopandas
   - numba
   - pygam
   - pysteps
diff --git a/envs/environment_b.yml b/envs/environment_b.yml
index c95b2db27a..4eac235c92 100644
--- a/envs/environment_b.yml
+++ b/envs/environment_b.yml
@@ -18,6 +18,7 @@ dependencies:
   - scipy=1.15
   - sphinx
   # Optional
+  - geopandas
   - lightgbm>=4.0.0
   - numba
   - python-stratify
diff --git a/envs/latest.yml b/envs/latest.yml
index 1fe7ed5332..f179b01513 100644
--- a/envs/latest.yml
+++ b/envs/latest.yml
@@ -18,6 +18,7 @@ dependencies:
   - sphinx
   # Optional
   - python-stratify
+  - geopandas
   - statsmodels
   - lightgbm
   - numba
diff --git a/improver/generate_ancillaries/generate_distance_to_feature.py b/improver/generate_ancillaries/generate_distance_to_feature.py
new file mode 100644
index 0000000000..c9e8c62f18
--- /dev/null
+++ b/improver/generate_ancillaries/generate_distance_to_feature.py
@@ -0,0 +1,285 @@
+# (C) Crown Copyright, Met Office. All rights reserved.
+#
+# This file is part of 'IMPROVER' and is released under the BSD 3-Clause license.
+# See LICENSE in the root of the repository for full licensing details.
+"""A module for generating a distance to feature ancillary cube."""
+
+from typing import List, Optional, Tuple
+
+import pyproj
+from geopandas import GeoDataFrame, GeoSeries, clip
+from iris.cube import Cube
+from numpy import array, min, round
+from shapely.geometry import Point
+
+from improver import BasePlugin
+
+
+class DistanceTo(BasePlugin):
+    """
+    Plugin to calculate the distance to the nearest feature in a geometry from
+    sites to the nearest metre.
+
+    Given a cube containing site locations and a GeoDataFrame, the distance to each site
+    from the nearest point of the feature geometry is calculated. This is done by converting
+    the feature geometry and sites to a common target orography that must be specified using a
+    European Petroleum Survey Group (EPSG) code that identifies the projection. For the
+    UK, EPSG code 3035 may be used to provide a Lambert Azimuthal Equal Area projection that is suitable for the region. The chosen projection should match the projection on which the
+    ancillary will be used. The distance method from Shapely is used to find the distance
+    from each site to every point in the feature geometry. The minimum of these distances is
+    returned as the distance to the nearest feature in the feature geometry and this is rounded
+    to the nearest metre.
+
+    If requested, the provided geometry can be clipped to the site bounds with a buffer to reduce computation.
+    locations with a buffer to improve performance. This is useful when the geometry is
+    large and it would be expensive to calculate the distance to all features in the
+    geometry but information may be lost at the edges of the domain.
+    """
+
+    def __init__(
+        self,
+        epsg_projection: int,
+        new_name: Optional[str] = None,
+        buffer: float = 30000,
+        clip_geometry_flag: bool = False,
+    ) -> None:
+        """
+        Initialise the DistanceTo plugin.
+
+        Args:
+            epsg_projection:
+                The EPSG code of the coordinate reference system on to which latitudes
+                and longitudes will be projected to calculate distances. This is
+                a projected coordinate system in which distances are measured in metres,
+                for example, EPSG code 3035, which defines a Lambert Azimuthal Equal
+                Areas projection suitable for the UK.
+            new_name:
+                The name of the output cube.
+            buffer:
+                A buffer distance in m. If the geometry is clipped, this distance will
+                be added onto the outermost site locations to define the domain to clip
+                the geometry to.
+            clip_geometry_flag:
+                A flag to indicate whether the geometry should be clipped to the bounds of
+                the site locations with a buffer distance added to the bounds. If set to
+                False, the full geometry will be used to calculate the distance to the
+                nearest feature.
+        """
+        self.epsg_projection = epsg_projection
+        self.new_name = new_name
+        self.buffer = buffer
+        self.clip_geometry_flag = clip_geometry_flag
+
+    @staticmethod
+    def get_clip_values(points: List[float], buffer: float) -> List[float]:
+        """Get the coordinates to use when clipping the geometry. This is determined by
+        finding the maximum and minimum coordinate points from a list. A buffer distance
+        may then be added/subtracted to the max/min.
+
+        Args:
+            points:
+                A list of points to find the min and max from.
+            buffer:
+                The buffer distance to add/subtract to the max/min points.
+        Returns:
+            A list containing the maximum and minimum points with the buffer added
+            or subtracted respectively."""
+
+        ordered_points = sorted(set(points))
+        min_point = ordered_points[0] - buffer
+        max_point = ordered_points[-1] + buffer
+
+        return [min_point, max_point]
+
+    def clip_geometry(
+        self, geometry: GeoDataFrame, bounds_x: List[float], bounds_y: List[float]
+    ) -> GeoDataFrame:
+        """Clip the geometry to the provided bounds.
+
+        Args:
+            geometry:
+                The geometry to clip.
+            bounds_x:
+                A list containing the minimum and maximum x coordinates.
+            bounds_y:
+                A list containing the minimum and maximum y coordinates.
+        Returns:
+            The clipped geometry.
+
+        Raises:
+            ValueError: If the clipped geometry is empty after clipping with the
+            provided bounds."""
+
+        clipped_geometry = clip(
+            geometry, mask=[bounds_x[0], bounds_y[0], bounds_x[1], bounds_y[1]]
+        )
+        if clipped_geometry.empty:
+            raise ValueError(
+                "Clipping the geometry with a buffer size of "
+                f"{self.buffer}m has produced an empty geometry. Either "
+                "increase the buffer size or set clip_geometry_flag to "
+                "False to use the full geometry."
+            )
+
+        return clipped_geometry
+
+    def check_target_crs(self, site_points: GeoSeries):
+        """Check that the provided target projection is suitable for the sites
+        being used.
+
+        Args:
+            site_points:
+                A GeoSeries containing the site points.
+
+        Raises:
+            ValueError: If the provided target coordinate reference system is not
+                        suitable for the site points.
+        """
+        x_bounds = self.get_clip_values(site_points.x, 0)
+        y_bounds = self.get_clip_values(site_points.y, 0)
+
+        target_crs = pyproj.CRS.from_epsg(self.epsg_projection)
+        x_min, y_min, x_max, y_max = target_crs.area_of_use.bounds
+
+        valid_target = (
+            x_bounds[0] > x_min
+            and x_bounds[1] < x_max
+            and y_bounds[0] > y_min
+            and y_bounds[1] < y_max
+        )
+        if not valid_target:
+            raise ValueError(
+                "The provided projection defined by EPSG code "
+                f"{self.epsg_projection} is not suitable for the site "
+                "locations provided. Limits of this domain are: "
+                f"x: {x_min} to {x_max}, y: {y_min} to {y_max}, whilst "
+                f"the site locations are bounded by x: {x_bounds[0]} to {x_bounds[1]}, "
+                f"y: {y_bounds[0]} to {y_bounds[1]}."
+            )
+
+    def project_geometry(
+        self, geometry: GeoDataFrame, site_cube: Cube
+    ) -> Tuple[GeoSeries, GeoDataFrame]:
+        """Project the geometry and site cube to the target projection.
+
+        Args:
+            geometry:
+                The geometry to reproject.
+            site_cube:
+                The cube containing the site locations. It is assumed that the site
+                coordinates are defined as latitude and longitude on a WGS84
+                coordinate system (EPSG:4326).
+        Returns:
+            A tuple containing the projected site locations and geometry."""
+
+        x_points = site_cube.coord(axis="x").points
+        y_points = site_cube.coord(axis="y").points
+
+        site_points_list = [Point(x, y) for x, y in zip(x_points, y_points)]
+        # Assumes site coordinates are defined as latitude and longitude
+        # defaulting to an EPSG:4326 coordinate system, which is WGS84.
+        site_points = GeoSeries(site_points_list, crs=4326)
+
+        # Check that the provided target projection is suitable for the site points
+        self.check_target_crs(site_points)
+
+        geometry_reprojection = geometry.to_crs(self.epsg_projection)
+        site_points = site_points.to_crs(self.epsg_projection)
+        return site_points, geometry_reprojection
+
+    def distance_to(self, site_points: GeoSeries, geometry: GeoDataFrame) -> List[int]:
+        """Calculate the distance from each site point to the nearest feature in the
+        geometry.
+
+        Args:
+            site_points:
+                A GeoSeries containing the site points in the target projection.
+            geometry:
+                A GeoDataFrame containing the geometry in the target projection.
+        Returns:
+            A list of distances from each site point to the nearest feature in the
+            geometry rounded to the nearest metre."""
+        distance_results = []
+        for point in site_points:
+            distance_to_nearest = min(point.distance(geometry.geometry))
+            distance_results.append(round(distance_to_nearest))
+
+        return distance_results
+
+    def create_output_cube(self, site_cube: Cube, data: List[int]) -> Cube:
+        """Create an output cube that will have the same metadata as the input site
+        cube except the units are changed to metres and, if requested, the name of the
+        output cube will be changed.
+
+        Args:
+           site_cube:
+               The input cube containing site locations that are defined by latitude
+               and longitude coordinates.
+           data:
+               A list of distances from each site point to the nearest feature in the
+               geometry.
+           Returns:
+               A new cube containing the distances with the same metadata as input site
+               cube but with updated units and name."""
+
+        output_cube = site_cube.copy(data=array(data))
+        if self.new_name:
+            output_cube.rename(self.new_name)
+        output_cube.units = "m"
+
+        return output_cube
+
+    def process(self, site_cube: Cube, geometry: GeoDataFrame) -> Cube:
+        """Generate a cube of the distance from sites in site_cube to the
+        nearest point in geometry.
+
+        The latitude, longitude coordinates in the site_cube are extracted
+        to define the location of the sites. The sites and feature geometry
+        are reprojected to the target projection.
+
+        If requested the feature geometry will be clipped to the smallest square possible
+        such that all sites in site_cube are included. A buffer distance is then added to each
+        edge of the square which defines the size the feature geometry will be clipped to. This
+        is useful when the feature geometry size is large and it would be expensive to calculate
+        the distance to all features in the geometry or where the domain of the feature geometry
+        is much larger than the area containing the site locations. Information may be lost at the edges of
+        the domain if the feature is sparsely located in the geometry.
+
+        The distance from each site to every point in the feature geometry is then calculated
+        and the minimum of these distances is returned. The distances are rounded to the
+        nearest metre.
+
+        The output cube will have the same metadata as the input site_cube except the
+        units will be changed to meters and, if requested, the name of the output cube
+        will be updated.
+
+        Args:
+            site_cube:
+                The input cube containing site locations. This cube must have x and y
+                axis which contain the site coordinates in latitude and longitude.
+            geometry:
+                The GeoDataFrame containing the geometry to calculate distances to.
+
+        Returns:
+            A new cube containing the distances from each site to the nearest feature
+            in the geometry rounded to the nearest metre."""
+
+        # Project the geometry and site cube coordinates to the target projection.
+        site_coords, geometry_projection = self.project_geometry(geometry, site_cube)
+
+        if self.clip_geometry_flag:
+            # Clip the geometry to the bounds of the site coordinates with a buffer if
+            # requested
+            x_bounds = self.get_clip_values(site_coords.x, self.buffer)
+            y_bounds = self.get_clip_values(site_coords.y, self.buffer)
+
+            clipped_geometry = self.clip_geometry(
+                geometry_projection, x_bounds, y_bounds
+            )
+        else:
+            clipped_geometry = geometry_projection
+
+        # Calculate the distance to the nearest feature in the geometry
+        distance_to_results = self.distance_to(site_coords, clipped_geometry)
+
+        return self.create_output_cube(site_cube, distance_to_results)
diff --git a/improver_tests/generate_ancillaries/test_DistanceTo.py b/improver_tests/generate_ancillaries/test_DistanceTo.py
new file mode 100644
index 0000000000..51d79bde24
--- /dev/null
+++ b/improver_tests/generate_ancillaries/test_DistanceTo.py
@@ -0,0 +1,504 @@
+# (C) Crown Copyright, Met Office. All rights reserved.
+#
+# This file is part of 'IMPROVER' and is released under the BSD 3-Clause license.
+# See LICENSE in the root of the repository for full licensing details.
+"""Unit tests for the DistanceTo plugin."""
+
+import numpy as np
+import pytest
+from geopandas import GeoDataFrame
+from shapely.geometry import LineString, Point, Polygon
+
+from improver.generate_ancillaries.generate_distance_to_feature import DistanceTo
+from improver.spotdata.build_spotdata_cube import build_spotdata_cube
+
+
+@pytest.fixture()
+def geometry_point_latlon():
+    """Create a geometry containing points on a latitude, longitude grid.
+    The location of the points is identical to geometry_point_laea, but in a different
+    CRS.
+
+    The points locations look like:
+             x        x
+
+
+             x        x
+    """
+    data = [
+        Point(-1.393298, 49.538352),
+        Point(-1.395401, 49.547221),
+        Point(-1.379621, 49.539742),
+        Point(-1.381721, 49.548611),
+    ]
+
+    return GeoDataFrame(geometry=data, crs="EPSG:4326")
+
+
+@pytest.fixture()
+def geometry_point_laea():
+    """Create a geometry containing points on a Lambert azimuthal equal-area grid.
+    The location of the points is identical to geometry_point_latlon, but in a
+    different CRS.
+
+    The points locations look like:
+             x        x
+
+
+             x        x
+    """
+    data = [
+        Point(3500000, 3000000),
+        Point(3500000, 3001000),
+        Point(3501000, 3000000),
+        Point(3501000, 3001000),
+    ]
+    return GeoDataFrame(geometry=data, crs="EPSG:3035")
+
+
+@pytest.fixture()
+def geometry_line_latlon():
+    """Create a simple line geometry on a latitude, longitude grid.
+    The line defined is identical to geometry_line_laea, but in a different CRS.
+
+    The line looks like:
+             x-------x
+             |       |
+             |       |
+             |       |
+             x-------x
+    """
+    data = [
+        LineString(
+            [
+                [-1.393298, 49.538352],
+                [-1.395401, 49.547221],
+                [-1.381721, 49.548611],
+                [-1.379621, 49.539742],
+                [-1.393298, 49.538352],
+            ]
+        )
+    ]
+
+    return GeoDataFrame(geometry=data, crs="EPSG:4326")
+
+
+@pytest.fixture()
+def geometry_line_laea():
+    """Create a simple line geometry on a Lambert azimuthal equal area grid.
+    The line defined is identical to geometry_point_latlon, but in a different CRS.
+
+    The line looks like:
+             x-------x
+             |       |
+             |       |
+             |       |
+             x-------x
+    """
+    data = [
+        LineString(
+            [
+                [3500000, 3000000],
+                [3500000, 3001000],
+                [3501000, 3001000],
+                [3501000, 3000000],
+                [3500000, 3000000],
+            ]
+        )
+    ]
+    return GeoDataFrame(geometry=data, crs="EPSG:3035")
+
+
+@pytest.fixture()
+def geometry_polygon_latlon():
+    """Create a simple polygon geometry on a latitude, longitude grid.
+    The polygon defined is identical to geometry_polygon_laea, but in a different CRS.
+
+    The polygon looks like:
+             x-------x
+             ---------
+             ---------
+             x-------x
+    """
+    data = [
+        Polygon(
+            [
+                [-1.393298, 49.538352],
+                [-1.395401, 49.547221],
+                [-1.381721, 49.548611],
+                [-1.379621, 49.539742],
+                [-1.393298, 49.538352],
+            ]
+        )
+    ]
+    return GeoDataFrame(geometry=data, crs="EPSG:4326")
+
+
+@pytest.fixture()
+def geometry_polygon_laea():
+    """Create a simple polygon geometry on a Lambert azimuthal equal area grid.
+    The polygon defined is identical to geometry_polygon_latlon, but in a different CRS.
+
+    The polygon looks like:
+             x-------x
+             ---------
+             ---------
+             x-------x
+    """
+    data = [
+        Polygon(
+            [
+                [3500000, 3000000],
+                [3500000, 3001000],
+                [3501000, 3001000],
+                [3501000, 3000000],
+                [3500000, 3000000],
+            ]
+        )
+    ]
+    return GeoDataFrame(geometry=data, crs="EPSG:3035")
+
+
+@pytest.fixture()
+def single_site_cube():
+    """Set up a site cube for a single site."""
+
+    latitude = 49.539047274  # This value is overridden in the test functions.
+    longitude = -1.386459578  # This value is overridden in the test functions.
+
+    altitude = -99999  # This value is not used but is required for cube creation.
+    data = -99999  # This value is not used but is required for cube creation.
+    wmo_id = ["00000"]  # This value is not used but is required for cube creation.
+
+    prob_cube = build_spotdata_cube(
+        data,
+        name="rain_rate",
+        units="1",
+        altitude=altitude,
+        wmo_id=wmo_id,
+        latitude=latitude,
+        longitude=longitude,
+    )
+    return prob_cube
+
+
+@pytest.fixture()
+def multiple_site_cube():
+    """Set up a site cube containing data at multiple sites."""
+
+    latitude = np.array([49.538352, 49.539047274, 49.543481633, 49.552350289])
+    longitude = np.array([-1.393298, -1.386459578, -1.387510304, -1.389612479])
+
+    altitude = np.array(
+        [-99999, -99999, -99999, -99999]
+    )  # These values are not used but are required for cube creation.
+    data = np.array(
+        [-99999, -99999, -99999, -99999]
+    )  # These values are not used but are required for cube creation.
+    wmo_id = [
+        "00000",
+        "00001",
+        "00002",
+        "00003",
+    ]  # These values are not used but are required for cube creation.
+
+    prob_cube = build_spotdata_cube(
+        data,
+        name="rain_rate",
+        units="1",
+        altitude=altitude,
+        wmo_id=wmo_id,
+        latitude=latitude,
+        longitude=longitude,
+    )
+    return prob_cube
+
+
+@pytest.mark.parametrize(
+    "target_projection, site_latitude, site_longitude, expected_distance",
+    [
+        (3035, 49.538352, -1.393298, 0),  # site is the same location as a point
+        (3035, 49.539047274, -1.386459578, 500),  # site is halfways between two points
+        (3035, 49.543481633, -1.387510304, 707),  # Site at centre of the 4 points
+        # Test a conic projection over Europe as well, which yields difference distances
+        # for the non-zero distance cases.
+        (9001, 49.538352, -1.393298, 0),  # site is the same location as a point
+        (9001, 49.539047274, -1.386459578, 498),  # site is halfways between two points
+        (9001, 49.543481633, -1.387510304, 603),  # Site at centre of the 4 points
+    ],
+)
+@pytest.mark.parametrize(
+    "shape_file_crs", ["geometry_point_laea", "geometry_point_latlon"]
+)
+def test_distance_to_with_points_geometry(
+    single_site_cube,
+    shape_file_crs,
+    target_projection,
+    site_latitude,
+    site_longitude,
+    expected_distance,
+    request,
+):
+    """Test the DistanceTo plugin with a single site and a geometry
+    of points."""
+
+    geometry = request.getfixturevalue(shape_file_crs)
+
+    single_site_cube.coord("latitude").points = site_latitude
+    single_site_cube.coord("longitude").points = site_longitude
+
+    output_cube = DistanceTo(target_projection)(single_site_cube, geometry)
+    assert output_cube.name() == "rain_rate"
+    assert output_cube.units == "m"
+    assert output_cube.coord("latitude").points == site_latitude
+    assert output_cube.coord("longitude").points == site_longitude
+    assert output_cube.data == expected_distance
+
+
+@pytest.mark.parametrize(
+    "target_projection, site_latitude, site_longitude, expected_distance",
+    [
+        (
+            3035,
+            49.538352,
+            -1.393298,
+            0,
+        ),  # site is the same location as a corner of the line
+        (
+            3035,
+            49.539047274,
+            -1.386459578,
+            0,
+        ),  # site is halfways between two points on the line
+        (
+            3035,
+            49.543481633,
+            -1.387510304,
+            500,
+        ),  # Site is at the exact centre of the square formed by the line
+        # Test a conic projection over Europe as well, which yields difference distances
+        # for the non-zero distance cases.
+        (
+            9001,
+            49.538352,
+            -1.393298,
+            0,
+        ),  # site is the same location as a corner of the line
+        (
+            9001,
+            49.539047274,
+            -1.386459578,
+            0,
+        ),  # site is halfways between two points on the line
+        (
+            9001,
+            49.543481633,
+            -1.387510304,
+            382,
+        ),  # Site is at the exact centre of the square formed by the line
+    ],
+)
+@pytest.mark.parametrize("geometry_crs", ["geometry_line_laea", "geometry_line_latlon"])
+def test_distance_to_with_line_geometry(
+    single_site_cube,
+    geometry_crs,
+    target_projection,
+    site_latitude,
+    site_longitude,
+    expected_distance,
+    request,
+):
+    """Test the DistanceTo plugin with a single site and a
+    single line geometry."""
+
+    geometry = request.getfixturevalue(geometry_crs)
+
+    single_site_cube.coord("latitude").points = site_latitude
+    single_site_cube.coord("longitude").points = site_longitude
+
+    output_cube = DistanceTo(target_projection)(single_site_cube, geometry)
+    assert output_cube.name() == "rain_rate"
+    assert output_cube.units == "m"
+    assert output_cube.coord("latitude").points == site_latitude
+    assert output_cube.coord("longitude").points == site_longitude
+    assert output_cube.data == expected_distance
+
+
+@pytest.mark.parametrize(
+    "target_projection, site_latitude, site_longitude, expected_distance",
+    [
+        (
+            3035,
+            49.538352,
+            -1.393298,
+            0,
+        ),  # site is the same location as a corner of the polygon
+        (
+            3035,
+            49.539047274,
+            -1.386459578,
+            0,
+        ),  # site is halfways between two points on the edge of the polygon
+        (
+            3035,
+            49.543481633,
+            -1.387510304,
+            0,
+        ),  # Site is at the exact centre of the polygon
+        (3035, 49.551655272, -1.3964531, 500),  # Site is outside the polygon
+        # Test a conic projection over Europe as well, which yields difference distances
+        # for the non-zero distance cases.
+        (
+            9001,
+            49.538352,
+            -1.393298,
+            0,
+        ),  # site is the same location as a corner of the polygon
+        (
+            9001,
+            49.539047274,
+            -1.386459578,
+            0,
+        ),  # site is halfways between two points on the edge of the polygon
+        (
+            9001,
+            49.543481633,
+            -1.387510304,
+            0,
+        ),  # Site is at the exact centre of the polygon
+        (9001, 49.551655272, -1.3964531, 383),  # Site is outside the polygon
+    ],
+)
+@pytest.mark.parametrize(
+    "geometry_crs", ["geometry_polygon_laea", "geometry_polygon_latlon"]
+)
+def test_distance_to_with_polygon_geometry(
+    single_site_cube,
+    geometry_crs,
+    target_projection,
+    site_latitude,
+    site_longitude,
+    expected_distance,
+    request,
+):
+    """Test the DistanceTo plugin with a single site and a simple polygon geometry."""
+
+    geometry = request.getfixturevalue(geometry_crs)
+
+    single_site_cube.coord("latitude").points = site_latitude
+    single_site_cube.coord("longitude").points = site_longitude
+
+    output_cube = DistanceTo(target_projection)(single_site_cube, geometry)
+    assert output_cube.name() == "rain_rate"
+    assert output_cube.units == "m"
+    assert output_cube.coord("latitude").points == site_latitude
+    assert output_cube.coord("longitude").points == site_longitude
+    assert output_cube.data == expected_distance
+
+
+@pytest.mark.parametrize(
+    "geometry_type,expected_distance",
+    [
+        ("point", [0, 500, 707, 707]),
+        ("line", [0, 0, 500, 500]),
+        ("polygon", [0, 0, 0, 500]),
+    ],
+)
+@pytest.mark.parametrize("geometry_crs", ("laea", "latlon"))
+def test_distance_to_with_multiple_sites(
+    multiple_site_cube,
+    geometry_type,
+    geometry_crs,
+    expected_distance,
+    request,
+):
+    """Test the DistanceTo plugin works when provided a site cube with multiple sites and
+    different types of geometry"""
+
+    geometry = request.getfixturevalue(f"geometry_{geometry_type}_{geometry_crs}")
+
+    output_cube = DistanceTo(3035)(multiple_site_cube, geometry)
+    assert output_cube.name() == "rain_rate"
+    assert output_cube.units == "m"
+
+    np.testing.assert_allclose(output_cube.data, expected_distance)
+
+
+def test_distance_to_with_new_name(single_site_cube, geometry_point_laea):
+    """Test the DistanceTo plugin correctly sets a new name."""
+
+    single_site_cube.coord("latitude").points = 49.539047274
+    single_site_cube.coord("longitude").points = -1.386459578
+
+    output_cube = DistanceTo(3035, new_name="distance_to_river")(
+        single_site_cube, geometry_point_laea
+    )
+    assert output_cube.name() == "distance_to_river"
+    assert output_cube.units == "m"
+    assert output_cube.coord("latitude").points == 49.539047274
+    assert output_cube.coord("longitude").points == -1.386459578
+
+
+@pytest.mark.parametrize(
+    "clip, buffer, expected",
+    [(True, 100, [100, 800]), (True, 3000, [100, 200]), (False, None, [100, 200])],
+)
+def test_distance_to_clipping_loss_of_data(
+    multiple_site_cube,
+    clip,
+    buffer,
+    expected,
+):
+    """Test the DistanceTo plugin with clipping and buffer. The test involves two sites
+    (o) and two features (x) configured as follows:
+
+        (-100) (0)     (800) (1000)
+            o   x        o    x
+
+    The numbers represent their relative distance to each other in metres.
+    """
+    site_cubes = multiple_site_cube[0:2].copy()  # Use only the first two sites
+    site_cubes.coord("latitude").points = [49.537465617, 49.545447324]
+    site_cubes.coord("longitude").points = [-1.393088146, -1.394980629]
+
+    data = [
+        Point(3500000, 3000000),
+        Point(3500000, 3001000),
+    ]
+    geometry = GeoDataFrame(geometry=data, crs="EPSG:3035")
+
+    if clip:
+        output_cube = DistanceTo(3035, clip_geometry_flag=True, buffer=buffer)(
+            site_cubes, geometry
+        )
+    else:
+        output_cube = DistanceTo(3035, clip_geometry_flag=False)(site_cubes, geometry)
+
+    assert output_cube.name() == "rain_rate"
+    assert output_cube.units == "m"
+    np.testing.assert_allclose(output_cube.data, expected)
+
+
+def test_distance_to_with_empty_geometry(single_site_cube, geometry_point_laea):
+    """Test the DistanceTo plugin raises a ValueError when clipping leads to an empty
+    geometry."""
+
+    with pytest.raises(
+        ValueError, match="Clipping the geometry with a buffer size of 100m"
+    ):
+        DistanceTo(3035, clip_geometry_flag=True, buffer=100)(
+            single_site_cube, geometry_point_laea
+        )
+
+
+def test_distance_to_with_unsuitable_projection(single_site_cube, geometry_point_laea):
+    """Test the DistanceTo plugin raises a ValueError when the projection is unsuitable."""
+
+    msg = (
+        "The provided projection defined by EPSG code 3112 is not suitable "
+        "for the site locations provided. Limits of this domain are: x: 112.85 "
+        "to 153.69, y: -43.7 to -9.86, whilst the site locations are bounded by "
+        "x: -1.386459578 to -1.386459578, y: 49.539047274 to 49.539047274."
+    )
+    with pytest.raises(ValueError, match=msg):
+        DistanceTo(3112)(single_site_cube, geometry_point_laea)