diff --git a/docs/source/howtos/configuration.rst b/docs/source/howtos/configuration.rst
index b7aa6e3f1..6a79acf10 100644
--- a/docs/source/howtos/configuration.rst
+++ b/docs/source/howtos/configuration.rst
@@ -76,6 +76,67 @@ Or for specific blocks of code:
 Similarly, if you need to access one of the values you can
 use the ``pyresample.config.get`` method.
 
+Cache Directory
+^^^^^^^^^^^^^^^
+
+* **Environment variable**: ``PYRESAMPLE_CACHE_DIR``
+* **YAML/Config Key**: ``cache_dir``
+* **Default**: See below
+
+Directory where any files cached by Pyresample will be stored. This
+directory is not necessarily cleared out by Pyresample, but is rarely used
+without explicitly being enabled by the user. This
+defaults to a different path depending on your operating system following
+the `platformdirs <https://github.com/platformdirs/platformdirs#example-output>`_
+"user cache dir".
+
+.. note::
+
+   Some resampling algorithms provide caching functionality when the user
+   provides a directory to cache to. These resamplers do not currently use this
+   configuration option.
+
+.. _config_cache_sensor_angles_setting:
+
+Cache Geometry Slices
+^^^^^^^^^^^^^^^^^^^^^
+
+* **Environment variable**: ``PYRESAMPLE_CACHE_GEOMETRY_SLICES``
+* **YAML/Config Key**: ``cache_geometry_slices``
+* **Default**: ``False``
+
+Whether or not generated slices for geometry objects are cached to disk.
+These slices are used in various parts of Pyresample like
+cropping or overlap calculations including those performed in some resampling
+algorithms. At the time of writing this is only performed on
+``AreaDefinition`` objects through their
+:meth:`~pyresample.geometry.AreaDefinition.get_area_slices` method.
+Slices are stored in ``cache_dir`` (see above).
+Unlike other caching performed in Pyresample where potentially large arrays
+are cached, this option saves a pair of ``slice`` objects that consist of
+only 3 integers each. This makes the amount of space used in the cache very
+small for many cached results.
+
+The slicing operations in Pyresample typically involve finding the intersection
+between two geometries. This requires generating bounding polygons for the
+geometries and doing polygon intersection calculations that can handle
+projection anti-meridians. At the time of writing these calculations can take
+as long as 15 seconds depending on number of vertices used in the bounding
+polygons. One use case for these slices is reducing input data to only the
+overlap of the target area. This can be done before or during resampling as
+part of the algorithm or as part of a third-party resampling interface
+(ex. Satpy). In the future as optimizations are made to the polygon
+intersection logic this caching option should hopefully not be needed.
+
+When setting this as an environment variable, this should be set with the
+string equivalent of the Python boolean values ``="True"`` or ``="False"``.
+
+.. warning::
+
+    This caching does not limit the number of entries nor does it expire old
+    entries. It is up to the user to manage the contents of the cache
+    directory.
+
 Feature Flags
 -------------
 
