diff --git a/pyresample/future/geometry/_subset.py b/pyresample/future/geometry/_subset.py
index bb04cf69..c6136067 100644
--- a/pyresample/future/geometry/_subset.py
+++ b/pyresample/future/geometry/_subset.py
@@ -5,35 +5,61 @@
 from typing import TYPE_CHECKING, Any
 
 import numpy as np
+from pyproj import Proj
 
 # this caching module imports the geometries so this subset module
 # must be imported inside functions in the geometry modules if needed
 # to avoid circular dependencies
 from pyresample._caching import cache_to_json_if
 from pyresample.boundary import Boundary
-from pyresample.geometry import get_geostationary_bounding_box_in_lonlats, logger
+from pyresample.geometry import (
+    DEFAULT_AREA_SLICE_SAMPLE_STEPS,
+    get_geostationary_bounding_box_in_lonlats,
+    logger,
+)
 from pyresample.utils import check_slice_orientation
 
 if TYPE_CHECKING:
     from pyresample import AreaDefinition
 
 
+MAX_POINTS_PER_CHUNK = 600_000
+
+
 @cache_to_json_if("cache_geometry_slices")
 def get_area_slices(
         src_area: AreaDefinition,
         area_to_cover: AreaDefinition,
         shape_divisible_by: int | None,
+        sample_steps: int | None = DEFAULT_AREA_SLICE_SAMPLE_STEPS,
+        sample_grid: bool = False,
 ) -> tuple[slice, slice]:
-    """Compute the slice to read based on an `area_to_cover`."""
+    """Compute the slice to read based on an `area_to_cover`.
+
+    For geostationary source areas in cross-projection mode:
+    - ``sample_steps`` >= 2 and ``sample_grid=False`` samples edge points only.
+    - ``sample_steps`` >= 2 and ``sample_grid=True`` samples an interior grid.
+    - ``sample_steps`` <= 0 or ``None`` and ``sample_grid=False`` samples all edge points.
+    - ``sample_steps`` <= 0 or ``None`` and ``sample_grid=True`` samples all destination points.
+    - ``sample_steps`` == 1 raises ``ValueError``.
+    """
     if not _is_area_like(src_area):
         raise NotImplementedError(f"Only AreaDefinitions are supported, not {type(src_area)}")
     if not _is_area_like(area_to_cover):
         raise NotImplementedError(f"Only AreaDefinitions are supported, not {type(area_to_cover)}")
 
-    # Intersection only required for two different projections
-    proj_def_to_cover = area_to_cover.crs
-    proj_def = src_area.crs
-    if proj_def_to_cover == proj_def:
+    normalized_sample_steps = _normalize_sample_steps(sample_steps)
+
+    # Intersection is only required for two different projections.
+    src_crs_wkt = getattr(src_area, "crs_wkt", None)
+    dst_crs_wkt = getattr(area_to_cover, "crs_wkt", None)
+    if src_crs_wkt is not None and src_crs_wkt == dst_crs_wkt:
+        proj_def_to_cover = src_crs_wkt
+    else:
+        proj_def_to_cover = area_to_cover.crs
+        if proj_def_to_cover != src_area.crs:
+            proj_def_to_cover = None
+    if proj_def_to_cover is not None:
         logger.debug('Projections for data and slice areas are identical: %s',
                      proj_def_to_cover)
         # Get slice parameters
