diff --git a/xarray/core/indexing.py b/xarray/core/indexing.py
index 9a29b63f4d0..07cfc86d687 100644
--- a/xarray/core/indexing.py
+++ b/xarray/core/indexing.py
@@ -22,6 +22,7 @@
 from .types import T_Xarray
 from .utils import (
     NDArrayMixin,
+    cached_property,
     either_dict_or_kwargs,
     get_valid_numpy_dtype,
     safe_cast_to_index,
@@ -472,7 +473,7 @@ def __getitem__(self, key):
 class LazilyIndexedArray(ExplicitlyIndexedNDArrayMixin):
     """Wrap an array to make basic and outer indexing lazy."""
 
-    __slots__ = ("array", "key")
+    __slots__ = ("array", "key", "__dict__")
 
     def __init__(self, array, key=None):
         """
@@ -509,7 +510,7 @@ def _updated_key(self, new_key):
             return BasicIndexer(full_key)
         return OuterIndexer(full_key)
 
-    @property
+    @cached_property
     def shape(self) -> tuple[int, ...]:
         shape = []
         for size, k in zip(self.array.shape, self.key.tuple):
@@ -552,7 +553,7 @@ def __repr__(self):
 class LazilyVectorizedIndexedArray(ExplicitlyIndexedNDArrayMixin):
     """Wrap an array to make vectorized indexing lazy."""
 
-    __slots__ = ("array", "key")
+    __slots__ = ("array", "key", "__dict__")
 
     def __init__(self, array, key):
         """
@@ -568,7 +569,7 @@ def __init__(self, array, key):
             self.key = _arrayize_vectorized_indexer(key, array.shape)
         self.array = as_indexable(array)
 
-    @property
+    @cached_property
     def shape(self) -> tuple[int, ...]:
         return np.broadcast(*self.key.tuple).shape
 
@@ -1367,7 +1368,7 @@ def transpose(self, order):
 class PandasIndexingAdapter(ExplicitlyIndexedNDArrayMixin):
     """Wrap a pandas.Index to preserve dtypes and handle explicit indexing."""
 
-    __slots__ = ("array", "_dtype")
+    __slots__ = ("array", "_dtype", "__dict__")
 
     def __init__(self, array: pd.Index, dtype: DTypeLike = None):
         self.array = safe_cast_to_index(array)
@@ -1391,7 +1392,7 @@ def __array__(self, dtype: DTypeLike = None) -> np.ndarray:
                 array = array.astype("object")
         return np.asarray(array.values, dtype=dtype)
 
-    @property
+    @cached_property
     def shape(self) -> tuple[int, ...]:
         return (len(self.array),)
 
diff --git a/xarray/core/utils.py b/xarray/core/utils.py
index ab3f8d3a282..67d463514f8 100644
--- a/xarray/core/utils.py
+++ b/xarray/core/utils.py
@@ -977,3 +977,11 @@ def contains_only_dask_or_numpy(obj) -> bool:
             for var in obj.variables.values()
         ]
     )
+
+
+class cached_property(functools.cached_property):
+    """Read only version of functools.cached_property."""
+
+    def __set__(self, instance, val):
+        """Raise an error when attempting to set a cached property."""
+        raise AttributeError("Can't set attribute")
diff --git a/xarray/core/variable.py b/xarray/core/variable.py
index 90edf652284..bbc3eacf829 100644
--- a/xarray/core/variable.py
+++ b/xarray/core/variable.py
@@ -48,6 +48,7 @@
     NdimSizeLenMixin,
     OrderedSet,
     _default,
+    cached_property,
     decode_numpy_dict_values,
     drop_dims_from_indexers,
     either_dict_or_kwargs,
@@ -292,7 +293,7 @@ class Variable(AbstractArray, NdimSizeLenMixin, VariableArithmetic):
     they can use more complete metadata in context of coordinate labels.
     """
 
-    __slots__ = ("_dims", "_data", "_attrs", "_encoding")
+    __slots__ = ("_dims", "_data", "_attrs", "_encoding", "__dict__")
 
     def __init__(self, dims, data, attrs=None, encoding=None, fastpath=False):
         """
@@ -327,7 +328,7 @@ def __init__(self, dims, data, attrs=None, encoding=None, fastpath=False):
     def dtype(self):
         return self._data.dtype
 
-    @property
+    @cached_property
     def shape(self):
         return self._data.shape