diff --git a/docs/_static/style.css b/docs/_static/style.css
index d749760..2d5f2ce 100644
--- a/docs/_static/style.css
+++ b/docs/_static/style.css
@@ -1,5 +1,5 @@
-.sidebar-log {
-  max-width: 70%;
+.sidebar-logo {
+  max-width: 50%;
 }
 
 .xr-wrap {
diff --git a/docs/heuristics.md b/docs/heuristics.md
deleted file mode 100644
index 567a261..0000000
--- a/docs/heuristics.md
+++ /dev/null
@@ -1,148 +0,0 @@
-# `assign_index` Heuristics
-
-This page documents the heuristics and decision-making logic used by rasterix when working with spatial data.
-
-## Affine Transform Detection in `assign_index`
-
-When calling {py:func}`assign_index()` to create a `RasterIndex` for an Xarray object, rasterix needs to determine the affine transformation that maps pixel coordinates to spatial coordinates. The function follows a specific priority order when searching for this information.
-
-### Priority Order
-
-Rasterix checks for transform information in the following order:
-
-1. `GeoTransform` attribute on a CF 'grid mapping' variable (commonly `spatial_ref`).
-1. STAC `proj:transform` attribute
-1. GeoTIFF metadata: (`model_tiepoint` and `model_pixel_scale` for now)
-
-Each method is tried in sequence, and the first successful match is used.
-
-If these fail, an index is constructed from explicit coordinate arrays if present.
-
-### 1. `GeoTransform`
-
-```{seealso}
-[GDAL GeoTransform tutorial](https://gdal.org/en/stable/tutorials/geotransforms_tut.html)
-```
-
-The function first looks for a CF conventions "grid mapping" variable (commonly named `spatial_ref`). The grid mapping variable is identified by:
-
-- For DataArrays: the `grid_mapping` attribute on the DataArray itself
-- For Datasets: the `grid_mapping` attribute on the first data variable
-- As a fallback: a coordinate variable named `spatial_ref`
-
-If found, rasterix checks for a `GeoTransform` attribute on this variable, which should be a GDAL-format geotransform string with 6 space-separated numbers.
-
-**Format**: `"c a b f d e"` (GDAL format)
-
-**Example**:
-
-```python
-da.coords["spatial_ref"].attrs["GeoTransform"] = "323400.0 30.0 0.0 4265400.0 0.0 30.0"
-```
-
-### 2. STAC `proj:transform`
-
-```{seealso}
-[STAC Projection Extension](https://github.com/stac-extensions/projection)
-```
-
-If no GeoTransform is found, rasterix checks the DataArray's attributes for a STAC projection extension `proj:transform` field. This represents the affine transformation as a flat array.
-
-**Format**: `[a, b, c, d, e, f]` or `[a, b, c, d, e, f, 0, 0, 1]`
-
-The transformation follows this matrix formula:
-
-```
-[Xp]   [a  b  c]   [Pixel]
-[Yp] = [d  e  f] * [Line ]
-[1 ]   [0  0  1]   [1    ]
-```
-
-Where:
-
-- `a` = pixel width (x-scale)
-- `b` = row rotation (typically 0)
-- `c` = x-coordinate of upper-left pixel corner
-- `d` = column rotation (typically 0)
-- `e` = pixel height (y-scale, negative if y decreases)
-- `f` = y-coordinate of upper-left pixel corner
-
-**Example**:
-
-```python
-da.attrs["proj:transform"] = [30.0, 0.0, 323400.0, 0.0, 30.0, 4268400.0]
-```
-
-### 3. GeoTIFF Metadata
-
-```{seealso}
-[GeoTIFF specification on coordinate transformations](https://docs.ogc.org/is/19-008r4/19-008r4.html#_geotiff_tags_for_coordinate_transformations)
-```
-
-If no STAC transform is found, rasterix checks for GeoTIFF-style metadata using model tiepoints and pixel scale.
-
-**Format**:
-
-- `model_tiepoint`: `[I, J, K, X, Y, Z]` - Maps pixel `(I, J, K)` to world coordinates `(X, Y, Z)`
-- `model_pixel_scale`: `[ScaleX, ScaleY, ScaleZ]` - Pixel dimensions in world coordinates
-
-**Constraints**:
-
-- `ScaleZ` must be 0 (only 2D rasters are supported)
-- If the tiepoint is at pixel `(I, J)`, the affine is computed as:
-  - `c = X - I * ScaleX`
-  - `f = Y - J * ScaleY`
-
-**Example**:
-
-```python
-da.attrs["model_tiepoint"] = [0.0, 0.0, 0.0, 323400.0, 4265400.0, 0.0]
-da.attrs["model_pixel_scale"] = [30.0, 30.0, 0.0]
-```
-
-**Trace log**: `"Creating affine from GeoTIFF model_tiepoint and model_pixel_scale attributes"`
-
-### 4. Coordinate Arrays
-
-If no metadata is found, rasterix falls back to computing the affine transformation from 1D coordinate arrays.
-
-**Requirements**:
-
-- Coordinate variables must exist for both x and y dimensions
-- Coordinates must be 1D arrays
-- Coordinates must have at least 2 values to compute spacing
-
-**Computation**:
-
-- Pixel spacing: `dx = x[1] - x[0]`, `dy = y[1] - y[0]`
-- Origin calculation accounts for whether y increases or decreases
-- Assumes pixel-centered coordinates, adjusts to pixel corners
-
-**Example**:
-
-```python
-da = xr.DataArray(
-    data,
-    coords={
-        "x": np.arange(0.5, 100.5),  # Pixel-centered
-        "y": np.arange(0.5, 100.5),
-    },
-)
-```
-
-### Logging
-
-To see which method is being used, enable trace-level logging:
-
-```python
-import logging
-
-logging.basicConfig(level=logging.TRACE)
-logging.getLogger("rasterix").setLevel(logging.TRACE)
-```
-
-This will output messages like:
-
-```
-TRACE:rasterix:Creating affine from STAC proj:transform attribute
-```
diff --git a/docs/index.md b/docs/index.md
index 24dbe07..28342ba 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -8,12 +8,14 @@ caption: Raster Index
 hidden:
 ---
 raster_index/intro
+raster_index/creating
 raster_index/indexing
 raster_index/crs
 raster_index/aligning
 raster_index/combining
 raster_index/tolerance
 raster_index/design_choices
+raster_index/heuristics
 ```
 
 ```{toctree}
diff --git a/docs/raster_index/aligning.md b/docs/raster_index/aligning.md
index bdaa8e1..c898a19 100644
--- a/docs/raster_index/aligning.md
+++ b/docs/raster_index/aligning.md
@@ -8,7 +8,7 @@ kernelspec:
 ---
 
 ```{eval-rst}
-.. currentmodule:: rasterix
+.. currentmodule:: rasterix.raster_index
 ```
 
 ```{code-cell} python
diff --git a/docs/raster_index/combining.md b/docs/raster_index/combining.md
index 45c1eed..0e535c7 100644
--- a/docs/raster_index/combining.md
+++ b/docs/raster_index/combining.md
@@ -8,7 +8,7 @@ kernelspec:
 ---
 
 ```{eval-rst}
-.. currentmodule:: rasterix
+.. currentmodule:: rasterix.raster_index
 ```
 
 ```{code-cell} python
diff --git a/docs/raster_index/creating.md b/docs/raster_index/creating.md
new file mode 100644
index 0000000..ce0b5fe
--- /dev/null
+++ b/docs/raster_index/creating.md
@@ -0,0 +1,185 @@
+---
+jupytext:
+  text_representation:
+    format_name: myst
+kernelspec:
+  display_name: Python 3
+  name: python
+---
+
+```{eval-rst}
+.. currentmodule:: rasterix.raster_index
+```
+
+```{code-cell} python
+---
+tags: [remove-cell]
+---
+%xmode minimal
+import xarray as xr
+xr.set_options(display_expand_indexes=True);
+```
+
+# Creating
+
+There are several ways to create a {py:class}`RasterIndex`, depending on the source of your data and the metadata conventions it uses.
+
+```{code-cell}
+import numpy as np
+import xarray as xr
+import rasterix
+```
+
+## Using `assign_index`
+
+The easiest way to create a RasterIndex is using {py:func}`assign_index`. This function automatically detects the affine transform from various metadata conventions, including:
+
+- [GDAL GeoTransform](https://gdal.org/en/stable/tutorials/geotransforms_tut.html) (from rioxarray/GeoTIFF files)
+- [GeoTIFF](https://docs.ogc.org/is/19-008r4/19-008r4.html) tiepoint and pixel scale attributes
+- [STAC projection extension](https://github.com/stac-extensions/projection) `proj:transform` metadata
+- [Zarr Spatial Convention](https://zarr-specs.readthedocs.io/en/latest/v3/conventions/spatial/v1.0.html) (`spatial:transform`)
+- 1D coordinate arrays (common in NetCDF files)
+
+See {doc}`/raster_index/heuristics` for the full priority order and detection logic.
+
+```{code-cell}
+import pyproj
+
+# Example: dataset with GDAL-style GeoTransform metadata
+ds = xr.Dataset(
+    {"temperature": (("y", "x"), np.random.rand(100, 100))},
+    coords={
+        "spatial_ref": (
+            (),
+            0,
+            pyproj.CRS.from_epsg(32610).to_cf() | {"GeoTransform": "400000.0 10.0 0.0 5000000.0 0.0 -10.0"},
+        )
+    },
+    attrs={"grid_mapping": "spatial_ref"},
+)
+
+# assign_index auto-detects the convention and creates the RasterIndex
+ds = rasterix.assign_index(ds)
+ds
+```
+
+## Direct construction with class methods
+
+For more control, you can create a {py:class}`RasterIndex` directly using class methods. This is useful when you have the transform parameters available directly, or when working with data that doesn't have embedded metadata.
+
+### From an Affine transform
+
+Use {py:meth}`RasterIndex.from_transform` to create from an [Affine](https://github.com/rasterio/affine) transform directly (corresponds to [GDAL GeoTransform](https://gdal.org/en/stable/tutorials/geotransforms_tut.html) convention):
+
+```{code-cell}
+from affine import Affine
+
+transform = Affine(10.0, 0.0, 400000.0, 0.0, -10.0, 5000000.0)
+index = rasterix.RasterIndex.from_transform(
+    transform,
+    width=100,
+    height=100,
+    crs="EPSG:32610",
+)
+index
+```
+
+```{code-cell}
+# Assign to data
+ds_manual = xr.Dataset(
+    {"data": (("y", "x"), np.random.rand(100, 100))},
+    coords=xr.Coordinates.from_xindex(index),
+)
+ds_manual
+```
+
+### From a GeoTransform
+
+Use {py:meth}`RasterIndex.from_geotransform` to create from a [GDAL GeoTransform](https://gdal.org/en/stable/tutorials/geotransforms_tut.html), either as a sequence or a space-separated string:
+
+```{code-cell}
+# From a tuple
+geotransform = (400000.0, 10.0, 0.0, 5000000.0, 0.0, -10.0)
+index = rasterix.RasterIndex.from_geotransform(
+    geotransform,
+    width=100,
+    height=100,
+    crs="EPSG:32610",
+)
+index
+```
+
+```{code-cell}
+# From a string (as commonly stored in netCDF attributes)
+geotransform = "400000.0 10.0 0.0 5000000.0 0.0 -10.0"
+index = rasterix.RasterIndex.from_geotransform(
+    geotransform,
+    width=100,
+    height=100,
+    crs="EPSG:32610",
+)
+index
+```
+
+### From GeoTIFF tiepoint and scale
+
+Use {py:meth}`RasterIndex.from_tiepoint_and_scale` to create from [GeoTIFF](https://docs.ogc.org/is/19-008r4/19-008r4.html)-style tiepoint and pixel scale attributes:
+
+```{code-cell}
+index = rasterix.RasterIndex.from_tiepoint_and_scale(
+    tiepoint=[0.0, 0.0, 0.0, 323400.0, 4265400.0, 0.0],
+    scale=[30.0, 30.0, 0.0],
+    width=100,
+    height=100,
+    crs="EPSG:32610",
+)
+index
+```
+
+### From STAC projection metadata
+
+Use {py:meth}`RasterIndex.from_stac_proj_metadata` to create from [STAC projection extension](https://github.com/stac-extensions/projection) metadata:
+
+```{code-cell}
+metadata = {"proj:transform": [30.0, 0.0, 323400.0, 0.0, -30.0, 4268400.0]}
+index = rasterix.RasterIndex.from_stac_proj_metadata(
+    metadata,
+    width=100,
+    height=100,
+    crs="EPSG:32610",
+)
+index
+```
+
+## Accessing transform information
+
+Once created, you can access the affine transform using methods on {py:class}`RasterIndex`:
+
+```{code-cell}
+ds = rasterix.assign_index(
+    xr.Dataset(
+        {"data": (("y", "x"), np.random.rand(100, 100))},
+        coords={
+            "spatial_ref": ((), 0, {"GeoTransform": "400000.0 10.0 0.0 5000000.0 0.0 -10.0"}),
+        },
+    )
+)
+
+# Top-left corner transform (GDAL convention) via transform()
+ds.xindexes["x"].transform()
+```
+
+```{code-cell}
+# Pixel center transform via center_transform()
+ds.xindexes["x"].center_transform()
+```
+
+```{code-cell}
+# Bounding box via the bbox property
+ds.xindexes["x"].bbox
+```
+
+```{code-cell}
+# As GeoTransform string (for saving) via as_geotransform()
+ds.xindexes["x"].as_geotransform()
+```
diff --git a/docs/raster_index/design_choices.md b/docs/raster_index/design_choices.md
index 0a63e63..2973cda 100644
--- a/docs/raster_index/design_choices.md
+++ b/docs/raster_index/design_choices.md
@@ -1,5 +1,5 @@
 ```{eval-rst}
-.. currentmodule:: rasterix
+.. currentmodule:: rasterix.raster_index
 ```
 
 # Design Choices
@@ -55,7 +55,7 @@ GDAL _chooses_ to track the affine transform using a `GeoTransform` attribute on
 the grid mapping variable. Thus, {py:class}`RasterIndex` _cannot_ keep the `"GeoTransform"` attribute on `"spatial_ref"` up-to-date
 because it does not control it.
 
-Thus, {py:func}`assign_index` will delete the `"GeoTransform"` attribute on the grid mapping variable if it is detected, after using it
+Thus, {py:func}`~rasterix.assign_index` will delete the `"GeoTransform"` attribute on the grid mapping variable if it is detected, after using it
 to construct the affine matrix.
 
 If you wish to extract the GeoTransform attribute to write it to a location of your choosing use {py:meth}`RasterIndex.as_geotransform`.
diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 89c7591..cc637fa 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -595,6 +595,89 @@ def from_transform(
 
         return cls(index, crs=crs)
 
+    @classmethod
+    def from_geotransform(
+        cls,
+        geotransform: Sequence[float] | str,
+        *,
+        width: int,
+        height: int,
+        x_dim: str = "x",
+        y_dim: str = "y",
+        x_coord_name: str = "xc",
+        y_coord_name: str = "yc",
+        crs: CRS | Any | None = None,
+    ) -> RasterIndex:
+        """Create a RasterIndex from a GDAL-style GeoTransform.
+
+        Parameters
+        ----------
+        geotransform : sequence of float or str
+            GDAL GeoTransform as a 6-element sequence (c, a, b, f, d, e) or a
+            space-separated string of 6 numbers. The elements are:
+            - c: x-coordinate of the upper-left corner of the upper-left pixel
+            - a: pixel width (x-direction resolution)
+            - b: row rotation (typically 0)
+            - f: y-coordinate of the upper-left corner of the upper-left pixel
+            - d: column rotation (typically 0)
+            - e: pixel height (y-direction resolution, typically negative)
+        width : int
+            Number of pixels in the x direction.
+        height : int
+            Number of pixels in the y direction.
+        x_dim : str, optional
+            Name for the x dimension.
+        y_dim : str, optional
+            Name for the y dimension.
+        x_coord_name : str, optional
+            Name for the x coordinate. For non-rectilinear transforms only.
+        y_coord_name : str, optional
+            Name for the y coordinate. For non-rectilinear transforms only.
+        crs : :class:`pyproj.crs.CRS` or any, optional
+            The coordinate reference system. Any value accepted by
+            :meth:`pyproj.crs.CRS.from_user_input`.
+
+        Returns
+        -------
+        RasterIndex
+            A new RasterIndex object with appropriate internal structure.
+
+        See Also
+        --------
+        from_transform : Create from an Affine transform.
+        as_geotransform : Convert RasterIndex back to GeoTransform string.
+
+        References
+        ----------
+        - `GDAL GeoTransform tutorial <https://gdal.org/en/stable/tutorials/geotransforms_tut.html>`_
+
+        Examples
+        --------
+        Create from a sequence:
+
+        >>> geotransform = (323400.0, 30.0, 0.0, 4265400.0, 0.0, -30.0)
+        >>> index = RasterIndex.from_geotransform(geotransform, width=100, height=100)
+
+        Create from a string (as stored in netCDF attributes):
+
+        >>> geotransform = "323400.0 30.0 0.0 4265400.0 0.0 -30.0"
+        >>> index = RasterIndex.from_geotransform(geotransform, width=100, height=100)
+        """
+        if isinstance(geotransform, str):
+            geotransform = tuple(map(float, geotransform.split()))
+
+        affine = Affine.from_gdal(*geotransform)
+        return cls.from_transform(
+            affine,
+            width=width,
+            height=height,
+            x_dim=x_dim,
+            y_dim=y_dim,
+            x_coord_name=x_coord_name,
+            y_coord_name=y_coord_name,
+            crs=crs,
+        )
+
     @classmethod
     def from_tiepoint_and_scale(
         cls,