@@ -45,6 +71,16 @@ def get_area_slices(
         y_slice = _ensure_integer_slice(y_slice)
         return x_slice, y_slice
 
+    if src_area.is_geostationary:
+        coverage_slices = _get_covered_source_slices(
+            src_area,
+            area_to_cover,
+            sample_steps=normalized_sample_steps,
+            sample_grid=sample_grid,
+        )
+        if coverage_slices is not None:
+            return _finalize_slices(src_area, coverage_slices[0], coverage_slices[1], shape_divisible_by)
+
     data_boundary = _get_area_boundary(src_area)
     area_boundary = _get_area_boundary(area_to_cover)
     intersection = data_boundary.contour_poly.intersection(
@@ -57,6 +93,33 @@ def get_area_slices(
         np.rad2deg(intersection.lon), np.rad2deg(intersection.lat))
     x_slice = slice(np.ma.min(x), np.ma.max(x) + 1)
     y_slice = slice(np.ma.min(y), np.ma.max(y) + 1)
+
+    return _finalize_slices(src_area, x_slice, y_slice, shape_divisible_by)
+
+
+def _normalize_sample_steps(sample_steps: int | None):
+    """Normalize sampling config to sampled mode or ALL-mode fallback.
+
+    ``None`` and values ``<= 0`` map to ``ALL`` mode, whose exact behavior is
+    determined later by ``sample_grid``.
+    A value of ``1`` is rejected because it does not provide meaningful sampled
+    coverage for edge or grid modes.
+    """
+    if sample_steps is None:
+        return None
+    try:
+        normalized_sample_steps = int(sample_steps)
+    except (TypeError, ValueError) as err:
+        raise ValueError(f"sample_steps must be an integer or None, got {sample_steps!r}") from err
+    if normalized_sample_steps <= 0:
+        return None
+    if normalized_sample_steps == 1:
+        raise ValueError("sample_steps=1 is not supported; use <= 0/None for ALL mode or >= 2 for sampled modes.")
+    return normalized_sample_steps
+
+
+def _finalize_slices(src_area: AreaDefinition, x_slice: slice, y_slice: slice, shape_divisible_by: int | None):
+    """Normalize slice bounds and apply orientation/divisibility rules."""
     x_slice = _ensure_integer_slice(x_slice)
     y_slice = _ensure_integer_slice(y_slice)
     if shape_divisible_by is not None:
@@ -105,6 +168,137 @@ def _get_area_boundary(area_to_cover: AreaDefinition) -> Boundary:
         raise NotImplementedError("Can't determine boundary of area to cover") from err
 
 
+def _get_covered_source_slices(
+        src_area: AreaDefinition,
+        area_to_cover: AreaDefinition,
+        sample_steps: int | None,
+        sample_grid: bool,
+):
+    """Estimate covering source slices from sampled destination points.
+
+    Returns ``None`` when sampled points do not produce any valid source
+    coverage, allowing the caller to fall back to boundary intersection.
+    """
+    min_col = None
+    max_col = None
+    min_row = None
+    max_row = None
+    try:
+        src_proj = Proj(src_area.crs)
+        for destination_lons, destination_lats in _iter_destination_lonlat_samples(
+            area_to_cover=area_to_cover,
+            sample_steps=sample_steps,
+            sample_grid=sample_grid,
+        ):
+            source_xs, source_ys = src_proj(
+                destination_lons,
+                destination_lats,
+            )
+            source_cols, source_rows = src_area.get_array_indices_from_projection_coordinates(
+                source_xs,
+                source_ys,
+            )
+            valid = ~np.ma.getmaskarray(source_cols) & ~np.ma.getmaskarray(source_rows)
+            if not valid.any():
+                continue
+            chunk_cols = np.ma.getdata(source_cols)[valid]
+            chunk_rows = np.ma.getdata(source_rows)[valid]
+            chunk_min_col = int(chunk_cols.min())
+            chunk_max_col = int(chunk_cols.max())
+            chunk_min_row = int(chunk_rows.min())
+            chunk_max_row = int(chunk_rows.max())
+            min_col = chunk_min_col if min_col is None else min(min_col, chunk_min_col)
+            max_col = chunk_max_col if max_col is None else max(max_col, chunk_max_col)
+            min_row = chunk_min_row if min_row is None else min(min_row, chunk_min_row)
+            max_row = chunk_max_row if max_row is None else max(max_row, chunk_max_row)
+    except (RuntimeError, TypeError, ValueError):
+        logger.debug("Failed to estimate covered source slices from sampled destination points.", exc_info=True)
+        return None
+    if min_col is None or min_row is None or max_col is None or max_row is None:
+        return None
+    col_start = max(0, min_col)
+    col_stop = min(src_area.width, max_col + 1)
+    row_start = max(0, min_row)
+    row_stop = min(src_area.height, max_row + 1)
+    if col_start >= col_stop or row_start >= row_stop:
+        return None
+    return slice(col_start, col_stop), slice(row_start, row_stop)
+
+
+def _iter_destination_lonlat_samples(
+        area_to_cover: AreaDefinition,
+        sample_steps: int | None,
+        sample_grid: bool,
+):
+    """Yield destination lon/lat samples for ALL, grid, or edge mode."""
+    if sample_steps is None:
+        if sample_grid:
+            yield from _iter_dense_lonlat_samples(area_to_cover)
+            return
+        yield _get_edge_lonlat_samples(area_to_cover, sample_steps=None)
+        return
+    if sample_grid:
+        yield _get_grid_lonlat_samples(area_to_cover, sample_steps)
+        return
+    yield _get_edge_lonlat_samples(area_to_cover, sample_steps)
+
+
+def _iter_dense_lonlat_samples(area_to_cover: AreaDefinition):
+    """Yield full destination lon/lat coverage in row chunks."""
+    row_block_size = max(1, MAX_POINTS_PER_CHUNK // area_to_cover.width)
+    for row_start in range(0, area_to_cover.height, row_block_size):
+        row_stop = min(area_to_cover.height, row_start + row_block_size)
+        yield area_to_cover.get_lonlats(
+            data_slice=(slice(row_start, row_stop), slice(None)),
+            dtype=np.float32,
+        )
+
+
+def _get_grid_lonlat_samples(area_to_cover: AreaDefinition, sample_steps: int):
+    """Return one evenly spaced interior destination sample grid."""
+    sample_rows = _get_sample_indices(area_to_cover.height, sample_steps)
+    sample_cols = _get_sample_indices(area_to_cover.width, sample_steps)
+    return area_to_cover.get_lonlats(
+        data_slice=(sample_rows[:, None], sample_cols[None, :]),
+        dtype=np.float32,
+    )
+
+
+def _get_edge_lonlat_samples(area_to_cover: AreaDefinition, sample_steps: int | None):
+    """Return perimeter destination samples for edge-only sampling mode."""
+    if area_to_cover.is_geostationary:
+        # Use limb-aware geostationary boundary sampling instead of raw array
+        # corners/edges. Projection corners may be off-earth and can under-cover.
+        if sample_steps is None:
+            nb_points = max(4, area_to_cover.width * 2 + area_to_cover.height * 2 - 4)
+        else:
+            nb_points = max(4, sample_steps * 4)
+        return get_geostationary_bounding_box_in_lonlats(
+            area_to_cover,
+            nb_points=nb_points,
+        )
+    sample_rows = _get_sample_indices(area_to_cover.height, sample_steps)
+    sample_cols = _get_sample_indices(area_to_cover.width, sample_steps)
+    top_rows = np.zeros(sample_cols.size, dtype=np.int64)
+    top_cols = sample_cols
+    right_rows = sample_rows[1:]
+    right_cols = np.full(right_rows.size, area_to_cover.width - 1, dtype=np.int64)
+    bottom_cols = sample_cols[-2::-1]
+    bottom_rows = np.full(bottom_cols.size, area_to_cover.height - 1, dtype=np.int64)
+    left_rows = sample_rows[-2:0:-1]
+    left_cols = np.zeros(left_rows.size, dtype=np.int64)
+    edge_rows = np.concatenate((top_rows, right_rows, bottom_rows, left_rows))
+    edge_cols = np.concatenate((top_cols, right_cols, bottom_cols, left_cols))
+    return area_to_cover.get_lonlats(data_slice=(edge_rows, edge_cols), dtype=np.float32)
+
+
+def _get_sample_indices(axis_size: int, sample_steps: int | None):
+    """Return evenly spaced integer sample indices including both endpoints."""
+    if sample_steps is None or sample_steps >= axis_size:
+        return np.arange(axis_size, dtype=np.int64)
+    return (np.arange(sample_steps, dtype=np.int64) * (axis_size - 1)) // (sample_steps - 1)
+
+
 def _make_slice_divisible(sli, max_size, factor=2):
     """Make the given slice even in size."""
     rem = (sli.stop - sli.start) % factor
@@ -121,6 +315,7 @@ def _make_slice_divisible(sli, max_size, factor=2):
 
 
 def _ensure_integer_slice(sli):
+    """Round slice bounds outward to conservatively preserve target coverage."""
     start = sli.start
     stop = sli.stop
     step = sli.step
diff --git a/pyresample/geometry.py b/pyresample/geometry.py
index be3f87e0..1c75a789 100644
--- a/pyresample/geometry.py
+++ b/pyresample/geometry.py
@@ -68,6 +68,8 @@
 
 logger = getLogger(__name__)
 
+DEFAULT_AREA_SLICE_SAMPLE_STEPS = 21
+
 
 class DimensionError(ValueError):
     """Wrap ValueError."""
@@ -619,7 +621,13 @@ def overlap_rate(self, other):
         inter_area = get_polygon_area(self.intersection(other))
         return inter_area / other_area
 
-    def get_area_slices(self, area_to_cover):
+    def get_area_slices(
+            self,
+            area_to_cover,
+            shape_divisible_by=None,
+            sample_steps=DEFAULT_AREA_SLICE_SAMPLE_STEPS,
+            sample_grid=False,
+    ):
         """Compute the slice to read based on an `area_to_cover`."""
         raise NotImplementedError
 
@@ -2651,10 +2659,32 @@ def proj4_string(self):
                       "instead.", DeprecationWarning, stacklevel=2)
         return proj4_dict_to_str(self.proj_dict)
 
-    def get_area_slices(self, area_to_cover, shape_divisible_by=None):
-        """Compute the slice to read based on an `area_to_cover`."""
+    def get_area_slices(
+            self,
+            area_to_cover,
+            shape_divisible_by=None,
+            sample_steps=DEFAULT_AREA_SLICE_SAMPLE_STEPS,
+            sample_grid=False,
+    ):
+        """Compute the slice to read based on an `area_to_cover`.
+
+        Args:
+            area_to_cover: Destination area to crop around.
+            shape_divisible_by: Optional factor to make the resulting shape divisible by.
+            sample_steps: Number of edge/grid samples used for geostationary cross-projection coverage checks.
+                Use ``None`` or ``<= 0`` for ``ALL`` mode: all edge points when ``sample_grid=False``
+                and all destination pixels when ``sample_grid=True``.
+            sample_grid: If ``True``, sample an interior grid with ``sample_steps`` density.
+                If ``False``, sample edge points only.
+        """
         from .future.geometry._subset import get_area_slices
-        return get_area_slices(self, area_to_cover, shape_divisible_by)
+        return get_area_slices(
+            self,
+            area_to_cover,
+            shape_divisible_by,
+            sample_steps,
+            sample_grid,
+        )
 
     def crop_around(self, other_area):
         """Crop this area around `other_area`."""
diff --git a/pyresample/test/test_geometry/test_area.py b/pyresample/test/test_geometry/test_area.py
index f55b5235..075fbdad 100644
--- a/pyresample/test/test_geometry/test_area.py
+++ b/pyresample/test/test_geometry/test_area.py
@@ -29,6 +29,7 @@
 from pyresample.future.geometry import AreaDefinition, SwathDefinition
 from pyresample.future.geometry.area import ignore_pyproj_proj_warnings
 from pyresample.future.geometry.base import get_array_hashable
+from pyresample.geometry import DEFAULT_AREA_SLICE_SAMPLE_STEPS
 from pyresample.geometry import AreaDefinition as LegacyAreaDefinition
 from pyresample.geometry import DynamicAreaDefinition
 from pyresample.test.utils import assert_future_geometry
@@ -1520,6 +1521,92 @@ def test_enclose_areas(create_test_area):
         enclose_areas()
 
 
+def _assert_slice_contains_destination_coverage(slice_x, slice_y, valid_cols, valid_rows):
+    assert valid_cols.size > 0
+    assert valid_rows.size > 0
+    assert slice_x.start <= valid_cols.min()
+    assert slice_y.start <= valid_rows.min()
+    assert slice_x.stop > valid_cols.max()
+    assert slice_y.stop > valid_rows.max()
+
+
+def _assert_default_matches_edge_sampling(source_area, area_to_cover):
+    default_slice = source_area.get_area_slices(area_to_cover)
+    edge_slice = source_area.get_area_slices(
+        area_to_cover,
+        sample_steps=DEFAULT_AREA_SLICE_SAMPLE_STEPS,
+        sample_grid=False,
+    )
+    assert default_slice == edge_slice
+    return default_slice
+
+
+def _assert_default_matches_edge_and_is_within_dense(source_area, area_to_cover):
+    default_slice = _assert_default_matches_edge_sampling(source_area, area_to_cover)
+    dense_slice = source_area.get_area_slices(
+        area_to_cover,
+        sample_steps=None,
+        sample_grid=True,
+    )
+    edge_all_slice = source_area.get_area_slices(
+        area_to_cover,
+        sample_steps=None,
+        sample_grid=False,
+    )
+    assert default_slice[0].start >= edge_all_slice[0].start
+    assert default_slice[0].stop <= edge_all_slice[0].stop
+    assert default_slice[1].start >= edge_all_slice[1].start
+    assert default_slice[1].stop <= edge_all_slice[1].stop
+    assert edge_all_slice[0].start >= dense_slice[0].start
+    assert edge_all_slice[0].stop <= dense_slice[0].stop
+    assert edge_all_slice[1].start >= dense_slice[1].start
+    assert edge_all_slice[1].stop <= dense_slice[1].stop
+
+
+@pytest.fixture
+def geos_cross_projection_areas():
+    source_area = LegacyAreaDefinition(
+        "src",
+        "src",
+        "src",
+        {
+            "proj": "geos",
+            "a": 6378169.0,
+            "b": 6356583.8,
+            "h": 35785831.0,
+            "lon_0": 0.0,
+            "units": "m",
+        },
+        3712,
+        1392,
+        (5568748.0, 5568748.0, -5568748.0, 1392187.0),
+    )
+    area_to_cover = LegacyAreaDefinition(
+        "dst",
+        "dst",
+        "dst",
+        {"proj": "latlong", "datum": "WGS84"},
+        768,
+        768,
+        (7.456086060029671, 43.98125198600542, 33.461915849571966, 59.24271064133026),
+    )
+    return source_area, area_to_cover
+
+
+@pytest.fixture
+def geos_cross_projection_valid_indices(geos_cross_projection_areas):
+    source_area, area_to_cover = geos_cross_projection_areas
+    destination_lons, destination_lats = area_to_cover.get_lonlats()
+    source_cols, source_rows = source_area.get_array_indices_from_lonlat(
+        destination_lons,
+        destination_lats,
+    )
+    valid = ~np.ma.getmaskarray(source_cols) & ~np.ma.getmaskarray(source_rows)
+    valid_cols = np.asarray(source_cols[valid], dtype=np.int64)
+    valid_rows = np.asarray(source_rows[valid], dtype=np.int64)
+    return valid_cols, valid_rows
+
+
 class TestAreaDefGetAreaSlices:
     """Test AreaDefinition's get_area_slices."""
 
@@ -1535,8 +1622,6 @@ def test_get_area_slices_geos_subset(self, geos_src_area, create_test_area):
                          area_def.area_extent[2] - 10000,
                          area_def.area_extent[3] - 10000))
         slice_x, slice_y = area_def.get_area_slices(area_to_cover)
