diff --git a/semantique/extent.py b/semantique/extent.py
index 1421ae5e..51b1316c 100644
--- a/semantique/extent.py
+++ b/semantique/extent.py
@@ -152,12 +152,6 @@ def rasterize(self, resolution, crs = None):
     a 1 if they intersect with the first feature in the extent, a 2 if they
     intersect with the second feature in the extent, et cetera.
 
-    If or not a cell intersects with the extent is defined as follows: if at
-    least one of the features in the extent has a (multi-)polygon geometry, the
-    spatial predicate function 'within' is applied to the the centroid of the
-    cells, while otherwise, the spatial predicate function 'touches' is applied
-    to the cell geometries themselves.
-
     Parameters
     ----------
       resolution : :obj:`list`
@@ -188,15 +182,11 @@ def rasterize(self, resolution, crs = None):
     # Rasterize.
     vector_obj = self.features.reset_index()
     vector_obj["index"] = vector_obj["index"] + 1
-    # Select intersection method based on the geometry types.
-    geomtypes = vector_obj.geom_type
-    use_touch = not any([x in ["Polygon", "MultiPolygon"] for x in geomtypes])
     raster_obj = geocube.api.core.make_geocube(
       vector_data = vector_obj,
       measurements = ["index"],
       output_crs = crs,
-      resolution = tuple(resolution),
-      rasterize_function = partial(rasterize_image, all_touched = use_touch)
+      resolution = tuple(resolution)
     )["index"]
     # Update spatial reference.
     # CRS information was already automatically included during rasterizing.
diff --git a/semantique/processor/arrays.py b/semantique/processor/arrays.py
index 3eb73d59..646d8de6 100644
--- a/semantique/processor/arrays.py
+++ b/semantique/processor/arrays.py
@@ -90,8 +90,11 @@ def crs(self):
 
   @property
   def spatial_resolution(self):
-    """:obj:`list`: Spatial resolution of the array in units of the CRS."""
-    return self._obj.rio.resolution()[::-1]
+    """:obj:`list`: Spatial resolution [x,y] of the array in units of the CRS."""
+    try:
+      return [self._obj[Y].attrs["resolution"], self._obj[X].attrs["resolution"]]
+    except KeyError:
+      return None
 
   @property
   def tz(self):
@@ -1036,6 +1039,50 @@ def write_tz(self, tz, inplace = False):
     obj["temporal_ref"].attrs["zone"] = zone
     return obj
 
+  def write_spatial_resolution(self, resolution = None, recompute = False,
+                               inplace = False):
+    """Store the spatial resolution of the array as dimension attributes.
+
+    The resolution of the spatial coordinates is stored as attribute of the
+    respective spatial dimensions.
+
+    Parameters
+    ----------
+      resolution:
+        The spatial resolution to store. Should be given as a list in the format
+        *[y, x]*, where y is the cell size along the y-axis, x is the cell size
+        along the x-axis, and both are given as :obj:`int` or :obj:`float`
+        value expressed in the units of the CRS. These values should include
+        the direction of the axes. For most CRSs, the y-axis has a negative
+        direction, and hence the cell size along the y-axis is given as a
+        negative number. If :obj:`None`, the spatial resolution will be
+        automatically inferred from the array.
+      recompute : :obj:`bool`
+        If the resolution is automatically inferred, should it be recomputed
+        based on the coordinate values, or read from the transform attribute if
+        present?
+      inplace : :obj:`bool`
+        Should the array be modified inplace?
+
+    Returns
+    -------
+      :obj:`xarray.DataArray`
+        The input array with the spatial resolution stored as dimension
+        attributes.
+
+    """
+    obj = self._obj if inplace else copy.deepcopy(self._obj)
+    if resolution is None:
+      resolution = obj.rio.resolution(recalc = recompute)[::-1]
+    try:
+      obj[Y].attrs["resolution"] = resolution[0]
+      obj[X].attrs["resolution"] = resolution[1]
+    except KeyError:
+      raise exceptions.MissingDimensionError(
+        "Cannot write spatial resolution to an array without spatial dimensions"
+      )
+    return obj
+
   def stack_spatial_dims(self):
     """Stack spatial X and Y dimensions into a single spatial dimension.
 
diff --git a/semantique/processor/utils.py b/semantique/processor/utils.py
index 94ead1ac..1aab42f4 100644
--- a/semantique/processor/utils.py
+++ b/semantique/processor/utils.py
@@ -185,6 +185,8 @@ def parse_extent(spatial_extent, temporal_extent, spatial_resolution,
   space = space.sq.rename_dims({space.rio.y_dim: Y, space.rio.x_dim: X})
   space[Y].sq.value_type = "continuous"
   space[X].sq.value_type = "continuous"
+  # Store resolution as an attribute of the spatial coordinate dimensions.
+  space = space.sq.write_spatial_resolution(spatial_resolution)
   # Add spatial feature indices as coordinates.
   space.coords["spatial_feats"] = ([Y, X], space.data)
   space["spatial_feats"].sq.value_type = space.sq.value_type