diff --git a/pyresample/_caching.py b/pyresample/_caching.py
new file mode 100644
index 000000000..08b326e0c
--- /dev/null
+++ b/pyresample/_caching.py
@@ -0,0 +1,128 @@
+"""Various tools for caching.
+
+These tools are rarely needed by users and are used where they make sense
+throughout pyresample.
+
+"""
+from __future__ import annotations
+
+import hashlib
+import json
+import os
+import warnings
+from functools import update_wrapper
+from glob import glob
+from pathlib import Path
+from typing import Any, Callable
+
+import pyresample
+
+
+class JSONCacheHelper:
+    """Decorator class to cache results to a JSON file on-disk."""
+
+    def __init__(
+            self,
+            func: Callable,
+            cache_config_key: str,
+            cache_version: int = 1,
+    ):
+        self._callable = func
+        self._cache_config_key = cache_config_key
+        self._cache_version = cache_version
+        self._uncacheable_arg_type_names = ("",)
+
+    @staticmethod
+    def cache_clear(cache_dir: str | None = None):
+        """Remove all on-disk files associated with this function.
+
+        Intended to mimic the :func:`functools.cache` behavior.
+        """
+        cache_path = _get_cache_dir_from_config(cache_dir=cache_dir, cache_version="*")
+        for json_file in glob(str(cache_path / "*.json")):
+            os.remove(json_file)
+
+    def __call__(self, *args):
+        """Call decorated function and cache the result to JSON."""
+        should_cache = pyresample.config.get(self._cache_config_key, False)
+        if not should_cache:
+            return self._callable(*args)
+
+        try:
+            arg_hash = _hash_args(args)
+        except TypeError as err:
+            warnings.warn("Cannot cache function due to unhashable argument: " + str(err),
+                          stacklevel=2)
+            return self._callable(*args)
+
+        return self._run_and_cache(arg_hash, args)
+
+    def _run_and_cache(self, arg_hash: str, args: tuple[Any]) -> Any:
+        base_cache_dir = _get_cache_dir_from_config(cache_version=self._cache_version)
+        json_path = base_cache_dir / f"{arg_hash}.json"
+        if not json_path.is_file():
+            res = self._callable(*args)
+            json_path.parent.mkdir(exist_ok=True)
+            with open(json_path, "w") as json_cache:
+                json.dump(res, json_cache, cls=_JSONEncoderWithSlice)
+
+        # for consistency, always load the cached result
+        with open(json_path, "r") as json_cache:
+            res = json.load(json_cache, object_hook=_object_hook)
+        return res
+
+
+def _get_cache_dir_from_config(cache_dir: str | None = None, cache_version: int | str = 1) -> Path:
+    cache_dir = cache_dir or pyresample.config.get("cache_dir")
+    if cache_dir is None:
+        raise RuntimeError("Can't use JSON caching. No 'cache_dir' configured.")
+    subdir = f"geometry_slices_v{cache_version}"
+    return Path(cache_dir) / subdir
+
+
+def _hash_args(args: tuple[Any]) -> str:
+    from pyresample.future.geometry import AreaDefinition, SwathDefinition
+    from pyresample.geometry import AreaDefinition as LegacyAreaDefinition
+    from pyresample.geometry import SwathDefinition as LegacySwathDefinition
+
+    hashable_args = []
+    for arg in args:
+        if isinstance(arg, (SwathDefinition, LegacySwathDefinition)):
+            raise TypeError(f"Unhashable type ({type(arg)})")
+        if isinstance(arg, (AreaDefinition, LegacyAreaDefinition)):
+            arg = hash(arg)
+        hashable_args.append(arg)
+    arg_hash = hashlib.sha1()  # nosec
+    arg_hash.update(json.dumps(tuple(hashable_args)).encode("utf8"))
+    return arg_hash.hexdigest()
+
+
+class _JSONEncoderWithSlice(json.JSONEncoder):
+    def default(self, obj: Any) -> Any:
+        if isinstance(obj, slice):
+            return {"__slice__": True, "start": obj.start, "stop": obj.stop, "step": obj.step}
+        return super().default(obj)
+
+
+def _object_hook(obj: object) -> Any:
+    if isinstance(obj, dict) and obj.get("__slice__", False):
+        return slice(obj["start"], obj["stop"], obj["step"])
+    return obj
+
+
+def cache_to_json_if(cache_config_key: str) -> Callable:
+    """Decorate a function and cache the results to a JSON file on disk.
+
+    This caching only happens if the ``pyresample.config`` boolean value for
+    the provided key is ``True`` as well as some other conditions. See
+    :class:`JSONCacheHelper` for more information. Most importantly this
+    decorator does not limit how many items can be cached and does not clear
+    out old entries. It is up to the user to manage the size of the cache.
+
+    """
+    def _decorator(func: Callable) -> Callable:
+        zarr_cacher = JSONCacheHelper(func, cache_config_key)
+        wrapper = update_wrapper(zarr_cacher, func)
+        return wrapper
+
+    return _decorator
diff --git a/pyresample/_config.py b/pyresample/_config.py
index ede6b33db..9a4a2b1ae 100644
--- a/pyresample/_config.py
+++ b/pyresample/_config.py
@@ -21,7 +21,6 @@
 from donfig import Config
 
 BASE_PATH = os.path.dirname(os.path.realpath(__file__))