-        assert isinstance(slice_x.start, int)
-        assert isinstance(slice_y.start, int)
         assert slice(3, 3709, None) == slice_x
         assert slice(3, 3709, None) == slice_y
 
@@ -1554,10 +1639,8 @@ def test_get_area_slices_geos_similar(self, geos_src_area, create_test_area):
             y_size,
             area_extent)
         slice_x, slice_y = area_def.get_area_slices(area_to_cover)
-        assert isinstance(slice_x.start, int)
-        assert isinstance(slice_y.start, int)
         assert slice(46, 3667, None) == slice_x
-        assert slice(56, 3659, None) == slice_y
+        assert slice(52, 3663, None) == slice_y
 
     def test_get_area_slices_geos_stereographic(self, geos_src_area, create_test_area):
         """Test slicing with a geos area and polar stereographic area."""
@@ -1573,6 +1656,157 @@ def test_get_area_slices_geos_stereographic(self, geos_src_area, create_test_are
         assert slice_x == slice(1610, 2343)
         assert slice_y == slice(158, 515, None)
 
+    def test_get_area_slices_geos_cross_projection_contains_destination_coverage(
+            self,
+            geos_cross_projection_areas,
+            geos_cross_projection_valid_indices,
+    ):
+        """Test geos cross-projection slicing includes all destination-covered source pixels in dense mode."""
+        source_area, area_to_cover = geos_cross_projection_areas
+        valid_cols, valid_rows = geos_cross_projection_valid_indices
+        slice_x, slice_y = source_area.get_area_slices(area_to_cover, sample_steps=None, sample_grid=True)
+
+        _assert_slice_contains_destination_coverage(slice_x, slice_y, valid_cols, valid_rows)
+
+    def test_get_area_slices_geos_cross_projection_shape_divisible(
+            self,
+            geos_cross_projection_areas,
+            geos_cross_projection_valid_indices,
+    ):
+        """Test geos cross-projection shape-divisible slicing keeps coverage in dense mode."""
+        source_area, area_to_cover = geos_cross_projection_areas
+        valid_cols, valid_rows = geos_cross_projection_valid_indices
+        slice_x, slice_y = source_area.get_area_slices(
+            area_to_cover,
+            shape_divisible_by=2,
+            sample_steps=None,
+            sample_grid=True,
+        )
+
+        _assert_slice_contains_destination_coverage(slice_x, slice_y, valid_cols, valid_rows)
+        assert (slice_x.stop - slice_x.start) % 2 == 0
+        assert (slice_y.stop - slice_y.start) % 2 == 0
+
+    def test_get_area_slices_geos_cross_projection_default_matches_edge_sampling(
+            self,
+            geos_cross_projection_areas,
+    ):
+        """Test default slicing uses edge-only sampling with default step count."""
+        source_area, area_to_cover = geos_cross_projection_areas
+        _assert_default_matches_edge_sampling(source_area, area_to_cover)
+
+    def test_get_area_slices_geos_cross_projection_sample_grid_expands_or_matches_edge(
+            self,
+            geos_cross_projection_areas,
+    ):
+        """Test interior grid sampling is a superset of edge-only sampling."""
+        source_area, area_to_cover = geos_cross_projection_areas
+        edge_slice_x, edge_slice_y = source_area.get_area_slices(
+            area_to_cover,
+            sample_steps=DEFAULT_AREA_SLICE_SAMPLE_STEPS,
+            sample_grid=False,
+        )
+        grid_slice_x, grid_slice_y = source_area.get_area_slices(
+            area_to_cover,
+            sample_steps=DEFAULT_AREA_SLICE_SAMPLE_STEPS,
+            sample_grid=True,
+        )
+
+        assert grid_slice_x.start <= edge_slice_x.start
+        assert grid_slice_x.stop >= edge_slice_x.stop
+        assert grid_slice_y.start <= edge_slice_y.start
+        assert grid_slice_y.stop >= edge_slice_y.stop
+
+    @pytest.mark.parametrize("sample_steps", [None, 0, -1, False])
+    def test_get_area_slices_geos_cross_projection_non_positive_steps_match_all_edge_mode(
+            self,
+            geos_cross_projection_areas,
+            sample_steps,
+    ):
+        """Test non-positive step values mean ALL edge points when grid sampling is disabled."""
+        source_area, area_to_cover = geos_cross_projection_areas
+        edge_all_slice_x, edge_all_slice_y = source_area.get_area_slices(
+            area_to_cover,
+            sample_steps=None,
+            sample_grid=False,
+        )
+        mode_slice_x, mode_slice_y = source_area.get_area_slices(
+            area_to_cover,
+            sample_steps=sample_steps,
+            sample_grid=False,
+        )
+
+        assert (mode_slice_x, mode_slice_y) == (edge_all_slice_x, edge_all_slice_y)
+
+    @pytest.mark.parametrize("sample_steps", [None, 0, -1, False])
+    def test_get_area_slices_geos_cross_projection_non_positive_steps_match_dense_grid_mode(
+            self,
+            geos_cross_projection_areas,
+            sample_steps,
+    ):
+        """Test non-positive step values mean dense sampling when grid sampling is enabled."""
+        source_area, area_to_cover = geos_cross_projection_areas
+        dense_slice_x, dense_slice_y = source_area.get_area_slices(
+            area_to_cover,
+            sample_steps=None,
+            sample_grid=True,
+        )
+        mode_slice_x, mode_slice_y = source_area.get_area_slices(
+            area_to_cover,
+            sample_steps=sample_steps,
+            sample_grid=True,
+        )
+
+        assert (mode_slice_x, mode_slice_y) == (dense_slice_x, dense_slice_y)
+
+    def test_get_area_slices_geos_cross_projection_step_1_raises_value_error(self, geos_cross_projection_areas):
+        """Test ``sample_steps=1`` is rejected as an invalid sampled-mode value."""
+        source_area, area_to_cover = geos_cross_projection_areas
+        with pytest.raises(ValueError, match="sample_steps=1 is not supported"):
+            source_area.get_area_slices(
+                area_to_cover,
+                sample_steps=1,
+                sample_grid=False,
+            )
+
+    def test_get_area_slices_geos_cross_projection_short_circuits_boundary_when_coverage_exists(
+            self,
+            geos_cross_projection_areas,
+            monkeypatch,
+    ):
+        """Test geos cross-projection slicing returns from sampled coverage before boundary fallback."""
+        from pyresample.future.geometry import _subset
+        source_area, area_to_cover = geos_cross_projection_areas
+
+        def _fail_if_called(_):
+            raise AssertionError("Boundary fallback should not run when sampled coverage exists.")
+
+        monkeypatch.setattr(_subset, "_get_area_boundary", _fail_if_called)
+        source_area.get_area_slices(area_to_cover)
+
+    def test_get_area_slices_geos_cross_projection_falls_back_to_boundary_when_coverage_missing(
+            self,
+            geos_cross_projection_areas,
+            monkeypatch,
+    ):
+        """Test geos cross-projection slicing can still use boundary fallback."""
+        from pyresample.future.geometry import _subset
+        source_area, area_to_cover = geos_cross_projection_areas
+
+        boundary_calls = 0
+        orig_boundary = _subset._get_area_boundary
+
+        def _boundary_spy(area):
+            nonlocal boundary_calls
+            boundary_calls += 1
+            return orig_boundary(area)
+
+        monkeypatch.setattr(_subset, "_get_covered_source_slices", lambda *args, **kwargs: None)
+        monkeypatch.setattr(_subset, "_get_area_boundary", _boundary_spy)
+        source_area.get_area_slices(area_to_cover)
+
+        assert boundary_calls >= 2
+
     def test_get_area_slices_geos_flipped_xy(self, geos_src_area, create_test_area):
         """Test slicing with two geos areas but one has flipped x/y dimensions."""
         area_def = geos_src_area
@@ -1602,12 +1836,7 @@ def test_get_area_slices_geos_epsg_lonlat(self, geos_src_area, create_test_area)
                 8192,
                 4096,
                 (-180.0, -90.0, 180.0, 90.0))
