diff --git a/deltametrics/plan.py b/deltametrics/plan.py
index 64dc4df9..017908a3 100644
--- a/deltametrics/plan.py
+++ b/deltametrics/plan.py
@@ -717,6 +717,89 @@ def compute_shoreline_length(shore_mask, origin=[0, 0], return_line=False):
         return length
 
 
+def compute_shoreline_distance(shore_mask, origin=[0, 0],
+                               return_distances=False):
+    """Compute mean and stddev distance from the delta apex to the shoreline.
+
+    Algorithm computes the mean distance from the delta apex/origin to all
+    shoreline points.
+
+    .. important::
+
+        This calculation is subtly different than the "mean delta radius",
+        because the measurements are not sampled evenly along the opening
+        angle of the delta.
+
+    .. note:: uses `np.nanmean` and `np.nanstd`.
+
+    Parameters
+    ----------
+    shore_mask : :obj:`~deltametrics.mask.ShorelineMask`, :obj:`ndarray`
+        Shoreline mask. Can be a :obj:`~deltametrics.mask.ShorelineMask`
+        object, or a binarized array.
+
+    origin : :obj:`list`, :obj:`np.ndarray`, optional
+        Determines the location from where the distance to all shoreline
+        points is computed.
+
+    return_distances : :obj:`bool`
+        Whether to return the sorted line as a second argument. If True, a
+        ``Nx2`` array of x-y points is returned. Default is `False`.
+
+    Returns
+    -------
+    mean : :obj:`float`
+        Mean shoreline distance.
+
+    stddev : :obj:`float`
+        Standard deviation of shoreline distance.
+
+    distances : :obj:`np.ndarray`
+        If :obj:`return_distances` is `True`, then distance to each point
+        along the shoreline is *also* returned as an array (i.e., 3 arguments
+        are returned).
+
+    Examples
+    --------
+
+    .. doctest::
+
+        golf = dm.sample_data.golf()
+
+        sm = dm.mask.ShorelineMask(
+            golf['eta'][-1, :, :],
+            elevation_threshold=0,
+            elevation_offset=-0.5)
+
+        # compute mean and stddev distance
+        mean, stddev = dm.plan.compute_shoreline_distance(
+            sm, origin=[golf.meta['CTR'].data, golf.meta['L0'].data])
+
+    """
+    # check if mask or already array
+    if isinstance(shore_mask, mask.ShorelineMask):
+        _sm = shore_mask.mask
+    else:
+        _sm = shore_mask
+
+    if not (np.sum(_sm) > 0):
+        raise ValueError('No pixels in shoreline mask.')
+
+    if _sm.ndim == 3:
+        _sm = _sm.squeeze()
+
+    # find where the mask is True (all x-y pairs along shore)
+    _y, _x = np.argwhere(_sm).T
+
+    # determine the distances
+    _dists = np.sqrt((_x - origin[0])**2 + (_y - origin[1])**2)
+
+    if return_distances:
+        return np.nanmean(_dists), np.nanstd(_dists), _dists
+    else:
+        return np.nanmean(_dists), np.nanstd(_dists)
+
+
 @njit
 def _compute_angles_between(c1, shoreandborder, Shallowsea, numviews):
     """Private helper for shaw_opening_angle_method.
diff --git a/docs/source/reference/plan/index.rst b/docs/source/reference/plan/index.rst
index d970d18c..fff69c5c 100644
--- a/docs/source/reference/plan/index.rst
+++ b/docs/source/reference/plan/index.rst
@@ -27,4 +27,5 @@ Functions
 
 	compute_shoreline_roughness
 	compute_shoreline_length
+    compute_shoreline_distance
 	shaw_opening_angle_method
\ No newline at end of file
diff --git a/tests/test_plan.py b/tests/test_plan.py
index 11e3e0ee..579fb5ab 100644
--- a/tests/test_plan.py
+++ b/tests/test_plan.py
@@ -252,3 +252,53 @@ def test_rcm8_defaults_opposite(self):
             plt.show()
             breakpoint()
         assert len_1 == pytest.approx(self.rcm8_expected, abs=5.0)
+
+
+class TestShorelineDistance:
+
+    golf_path = _get_golf_path()
+    golf = cube.DataCube(golf_path)
+
+    sm = mask.ShorelineMask(
+        golf['eta'][-1, :, :],
+        elevation_threshold=0,
+        elevation_offset=-0.5)
+
+    def test_empty(self):
+        _arr = np.zeros((10, 10))
+        with pytest.raises(ValueError):
+            _, _ = plan.compute_shoreline_distance(_arr)
+
+    def test_single_point(self):
+        _arr = np.zeros((10, 10))
+        _arr[7, 5] = 1
+        mean00, stddev00 = plan.compute_shoreline_distance(
+            _arr)
+        mean05, stddev05 = plan.compute_shoreline_distance(
+            _arr, origin=[5, 0])
+        assert mean00 == np.sqrt(49 + 25)
+        assert mean05 == 7
+        assert stddev00 == 0
+        assert stddev05 == 0
+
+    def test_simple_case(self):
+        mean, stddev = plan.compute_shoreline_distance(
+            self.sm, origin=[self.golf.meta['CTR'].data,
+                             self.golf.meta['L0'].data])
+
+        assert mean > stddev
+        assert stddev > 0
+
+    def test_simple_case_distances(self):
+        m, s = plan.compute_shoreline_distance(
+            self.sm, origin=[self.golf.meta['CTR'].data,
+                             self.golf.meta['L0'].data])
+        m2, s2, dists = plan.compute_shoreline_distance(
+            self.sm, origin=[self.golf.meta['CTR'].data,
+                             self.golf.meta['L0'].data],
+            return_distances=True)
+
+        assert len(dists) > 0
+        assert np.mean(dists) == m
+        assert m2 == m
+        assert s2 == s