-# FIXME: Use package_resources?
 PACKAGE_CONFIG_PATH = os.path.join(BASE_PATH, 'etc')
 
 _user_config_dir = platformdirs.user_config_dir("pyresample", "pytroll")
@@ -36,6 +35,8 @@
 config = Config(
     "pyresample",
     defaults=[{
+        "cache_dir": platformdirs.user_cache_dir("pyresample", "pytroll"),
+        "cache_geometry_slices": False,
         "features": {
             "future_geometries": False,
         },
diff --git a/pyresample/future/geometry/_subset.py b/pyresample/future/geometry/_subset.py
new file mode 100644
index 000000000..bb04cf69b
--- /dev/null
+++ b/pyresample/future/geometry/_subset.py
@@ -0,0 +1,131 @@
+"""Functions and tools for subsetting a geometry object."""
+from __future__ import annotations
+
+import math
+from typing import TYPE_CHECKING, Any
+
+import numpy as np
+
+# 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.utils import check_slice_orientation
+
+if TYPE_CHECKING:
+    from pyresample import AreaDefinition
+
+
+@cache_to_json_if("cache_geometry_slices")
+def get_area_slices(
+        src_area: AreaDefinition,
+        area_to_cover: AreaDefinition,
+        shape_divisible_by: int | None,
+) -> tuple[slice, slice]:
+    """Compute the slice to read based on an `area_to_cover`."""
+    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:
+        logger.debug('Projections for data and slice areas are identical: %s',
+                     proj_def_to_cover)
+        # Get slice parameters
+        xstart, xstop, ystart, ystop = _get_slice_starts_stops(src_area, area_to_cover)
+
+        x_slice = check_slice_orientation(slice(xstart, xstop))
+        y_slice = check_slice_orientation(slice(ystart, ystop))
+        x_slice = _ensure_integer_slice(x_slice)
+        y_slice = _ensure_integer_slice(y_slice)
+        return x_slice, y_slice
+
+    data_boundary = _get_area_boundary(src_area)
+    area_boundary = _get_area_boundary(area_to_cover)
+    intersection = data_boundary.contour_poly.intersection(
+        area_boundary.contour_poly)
+    if intersection is None:
+        logger.debug('Cannot determine appropriate slicing. '
+                     "Data and projection area do not overlap.")
+        raise NotImplementedError
+    x, y = src_area.get_array_indices_from_lonlat(
+        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)
+    x_slice = _ensure_integer_slice(x_slice)
+    y_slice = _ensure_integer_slice(y_slice)
+    if shape_divisible_by is not None:
+        x_slice = _make_slice_divisible(x_slice, src_area.width,
+                                        factor=shape_divisible_by)
+        y_slice = _make_slice_divisible(y_slice, src_area.height,
+                                        factor=shape_divisible_by)
+
+    return (check_slice_orientation(x_slice),
+            check_slice_orientation(y_slice))
+
+
+def _is_area_like(area_obj: Any) -> bool:
+    return hasattr(area_obj, "crs") and hasattr(area_obj, "area_extent")
+
+
+def _get_slice_starts_stops(src_area, area_to_cover):
+    """Get x and y start and stop points for slicing."""
+    llx, lly, urx, ury = area_to_cover.area_extent
+    x, y = src_area.get_array_coordinates_from_projection_coordinates([llx, urx], [lly, ury])
+
+    # we use `round` because we want the *exterior* of the pixels to contain the area_to_cover's area extent.
+    if (src_area.area_extent[0] > src_area.area_extent[2]) ^ (llx > urx):
+        xstart = max(0, round(x[1]))
+        xstop = min(src_area.width, round(x[0]) + 1)
+    else:
+        xstart = max(0, round(x[0]))
+        xstop = min(src_area.width, round(x[1]) + 1)
+    if (src_area.area_extent[1] > src_area.area_extent[3]) ^ (lly > ury):
+        ystart = max(0, round(y[0]))
+        ystop = min(src_area.height, round(y[1]) + 1)
+    else:
+        ystart = max(0, round(y[1]))
+        ystop = min(src_area.height, round(y[0]) + 1)
+
+    return xstart, xstop, ystart, ystop
+
+
+def _get_area_boundary(area_to_cover: AreaDefinition) -> Boundary:
+    try:
+        if area_to_cover.is_geostationary:
+            return Boundary(*get_geostationary_bounding_box_in_lonlats(area_to_cover))
+        boundary_shape = max(max(*area_to_cover.shape) // 100 + 1, 3)
+        return area_to_cover.boundary(frequency=boundary_shape, force_clockwise=True)
+    except ValueError as err:
+        raise NotImplementedError("Can't determine boundary of area to cover") from err
+
+
+def _make_slice_divisible(sli, max_size, factor=2):
+    """Make the given slice even in size."""
+    rem = (sli.stop - sli.start) % factor
+    if rem != 0:
+        adj = factor - rem
+        if sli.stop + 1 + rem < max_size:
+            sli = slice(sli.start, sli.stop + adj)
+        elif sli.start > 0:
+            sli = slice(sli.start - adj, sli.stop)
+        else:
+            sli = slice(sli.start, sli.stop - rem)
+
+    return sli
+
+
+def _ensure_integer_slice(sli):
+    start = sli.start
+    stop = sli.stop
+    step = sli.step
+    return slice(
+        math.floor(start) if start is not None else None,
+        math.ceil(stop) if stop is not None else None,
+        math.floor(step) if step is not None else None
+    )
diff --git a/pyresample/geometry.py b/pyresample/geometry.py
index 91baa1b23..2c2b52254 100644
--- a/pyresample/geometry.py
+++ b/pyresample/geometry.py
@@ -37,8 +37,8 @@
 from pyresample import CHUNK_SIZE
 from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj_MP
 from pyresample.area_config import create_area_def
-from pyresample.boundary import Boundary, SimpleBoundary
-from pyresample.utils import check_slice_orientation, load_cf_area
+from pyresample.boundary import SimpleBoundary
+from pyresample.utils import load_cf_area
 from pyresample.utils.proj4 import (
     get_geodetic_crs_with_no_datum_shift,
     get_geostationary_height,
@@ -2579,70 +2579,10 @@ def proj4_string(self):
                       "instead.", DeprecationWarning, stacklevel=2)
         return proj4_dict_to_str(self.proj_dict)
 
-    def _get_slice_starts_stops(self, area_to_cover):
-        """Get x and y start and stop points for slicing."""
-        llx, lly, urx, ury = area_to_cover.area_extent
-        x, y = self.get_array_coordinates_from_projection_coordinates([llx, urx], [lly, ury])
-
-        # we use `round` because we want the *exterior* of the pixels to contain the area_to_cover's area extent.
-        if (self.area_extent[0] > self.area_extent[2]) ^ (llx > urx):
-            xstart = max(0, round(x[1]))
-            xstop = min(self.width, round(x[0]) + 1)
-        else:
-            xstart = max(0, round(x[0]))
-            xstop = min(self.width, round(x[1]) + 1)
-        if (self.area_extent[1] > self.area_extent[3]) ^ (lly > ury):
-            ystart = max(0, round(y[0]))
-            ystop = min(self.height, round(y[1]) + 1)
-        else:
-            ystart = max(0, round(y[1]))
-            ystop = min(self.height, round(y[0]) + 1)
-
-        return xstart, xstop, ystart, ystop
-
     def get_area_slices(self, area_to_cover, shape_divisible_by=None):
         """Compute the slice to read based on an `area_to_cover`."""
-        if not isinstance(area_to_cover, AreaDefinition):
-            raise NotImplementedError('Only AreaDefinitions can be used')
-
-        # Intersection only required for two different projections
-        proj_def_to_cover = area_to_cover.crs
-        proj_def = self.crs
-        if proj_def_to_cover == proj_def:
-            logger.debug('Projections for data and slice areas are'
-                         ' identical: %s',
-                         proj_def_to_cover)
-            # Get slice parameters
-            xstart, xstop, ystart, ystop = self._get_slice_starts_stops(area_to_cover)
-
-            x_slice = check_slice_orientation(slice(xstart, xstop))
-            y_slice = check_slice_orientation(slice(ystart, ystop))
-            x_slice = _ensure_integer_slice(x_slice)
-            y_slice = _ensure_integer_slice(y_slice)
-            return x_slice, y_slice
-
-        data_boundary = _get_area_boundary(self)
-        area_boundary = _get_area_boundary(area_to_cover)
-        intersection = data_boundary.contour_poly.intersection(
-            area_boundary.contour_poly)
-        if intersection is None:
-            logger.debug('Cannot determine appropriate slicing. '
-                         "Data and projection area do not overlap.")
-            raise NotImplementedError
-        x, y = self.get_array_indices_from_lonlat(
-            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)
-        x_slice = _ensure_integer_slice(x_slice)
-        y_slice = _ensure_integer_slice(y_slice)
-        if shape_divisible_by is not None:
-            x_slice = _make_slice_divisible(x_slice, self.width,
-                                            factor=shape_divisible_by)
-            y_slice = _make_slice_divisible(y_slice, self.height,
-                                            factor=shape_divisible_by)
-
-        return (check_slice_orientation(x_slice),
-                check_slice_orientation(y_slice))
+        from .future.geometry._subset import get_area_slices
+        return get_area_slices(self, area_to_cover, shape_divisible_by)
 
     def crop_around(self, other_area):
         """Crop this area around `other_area`."""
@@ -2743,42 +2683,6 @@ def geocentric_resolution(self, ellps='WGS84', radius=None):
         return res
 
 
-def _get_area_boundary(area_to_cover: AreaDefinition) -> Boundary:
-    try:
-        if area_to_cover.is_geostationary:
-            return Boundary(*get_geostationary_bounding_box_in_lonlats(area_to_cover))
-        boundary_shape = max(max(*area_to_cover.shape) // 100 + 1, 3)
-        return area_to_cover.boundary(frequency=boundary_shape, force_clockwise=True)
-    except ValueError as err:
-        raise NotImplementedError("Can't determine boundary of area to cover") from err
-
-
-def _make_slice_divisible(sli, max_size, factor=2):
-    """Make the given slice even in size."""
-    rem = (sli.stop - sli.start) % factor
-    if rem != 0:
-        adj = factor - rem
-        if sli.stop + 1 + rem < max_size:
-            sli = slice(sli.start, sli.stop + adj)
-        elif sli.start > 0:
-            sli = slice(sli.start - adj, sli.stop)
-        else:
-            sli = slice(sli.start, sli.stop - rem)
-
-    return sli
-
-
-def _ensure_integer_slice(sli):
-    start = sli.start
-    stop = sli.stop
-    step = sli.step
-    return slice(
-        math.floor(start) if start is not None else None,
-        math.ceil(stop) if stop is not None else None,
-        math.floor(step) if step is not None else None
-    )
-
-
 def get_geostationary_angle_extent(geos_area):
     """Get the max earth (vs space) viewing angles in x and y."""
     # get some projection parameters
diff --git a/pyresample/test/conftest.py b/pyresample/test/conftest.py
index 992fa5a02..2643caad7 100644
--- a/pyresample/test/conftest.py
+++ b/pyresample/test/conftest.py
@@ -42,6 +42,7 @@
 def reset_pyresample_config(tmpdir):
     """Set pyresample config to logical defaults for tests."""
     test_config = {
+        "cache_geometry_slices": False,
         "features": {
             "future_geometries": False,
         },
diff --git a/pyresample/test/test_geometry/test_area.py b/pyresample/test/test_geometry/test_area.py
index 4ed0b283c..d1aeda7ee 100644
--- a/pyresample/test/test_geometry/test_area.py
+++ b/pyresample/test/test_geometry/test_area.py
@@ -15,6 +15,7 @@
 """Test AreaDefinition objects."""
 import io
 import sys
+from glob import glob
 from unittest.mock import MagicMock, patch
 
 import dask.array as da
@@ -879,41 +880,31 @@ def test_get_array_indices_from_lonlat(self, create_test_area):
         assert (x__.data == x_expects).all()
         assert (y__.data == y_expects).all()
 
-    def test_get_slice_starts_stops(self, create_test_area):
+    @pytest.mark.parametrize(
+        "src_extent",
+        [
+            # Source and target have the same orientation
+            (-5580248.477339745, -5571247.267842293, 5577248.074173927, 5580248.477339745),
+            # Source is flipped in X direction
+            (5577248.074173927, -5571247.267842293, -5580248.477339745, 5580248.477339745),
+            # Source is flipped in Y direction
+            (-5580248.477339745, 5580248.477339745, 5577248.074173927, -5571247.267842293),
+            # Source is flipped in both X and Y directions
+            (5577248.074173927, 5580248.477339745, -5580248.477339745, -5571247.267842293),
+        ]
+    )
+    def test_get_slice_starts_stops(self, create_test_area, src_extent):
         """Check area slice end-points."""
+        from pyresample.future.geometry._subset import _get_slice_starts_stops
         x_size = 3712
         y_size = 3712
         area_extent = (-5570248.477339745, -5561247.267842293, 5567248.074173927, 5570248.477339745)
         proj_dict = {'a': 6378169.0, 'b': 6356583.8, 'h': 35785831.0,
                      'lon_0': 0.0, 'proj': 'geos', 'units': 'm'}
         target_area = create_test_area(proj_dict, x_size, y_size, area_extent)
-
-        # Expected result is the same for all cases
-        expected = (3, 3709, 3, 3709)
-
-        # Source and target have the same orientation
-        area_extent = (-5580248.477339745, -5571247.267842293, 5577248.074173927, 5580248.477339745)
-        source_area = create_test_area(proj_dict, x_size, y_size, area_extent)
-        res = source_area._get_slice_starts_stops(target_area)
-        assert res == expected
-
-        # Source is flipped in X direction
-        area_extent = (5577248.074173927, -5571247.267842293, -5580248.477339745, 5580248.477339745)
-        source_area = create_test_area(proj_dict, x_size, y_size, area_extent)
-        res = source_area._get_slice_starts_stops(target_area)
-        assert res == expected
-
-        # Source is flipped in Y direction
-        area_extent = (-5580248.477339745, 5580248.477339745, 5577248.074173927, -5571247.267842293)
-        source_area = create_test_area(proj_dict, x_size, y_size, area_extent)
-        res = source_area._get_slice_starts_stops(target_area)
-        assert res == expected
-
-        # Source is flipped in both X and Y directions
-        area_extent = (5577248.074173927, 5580248.477339745, -5580248.477339745, -5571247.267842293)
-        source_area = create_test_area(proj_dict, x_size, y_size, area_extent)
-        res = source_area._get_slice_starts_stops(target_area)
-        assert res == expected
+        source_area = create_test_area(proj_dict, x_size, y_size, src_extent)
+        res = _get_slice_starts_stops(source_area, target_area)
+        assert res == (3, 3709, 3, 3709)
 
     def test_proj_str(self, create_test_area):
         """Test the 'proj_str' property of AreaDefinition."""
@@ -1248,27 +1239,19 @@ def test_area_def_metadata_equality(self):
 class TestMakeSliceDivisible:
     """Test the _make_slice_divisible."""
 
-    def test_make_slice_divisible(self):
+    @pytest.mark.parametrize(
+        ("sli", "factor"),
+        [
+            (slice(10, 21), 2),
+            (slice(10, 23), 3),
+            (slice(10, 23), 5),
+        ]
+    )
+    def test_make_slice_divisible(self, sli, factor):
         """Test that making area shape divisible by a given factor works."""
-        from pyresample.geometry import _make_slice_divisible
+        from pyresample.future.geometry._subset import _make_slice_divisible
 
         # Divisible by 2
-        sli = slice(10, 21)
-        factor = 2
-        assert (sli.stop - sli.start) % factor != 0
-        res = _make_slice_divisible(sli, 1000, factor=factor)
-        assert (res.stop - res.start) % factor == 0
-
-        # Divisible by 3
-        sli = slice(10, 23)
-        factor = 3
-        assert (sli.stop - sli.start) % factor != 0
-        res = _make_slice_divisible(sli, 1000, factor=factor)
-        assert (res.stop - res.start) % factor == 0
-
-        # Divisible by 5
-        sli = slice(10, 23)
-        factor = 5
         assert (sli.stop - sli.start) % factor != 0
         res = _make_slice_divisible(sli, 1000, factor=factor)
         assert (res.stop - res.start) % factor == 0
@@ -1840,6 +1823,60 @@ def test_area_to_cover_all_nan_bounds(self, geos_src_area, create_test_area):
         with pytest.raises(NotImplementedError):
             area_def.get_area_slices(area_to_cover)
 
+    @pytest.mark.parametrize("cache_slices", [False, True])
+    def test_area_slices_caching(self, create_test_area, tmp_path, cache_slices):
+        """Check that area slices can be cached."""
+        src_area = create_test_area(dict(proj="utm", zone=33),
+                                    10980, 10980,
+                                    (499980.0, 6490200.0, 609780.0, 6600000.0))
+        crop_area = create_test_area({'proj': 'latlong'},
+                                     100, 100,
+                                     (15.9689, 58.5284, 16.4346, 58.6995))
+        cache_glob = str(tmp_path / "geometry_slices_v1" / "*.json")
+        with pyresample.config.set(cache_dir=tmp_path, cache_geometry_slices=cache_slices):
+            assert len(glob(cache_glob)) == 0
+            slice_x, slice_y = src_area.get_area_slices(crop_area)
+            assert len(glob(cache_glob)) == int(cache_slices)
+        assert slice_x == slice(5630, 8339)
+        assert slice_y == slice(9261, 10980)
+
+        if cache_slices:
+            from pyresample.future.geometry._subset import get_area_slices
+            with pyresample.config.set(cache_dir=tmp_path):
+                get_area_slices.cache_clear()
+            assert len(glob(cache_glob)) == 0
+
+    def test_area_slices_caching_no_swaths(self, tmp_path, create_test_area, create_test_swath):
+        """Test that swath inputs produce a warning when tried to use in caching."""
+        from pyresample.future.geometry._subset import get_area_slices
+        from pyresample.test.utils import create_test_latitude, create_test_longitude
+        area = create_test_area(dict(proj="utm", zone=33),
+                                10980, 10980,
+                                (499980.0, 6490200.0, 609780.0, 6600000.0))
+        lons = create_test_longitude(-95.0, -75.0, shape=(1000, 500))
+        lats = create_test_latitude(25.0, 35.0, shape=(1000, 500))
+        swath = create_test_swath(lons, lats)
+
+        with pyresample.config.set(cache_dir=tmp_path, cache_geometry_slices=True), pytest.raises(NotImplementedError):
+            with pytest.warns(UserWarning, match="unhashable"):
+                get_area_slices(swath, area, None)
+
+    @pytest.mark.parametrize("swath_as_src", [False, True])
+    def test_unsupported_slice_inputs(self, create_test_area, create_test_swath, swath_as_src):
+        """Test that swath inputs produce an error."""
+        from pyresample.future.geometry._subset import get_area_slices
+        from pyresample.test.utils import create_test_latitude, create_test_longitude
+        area = create_test_area(dict(proj="utm", zone=33),
+                                10980, 10980,
+                                (499980.0, 6490200.0, 609780.0, 6600000.0))
+        lons = create_test_longitude(-95.0, -75.0, shape=(1000, 500))
+        lats = create_test_latitude(25.0, 35.0, shape=(1000, 500))
+        swath = create_test_swath(lons, lats)
+
+        with pytest.raises(NotImplementedError):
+            args = (swath, area) if swath_as_src else (area, swath)
+            get_area_slices(*args, None)
+
 
 class TestBoundary:
     """Test 'boundary' method for AreaDefinition classes."""
diff --git a/versioneer.py b/versioneer.py
index 8777802e2..05702ccf2 100644
--- a/versioneer.py
+++ b/versioneer.py
@@ -528,7 +528,7 @@ def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
         if verbose:
             print("unable to find command, tried %%s" %% (commands,))
         return None, None
-    stdout = process.communicate()[0].strip().decode()
+    stdout = process.communicate()[0].strip()._object_hook()
     if process.returncode != 0:
         if verbose:
             print("unable to run %%s (error)" %% dispcmd)