-
-            slice_x, slice_y = area_def.get_area_slices(area_to_cover)
-            assert isinstance(slice_x.start, int)
-            assert isinstance(slice_y.start, int)
-            assert slice_x == slice(46, 3667, None)
-            assert slice_y == slice(56, 3659, None)
+            _assert_default_matches_edge_sampling(area_def, area_to_cover)
 
     def test_get_area_slices_nongeos(self, create_test_area):
         """Check area slicing for non-geos projections."""
@@ -1663,15 +1892,14 @@ def test_non_geos_can_be_cropped(self, create_test_area):
         assert slice_y == slice(9261, 10980)
 
     def test_area_to_cover_all_nan_bounds(self, geos_src_area, create_test_area):
-        """Check area slicing when the target doesn't have a valid boundary."""
+        """Check dense slicing still works when target boundary extraction is problematic."""
         area_def = geos_src_area
         # An area that is a subset of the original one
         area_to_cover = create_test_area(
             {"proj": "moll"},
             1000, 1000,
             area_extent=(-18000000.0, -9000000.0, 18000000.0, 9000000.0))
-        with pytest.raises(NotImplementedError):
-            area_def.get_area_slices(area_to_cover)
+        _assert_default_matches_edge_and_is_within_dense(area_def, area_to_cover)
 
     @pytest.mark.parametrize("cache_slices", [False, True])
     def test_area_slices_caching(self, create_test_area, tmp_path, cache_slices):