diff --git a/Cargo.lock b/Cargo.lock
index 1d53519..8b58e75 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -3375,7 +3375,7 @@ checksum = "ea2f10b9bb0928dfb1b42b65e1f9e36f7f54dbdf08457afefb38afcdec4fa2bb"
 
 [[package]]
 name = "xarray_sql"
-version = "0.2.2"
+version = "0.2.3"
 dependencies = [
  "arrow",
  "async-stream",
diff --git a/README.md b/README.md
index ee512dc..ddbd497 100644
--- a/README.md
+++ b/README.md
@@ -57,52 +57,37 @@ ctx.sql('''
 # 0  8.640069
 
 # Average temperature per pressure level, globally. 
-ctx.sql('''
+result = ctx.sql('''
   SELECT level, AVG(temperature) - 273.15 AS avg_c
   FROM era5.atmosphere
   WHERE time BETWEEN TIMESTAMP '2020-01-01'
                  AND TIMESTAMP '2020-01-01 05:00:00'
   GROUP BY level
   ORDER BY level DESC
-''').to_pandas()
-#     level      avg_c
-# 0    1000   6.621012   ← surface
-# 1     975   5.185638
-# 2     950   4.028429
-# 3     925   3.082812
-# 4     900   2.210917
-# 5     875   1.395018
-# 6     850   0.634267
-# 7     825  -0.210372
-# 8     800  -1.181075
-# 9     775  -2.306465
-# 10    750  -3.535534
-# 11    700  -6.241685
-# 12    650  -9.236364
-# 13    600 -12.580938
-# 14    550 -16.335386
-# 15    500 -20.643604
-# 16    450 -25.573401
-# 17    400 -31.156920
-# 18    350 -37.400552
-# 19    300 -43.852607
-# 20    250 -49.322132
-# 21    225 -51.569113
-# 22    200 -53.693248
-# 23    175 -55.890484
-# 24    150 -58.382290
-# 25    125 -61.091916
-# 26    100 -63.624885   ← tropopause
-# 27     70 -63.182300
-# 28     50 -60.124845
-# 29     30 -55.986327
-# 30     20 -52.433089
-# 31     10 -44.140750
-# 32      7 -38.707350
-# 33      5 -32.621999
-# 34      3 -21.509175
-# 35      2 -13.355764
-# 36      1  -9.020513   ← top of atmosphere
+''')
+# DataFrame()
+# +-------+----------------------+
+# | level | avg_c                |
+# +-------+----------------------+
+# | 1000  | 6.6210120796502565   |
+# | 975   | 5.185637919348153    |
+# | 950   | 4.028428657263021    |
+# | 925   | 3.0828117974912743   |
+# | 900   | 2.2109172992531967   |
+# | 875   | 1.395017610194202    |
+# | 850   | 0.6342670572626616   |
+# | 825   | -0.21037158786759846 |
+# | 800   | -1.1810754318269687  |
+# | 775   | -2.3064649711534457  |
+# +-------+----------------------+
+
+avg_temp_ds = result.to_dataset(dims=["level"])
+# <xarray.Dataset> Size: 592B
+# Dimensions:  (level: 37)
+# Coordinates:
+#   * level    (level) int64 296B 1000 975 950 925 900 875 850 ... 20 10 7 5 3 2 1
+# Data variables:
+#     avg_c    (level) float64 296B 6.621 5.186 4.028 ... -21.51 -13.36 -9.021
 ```
 
 _(A runnable version of this example lives at
@@ -211,6 +196,8 @@ I want to give a special thanks to the following folks and institutions:
   who are working to make this library better.
 - Andrew Huang for the sense of taste he brings to the project and consummate code
   changes.
+- Aman Kumar for spending a considerable amount of his GSoC internship 
+  contributing to this project. 
 
 
 ## License
diff --git a/pyproject.toml b/pyproject.toml
index d9253bd..409da03 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -90,6 +90,7 @@ module = [
     "pyarrow.*",
     "datafusion.*",
     "xarray.*",
+    "pandas.*",
 ]
 ignore_missing_imports = true
 
diff --git a/tests/test_ds.py b/tests/test_ds.py
new file mode 100644
index 0000000..b674239
--- /dev/null
+++ b/tests/test_ds.py
@@ -0,0 +1,571 @@
+"""Tests for the SQL -> xarray reverse path.
+
+Covers the user-facing contract of ``ctx.sql(...).to_dataset(...)``:
+
+* Wrapper behavior on the object returned by ``ctx.sql`` and DataFusion
+  method passthrough.
+* Round-trip identity across varied source Datasets (one parametrized
+  ``assert_identical`` test, not eight per-aspect checks).
+* Aggregation, ``dimension_columns`` inference, and the template /
+  ``template`` resolution rules (name or Dataset) with their error paths.
+* Sparsity handling and ``fill_value`` dtype behavior.
+* The vectorized-indexer fallback through xarray's adapter.
+
+The tests favor checking the user-visible contract (values, dims,
+attrs) over the implementation path (call counts, internal class
+identity), so the suite stays useful as the lazy backend evolves.
+"""
+
+import numpy as np
+import pandas as pd
+import pytest
+import xarray as xr
+
+from xarray_sql import XarrayContext
+from xarray_sql.ds import XarrayDataFrame
+
+
+# ---------------------------------------------------------------------------
+# Wrapper: ctx.sql(...) returns XarrayDataFrame
+# ---------------------------------------------------------------------------
+
+
+def test_ctx_sql_returns_xarray_dataframe(air_dataset_small):
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    result = ctx.sql("SELECT * FROM air LIMIT 5")
+    assert isinstance(result, XarrayDataFrame)
+
+
+def test_to_pandas_unchanged_behavior(air_dataset_small):
+    """Wrapped ``.to_pandas()`` is bit-for-bit equal to the un-wrapped path."""
+    from datafusion import SessionContext
+
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    wrapped = ctx.sql("SELECT * FROM air LIMIT 7").to_pandas()
+    raw = SessionContext.sql(ctx, "SELECT * FROM air LIMIT 7").to_pandas()
+    pd.testing.assert_frame_equal(wrapped, raw)
+
+
+def test_passthrough_methods(air_dataset_small):
+    """DataFusion methods we did not override forward via ``__getattr__``."""
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    result = ctx.sql("SELECT * FROM air LIMIT 5")
+    names = [f.name for f in result.schema()]
+    assert {"lat", "lon", "time", "air"}.issubset(set(names))
+
+
+# ---------------------------------------------------------------------------
+# Round-trip identity (parametrized over local + tutorial datasets)
+# ---------------------------------------------------------------------------
+
+
+def _clear_encoding(ds: xr.Dataset) -> xr.Dataset:
+    """Strip ``encoding`` from a Dataset and all its variables.
+
+    Round-trip identity tests should not be coupled to encoding choices,
+    since template-recovery deliberately drops dtype-bound keys.
+    """
+    ds = ds.copy()
+    for v in ds.variables.values():
+        v.encoding.clear()
+    ds.encoding.clear()
+    return ds
+
+
+def _load_tutorial(name: str) -> xr.Dataset | None:
+    """Return a small xarray tutorial Dataset, or None when unavailable.
+
+    Used to widen round-trip coverage beyond the conftest fixtures without
+    requiring network in CI. Pooch caches downloads locally on first run.
+    """
+    try:
+        return xr.tutorial.open_dataset(name)
+    except (OSError, ValueError, ImportError):
+        return None
+
+
+@pytest.mark.parametrize(
+    "fixture_name",
+    ["air_dataset_small", "weather_dataset", "synthetic_dataset", "eraint_uvz"],
+)
+def test_round_trip_identity(request, fixture_name):
+    """``SELECT *`` round-trips to a Dataset that is ``assert_identical``
+    to the source: values, dims, coord values, dtypes, non-dim coords,
+    and attrs all match (modulo coord ordering, normalized on both
+    sides). One test covers what was previously a fan of narrow checks,
+    parametrized over local fixtures and one xarray tutorial dataset.
+    """
+    if fixture_name == "eraint_uvz":
+        source = _load_tutorial("eraint_uvz")
+        if source is None:
+            pytest.skip("eraint_uvz tutorial dataset unavailable")
+        source = source.chunk()
+    else:
+        source = request.getfixturevalue(fixture_name).copy()
+        source.attrs["round_trip_marker"] = "yes"
+        first_var = next(iter(source.data_vars))
+        source[first_var].attrs["units"] = "test_units"
+
+    ctx = XarrayContext()
+    ctx.from_dataset("t", source)
+    out = ctx.sql("SELECT * FROM t").to_dataset().compute()
+
+    sort_keys = list(out.dims)
+    actual = _clear_encoding(out.sortby(sort_keys))
+    expected = _clear_encoding(source.compute().sortby(sort_keys))
+    xr.testing.assert_identical(actual, expected)
+
+
+def test_aggregation_drops_dim(air_dataset_small):
+    """``GROUP BY lat, lon`` over time -> 2D Dataset with the alias."""
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    out = ctx.sql(
+        "SELECT lat, lon, AVG(air) AS air_avg FROM air GROUP BY lat, lon"
+    ).to_dataset(dims=["lat", "lon"])
+    assert set(out.dims) == {"lat", "lon"}
+    assert "air_avg" in out.data_vars
+    assert "air" not in out.data_vars
+    expected = (
+        air_dataset_small.compute()
+        .sortby(["lat", "lon"])
+        .mean(dim="time")["air"]
+        .values
+    )
+    actual = out.sortby(["lat", "lon"])["air_avg"].values
+    np.testing.assert_allclose(actual, expected)
+
+
+def test_barrier_query_scans_source_once(air_dataset_small):
+    """A barrier plan (aggregation) executes the source exactly once.
+
+    The lazy scan path re-runs the whole upstream plan for every coordinate
+    discovery and every variable access; for an aggregation -- which cannot push
+    an indexer filter below the GROUP BY -- that is pure re-computation of an
+    expensive scan. ``to_dataset()`` on a barrier plan must instead make a
+    single streamed pass over the source, and ``.compute()`` must trigger no
+    further reads.
+    """
+    from xarray_sql.df import block_slices
+    from xarray_sql.reader import read_xarray_table
+
+    reads: list = []
+    table = read_xarray_table(
+        air_dataset_small,
+        chunks={"time": 6},
+        _iteration_callback=lambda block, proj: reads.append(block),
+    )
+    n_partitions = len(list(block_slices(air_dataset_small, {"time": 6})))
+
+    ctx = XarrayContext()
+    ctx.register_table("air", table)
+    ctx._registered_datasets["air"] = air_dataset_small
+
+    out = ctx.sql(
+        "SELECT lat, lon, AVG(air) AS air_avg FROM air GROUP BY lat, lon"
+    ).to_dataset(dims=["lat", "lon"])
+    reads_after_construct = len(reads)
+    out.compute()
+    reads_after_compute = len(reads)
+
+    # Exactly one pass over the source (each partition read once) ...
+    assert reads_after_construct == n_partitions
+    # ... and computing the materialized result re-reads nothing.
+    assert reads_after_compute == reads_after_construct
+
+
+def test_order_by_direction_sets_dim_order(air_dataset_small):
+    """A barrier query's ORDER BY direction carries through to the Dataset
+    dimension order, rather than being force-sorted ascending.
+
+    ``ORDER BY lat DESC`` must yield a strictly descending ``lat`` dimension,
+    with data still correctly aligned to those (descending) coordinates.
+    """
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    out = ctx.sql(
+        "SELECT lat, AVG(air) AS air_avg FROM air GROUP BY lat ORDER BY lat DESC"
+    ).to_dataset(dims=["lat"])
+
+    lat = out["lat"].values
+    assert (np.diff(lat) < 0).all(), f"expected descending lat, got {lat}"
+
+    # Values stay aligned to the descending coordinate (scatter handles order).
+    expected = (
+        air_dataset_small.compute()
+        .mean(dim=["time", "lon"])["air"]
+        .sortby("lat", ascending=False)
+    )
+    np.testing.assert_allclose(out["air_avg"].values, expected.values)
+
+
+def test_chunks_argument_controls_partitioning(synthetic_dataset):
+    """``chunks`` controls eager-vs-chunked and inherits the source grid.
+
+    The default ``"inherit"`` reuses the source's genuinely multi-chunk
+    dimensions, so the output chunk grid maps onto the source partitions;
+    ``chunks=None`` forces an eager, in-memory result. Both reproduce the source.
+    """
+    import dask.array as da
+
+    ctx = XarrayContext()
+    ctx.from_dataset("t", synthetic_dataset)
+    var = next(iter(synthetic_dataset.data_vars))
+
+    inherited = ctx.sql("SELECT * FROM t").to_dataset()
+    assert isinstance(inherited[var].data, da.Array)
+    # Output time chunks align to the source's time partitions.
+    assert inherited.chunksizes["time"] == synthetic_dataset.chunksizes["time"]
+
+    eager = ctx.sql("SELECT * FROM t").to_dataset(chunks=None)
+    assert not isinstance(eager[var].data, da.Array)
+
+    xr.testing.assert_allclose(
+        inherited.compute().sortby(["time", "lat", "lon"]),
+        synthetic_dataset.compute().sortby(["time", "lat", "lon"]),
+    )
+
+
+def test_chunks_auto_snaps_to_source_partitions():
+    """``chunks="auto"`` coarsens to the byte budget but snaps chunk boundaries
+    to whole source partitions (so no chunk splits a source partition)."""
+    import dask
+
+    # 12 source partitions of size 2 along time.
+    ds = xr.Dataset(
+        {
+            "v": (
+                ("time", "x"),
+                np.arange(24 * 4, dtype="float64").reshape(24, 4),
+            )
+        },
+        coords={"time": np.arange(24), "x": np.arange(4)},
+    ).chunk({"time": 2})
+    ctx = XarrayContext()
+    ctx.from_dataset("t", ds)
+
+    # block bytes = 8 * 2(time) * 4(x) = 64; target 192 -> merge 3 partitions.
+    with dask.config.set({"array.chunk-size": "192B"}):
+        out = ctx.sql("SELECT * FROM t").to_dataset(chunks="auto")
+
+    time_chunks = out.chunksizes["time"]
+    assert all(c % 2 == 0 for c in time_chunks)  # aligned to source size 2
+    assert time_chunks[0] > 2  # genuinely coarsened
+    assert len(time_chunks) < 12  # fewer chunks than source partitions
+
+    xr.testing.assert_allclose(
+        out.compute().sortby(["time", "x"]),
+        ds.compute().sortby(["time", "x"]),
+    )
+
+
+# ---------------------------------------------------------------------------
+# dimension_columns / template resolution rules
+# ---------------------------------------------------------------------------
+
+
+def test_to_dataset_multi_registered_requires_explicit_template(
+    air_dataset_small,
+):
+    """With more than one registered Dataset, the caller disambiguates by
+    passing a registered table name as ``template=``."""
+    ctx = XarrayContext()
+    ctx.from_dataset("air1", air_dataset_small)
+    ctx.from_dataset("air2", air_dataset_small)
+    out = ctx.sql("SELECT * FROM air1").to_dataset(template="air1")
+    assert set(out.dims) == {"time", "lat", "lon"}
+
+
+def test_to_dataset_infer_fails_when_no_template_fits(air_dataset_small):
+    """If no registered Dataset's dims fit the result -> clear error."""
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    with pytest.raises(ValueError, match="dims cannot be inferred"):
+        ctx.sql(
+            "SELECT lat, lon, AVG(air) AS air_avg FROM air GROUP BY lat, lon"
+        ).to_dataset()
+
+
+def test_template_accepts_name_or_dataset(air_dataset_small):
+    """``template=`` accepts either a registered table name or a Dataset
+    object, with equivalent metadata recovery."""
+    other = air_dataset_small.copy()
+    other.attrs = {"flag": "other"}
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    ctx.from_dataset("other", other)
+
+    by_name = ctx.sql("SELECT * FROM air").to_dataset(
+        dims=["time", "lat", "lon"], template="other"
+    )
+    by_object = ctx.sql("SELECT * FROM air").to_dataset(
+        dims=["time", "lat", "lon"], template=other
+    )
+    assert by_name.attrs == {"flag": "other"}
+    assert by_object.attrs == {"flag": "other"}
+
+
+def test_template_unknown_name_raises(air_dataset_small):
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    with pytest.raises(ValueError, match="not a registered table"):
+        ctx.sql("SELECT * FROM air").to_dataset(
+            dims=["time", "lat", "lon"], template="missing"
+        )
+
+
+def test_template_recovers_var_encoding_strips_dtype(air_dataset_small):
+    """``zlib`` survives; dtype-bound keys are stripped (SQL may have cast)."""
+    ds = air_dataset_small.copy()
+    ds["air"].encoding = {
+        "zlib": True,
+        "dtype": "int16",
+        "_FillValue": -999,
+        "missing_value": -999,
+    }
+    ctx = XarrayContext()
+    ctx.from_dataset("air", ds)
+    out = ctx.sql("SELECT * FROM air").to_dataset(dims=["time", "lat", "lon"])
+    assert out["air"].encoding.get("zlib") is True
+    assert "dtype" not in out["air"].encoding
+    assert "_FillValue" not in out["air"].encoding
+    assert "missing_value" not in out["air"].encoding
+
+
+def test_template_aggregation_alias_no_attrs(air_dataset_small):
+    """``air_avg`` from ``AVG(air)`` does NOT inherit attrs from ``air``."""
+    ds = air_dataset_small.copy()
+    ds["air"].attrs = {"units": "K"}
+    ctx = XarrayContext()
+    ctx.from_dataset("air", ds)
+    out = ctx.sql(
+        "SELECT lat, lon, AVG(air) AS air_avg FROM air GROUP BY lat, lon"
+    ).to_dataset(dims=["lat", "lon"])
+    assert "air_avg" in out.data_vars
+    assert out["air_avg"].attrs == {}
+
+
+def test_to_dataset_explicit_template_overrides_auto_resolve(
+    air_dataset_small,
+):
+    """Explicit template= wins over the auto-resolved FROM-clause table."""
+    other = air_dataset_small.copy()
+    other.attrs = {"flag": "explicit"}
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    out = ctx.sql("SELECT * FROM air").to_dataset(
+        dims=["time", "lat", "lon"], template=other
+    )
+    assert out.attrs == {"flag": "explicit"}
+
+
+# ---------------------------------------------------------------------------
+# Lazy backend: value-level contract (not call counts)
+# ---------------------------------------------------------------------------
+
+
+def test_lazy_isel_int_round_trip(air_dataset_small):
+    """``isel(time=0)`` on the lazy result matches the eager equivalent."""
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    lazy = ctx.sql("SELECT * FROM air").to_dataset()
+    eager = lazy.compute()
+    actual = lazy["air"].isel(time=0).sortby(["lat", "lon"]).values
+    expected = eager["air"].isel(time=0).sortby(["lat", "lon"]).values
+    np.testing.assert_array_equal(actual, expected)
+
+
+def test_lazy_isel_slice_round_trip(air_dataset_small):
+    """isel(time=slice(0, 3)) round-trip matches the source."""
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    out = ctx.sql("SELECT * FROM air").to_dataset()
+    actual = out["air"].isel(time=slice(0, 3)).sortby(["lat", "lon"]).values
+    expected = (
+        air_dataset_small["air"]
+        .compute()
+        .isel(time=slice(0, 3))
+        .sortby(["lat", "lon"])
+        .values
+    )
+    np.testing.assert_array_equal(actual, expected)
+
+
+def test_lazy_outer_indexer_array(air_dataset_small):
+    """Fancy index along one dim works (IN-equivalent pushdown)."""
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    lazy = ctx.sql("SELECT * FROM air").to_dataset()
+    eager = lazy.compute()
+    indices = [0, 3, 5]
+    np.testing.assert_array_equal(
+        lazy["air"].isel(lat=indices).values,
+        eager["air"].isel(lat=indices).values,
+    )
+
+
+def test_lazy_compute_returns_eager(air_dataset_small):
+    """``.compute()`` returns an in-memory Dataset matching the source."""
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    out = ctx.sql("SELECT * FROM air").to_dataset().compute()
+    np.testing.assert_array_equal(
+        out.sortby(["time", "lat", "lon"])["air"].values,
+        air_dataset_small.compute()
+        .sortby(["time", "lat", "lon"])["air"]
+        .values,
+    )
+
+
+def test_vectorized_indexer_falls_back_via_xarray_adapter(
+    air_dataset_small,
+):
+    """VectorizedIndexer paths through xarray's adapter to outer + gather.
+
+    Our SQLBackendArray declares ``IndexingSupport.OUTER``, so xarray's
+    ``explicit_indexing_adapter`` converts vectorized indexers into a
+    series of outer reads followed by an in-memory numpy gather. The
+    public contract: values match the eager-computed equivalent.
+    """
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    lazy = ctx.sql("SELECT * FROM air").to_dataset()
+    eager = lazy.compute()
+
+    points_t = xr.DataArray([0, 3, 1], dims="point")
+    points_lat = xr.DataArray([2, 0, 5], dims="point")
+    np.testing.assert_array_equal(
+        lazy["air"].isel(time=points_t, lat=points_lat).values,
+        eager["air"].isel(time=points_t, lat=points_lat).values,
+    )
+
+
+# ---------------------------------------------------------------------------
+# Sparsity handling and fill_value
+# ---------------------------------------------------------------------------
+
+
+def test_sparsity_result_default_filters_lazy(air_dataset_small):
+    """Default sparsity='result' keeps only filtered coords (lazy path)."""
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    threshold = float(air_dataset_small["lat"].values[5])
+    out = ctx.sql(f"SELECT * FROM air WHERE lat > {threshold}").to_dataset()
+    assert (out["lat"].values > threshold).all()
+    assert out.sizes["lat"] < air_dataset_small.sizes["lat"]
+
+
+def test_sparsity_template_full_grid(air_dataset_small):
+    """sparsity='template' reindexes to the full grid with NaN fills."""
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    threshold = float(air_dataset_small["lat"].values[5])
+    out = ctx.sql(f"SELECT * FROM air WHERE lat > {threshold}").to_dataset(
+        sparsity="template"
+    )
+    assert out.sizes["lat"] == air_dataset_small.sizes["lat"]
+    lat_vals = out["lat"].values
+    below_mask = lat_vals <= threshold
+    above_mask = lat_vals > threshold
+    below = out["air"].isel(lat=np.where(below_mask)[0])
+    above = out["air"].isel(lat=np.where(above_mask)[0])
+    assert np.isnan(below.values).all()
+    assert not np.isnan(above.values).any()
+
+
+def test_sparsity_template_requires_template(air_dataset_small):
+    """No resolvable template -> sparsity='template' raises."""
+    other = air_dataset_small.copy()
+    ctx = XarrayContext()
+    ctx.from_dataset("a", air_dataset_small)
+    ctx.from_dataset("b", other)
+    with pytest.raises(ValueError, match="requires template= to be supplied"):
+        ctx.sql("SELECT * FROM a").to_dataset(
+            dims=["time", "lat", "lon"],
+            sparsity="template",
+        )
+
+
+def test_sparsity_invalid_value_raises(air_dataset_small):
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    with pytest.raises(ValueError, match="sparsity must be"):
+        ctx.sql("SELECT * FROM air").to_dataset(
+            dims=["time", "lat", "lon"],
+            sparsity="bogus",  # type: ignore[arg-type]
+        )
+
+
+def test_sparsity_template_with_aggregation(air_dataset_small):
+    """sparsity='template' on an aggregation respects dimension_columns subset."""
+    ctx = XarrayContext()
+    ctx.from_dataset("air", air_dataset_small)
+    threshold = float(air_dataset_small["lat"].values[5])
+    out = ctx.sql(
+        f"""
+        SELECT lat, lon, AVG(air) AS air_avg
+        FROM air
+        WHERE lat > {threshold}
+        GROUP BY lat, lon
+        """
+    ).to_dataset(dims=["lat", "lon"], sparsity="template")
+    assert out.sizes["lat"] == air_dataset_small.sizes["lat"]
+    assert "time" not in out.dims
+    below_mask = out["lat"].values <= threshold
+    below = out["air_avg"].isel(lat=np.where(below_mask)[0])
+    assert np.isnan(below.values).all()
+
+
+def test_fill_value_int_upcasts_to_float():
+    """fill_value=NaN forces float upcast on int columns -- documented."""
+    ds = xr.Dataset(
+        {"v": (("lat", "lon"), np.arange(6, dtype=np.int64).reshape(3, 2))},
+        coords={"lat": [0, 1, 2], "lon": [10, 11]},
+    ).chunk({"lat": 3})
+    ctx = XarrayContext()
+    ctx.from_dataset("t", ds)
+    out = ctx.sql("SELECT * FROM t WHERE lat > 0").to_dataset(
+        sparsity="template"
+    )
+    assert np.issubdtype(out["v"].dtype, np.floating)
+    assert np.isnan(out["v"].sel(lat=0).values).all()
+
+
+def test_fill_value_custom_preserves_int(air_dataset_small):
+    """Passing a typed sentinel preserves the data var's int dtype."""
+    source = xr.Dataset(
+        {
+            "v": (
+                ("lat", "lon"),
+                np.arange(6, dtype=np.int64).reshape(3, 2) + 1,
+            ),
+        },
+        coords={"lat": [0, 1, 2], "lon": [10, 11]},
+    ).chunk({"lat": 3})
+    ctx = XarrayContext()
+    ctx.from_dataset("t", source)
+    out = ctx.sql("SELECT * FROM t WHERE lat > 0").to_dataset(
+        sparsity="template", fill_value=-1
+    )
+    assert np.issubdtype(out["v"].dtype, np.integer)
+    assert (out["v"].sel(lat=0).values == -1).all()
+    assert out["v"].sel(lat=2, lon=11).item() == 6
+
+
+def test_sparsity_template_then_metadata(air_dataset_small):
+    """sparsity='template' composes with template metadata recovery."""
+    ds = air_dataset_small.copy()
+    ds.attrs = {"src": "tmpl"}
+    ds["air"].attrs = {"units": "K"}
+    ctx = XarrayContext()
+    ctx.from_dataset("air", ds)
+    threshold = float(ds["lat"].values[5])
+    out = ctx.sql(f"SELECT * FROM air WHERE lat > {threshold}").to_dataset(
+        sparsity="template"
+    )
+    assert out.attrs == {"src": "tmpl"}
+    assert out["air"].attrs == {"units": "K"}
+    assert out.sizes["lat"] == ds.sizes["lat"]
diff --git a/uv.lock b/uv.lock
index fa753dc..1e80bb0 100644
--- a/uv.lock
+++ b/uv.lock
@@ -168,40 +168,6 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/50/cd/30110dc0ffcf3b131156077b90e9f60ed75711223f306da4db08eff8403b/beautifulsoup4-4.13.4-py3-none-any.whl", hash = "sha256:9bbbb14bfde9d79f38b8cd5f8c7c85f4b8f2523190ebed90e950a8dea4cb1c4b", size = 187285, upload-time = "2025-04-15T17:05:12.221Z" },
 ]
 
-[[package]]
-name = "black"
-version = "24.10.0"
-source = { registry = "https://pypi.org/simple" }
-dependencies = [
-    { name = "click" },
-    { name = "mypy-extensions" },
-    { name = "packaging" },
-    { name = "pathspec" },
-    { name = "platformdirs" },
-    { name = "tomli", marker = "python_full_version < '3.11'" },
-    { name = "typing-extensions", marker = "python_full_version < '3.11'" },
-]
-sdist = { url = "https://files.pythonhosted.org/packages/d8/0d/cc2fb42b8c50d80143221515dd7e4766995bd07c56c9a3ed30baf080b6dc/black-24.10.0.tar.gz", hash = "sha256:846ea64c97afe3bc677b761787993be4991810ecc7a4a937816dd6bddedc4875", size = 645813, upload-time = "2024-10-07T19:20:50.361Z" }
-wheels = [
-    { url = "https://files.pythonhosted.org/packages/a3/f3/465c0eb5cddf7dbbfe1fecd9b875d1dcf51b88923cd2c1d7e9ab95c6336b/black-24.10.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:e6668650ea4b685440857138e5fe40cde4d652633b1bdffc62933d0db4ed9812", size = 1623211, upload-time = "2024-10-07T19:26:12.43Z" },
-    { url = "https://files.pythonhosted.org/packages/df/57/b6d2da7d200773fdfcc224ffb87052cf283cec4d7102fab450b4a05996d8/black-24.10.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:1c536fcf674217e87b8cc3657b81809d3c085d7bf3ef262ead700da345bfa6ea", size = 1457139, upload-time = "2024-10-07T19:25:06.453Z" },
-    { url = "https://files.pythonhosted.org/packages/6e/c5/9023b7673904a5188f9be81f5e129fff69f51f5515655fbd1d5a4e80a47b/black-24.10.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:649fff99a20bd06c6f727d2a27f401331dc0cc861fb69cde910fe95b01b5928f", size = 1753774, upload-time = "2024-10-07T19:23:58.47Z" },
-    { url = "https://files.pythonhosted.org/packages/e1/32/df7f18bd0e724e0d9748829765455d6643ec847b3f87e77456fc99d0edab/black-24.10.0-cp310-cp310-win_amd64.whl", hash = "sha256:fe4d6476887de70546212c99ac9bd803d90b42fc4767f058a0baa895013fbb3e", size = 1414209, upload-time = "2024-10-07T19:24:42.54Z" },
-    { url = "https://files.pythonhosted.org/packages/c2/cc/7496bb63a9b06a954d3d0ac9fe7a73f3bf1cd92d7a58877c27f4ad1e9d41/black-24.10.0-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:5a2221696a8224e335c28816a9d331a6c2ae15a2ee34ec857dcf3e45dbfa99ad", size = 1607468, upload-time = "2024-10-07T19:26:14.966Z" },
-    { url = "https://files.pythonhosted.org/packages/2b/e3/69a738fb5ba18b5422f50b4f143544c664d7da40f09c13969b2fd52900e0/black-24.10.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:f9da3333530dbcecc1be13e69c250ed8dfa67f43c4005fb537bb426e19200d50", size = 1437270, upload-time = "2024-10-07T19:25:24.291Z" },
-    { url = "https://files.pythonhosted.org/packages/c9/9b/2db8045b45844665c720dcfe292fdaf2e49825810c0103e1191515fc101a/black-24.10.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:4007b1393d902b48b36958a216c20c4482f601569d19ed1df294a496eb366392", size = 1737061, upload-time = "2024-10-07T19:23:52.18Z" },
-    { url = "https://files.pythonhosted.org/packages/a3/95/17d4a09a5be5f8c65aa4a361444d95edc45def0de887810f508d3f65db7a/black-24.10.0-cp311-cp311-win_amd64.whl", hash = "sha256:394d4ddc64782e51153eadcaaca95144ac4c35e27ef9b0a42e121ae7e57a9175", size = 1423293, upload-time = "2024-10-07T19:24:41.7Z" },
-    { url = "https://files.pythonhosted.org/packages/90/04/bf74c71f592bcd761610bbf67e23e6a3cff824780761f536512437f1e655/black-24.10.0-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:b5e39e0fae001df40f95bd8cc36b9165c5e2ea88900167bddf258bacef9bbdc3", size = 1644256, upload-time = "2024-10-07T19:27:53.355Z" },
-    { url = "https://files.pythonhosted.org/packages/4c/ea/a77bab4cf1887f4b2e0bce5516ea0b3ff7d04ba96af21d65024629afedb6/black-24.10.0-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:d37d422772111794b26757c5b55a3eade028aa3fde43121ab7b673d050949d65", size = 1448534, upload-time = "2024-10-07T19:26:44.953Z" },
-    { url = "https://files.pythonhosted.org/packages/4e/3e/443ef8bc1fbda78e61f79157f303893f3fddf19ca3c8989b163eb3469a12/black-24.10.0-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:14b3502784f09ce2443830e3133dacf2c0110d45191ed470ecb04d0f5f6fcb0f", size = 1761892, upload-time = "2024-10-07T19:24:10.264Z" },
-    { url = "https://files.pythonhosted.org/packages/52/93/eac95ff229049a6901bc84fec6908a5124b8a0b7c26ea766b3b8a5debd22/black-24.10.0-cp312-cp312-win_amd64.whl", hash = "sha256:30d2c30dc5139211dda799758559d1b049f7f14c580c409d6ad925b74a4208a8", size = 1434796, upload-time = "2024-10-07T19:25:06.239Z" },
-    { url = "https://files.pythonhosted.org/packages/d0/a0/a993f58d4ecfba035e61fca4e9f64a2ecae838fc9f33ab798c62173ed75c/black-24.10.0-cp313-cp313-macosx_10_13_x86_64.whl", hash = "sha256:1cbacacb19e922a1d75ef2b6ccaefcd6e93a2c05ede32f06a21386a04cedb981", size = 1643986, upload-time = "2024-10-07T19:28:50.684Z" },
-    { url = "https://files.pythonhosted.org/packages/37/d5/602d0ef5dfcace3fb4f79c436762f130abd9ee8d950fa2abdbf8bbc555e0/black-24.10.0-cp313-cp313-macosx_11_0_arm64.whl", hash = "sha256:1f93102e0c5bb3907451063e08b9876dbeac810e7da5a8bfb7aeb5a9ef89066b", size = 1448085, upload-time = "2024-10-07T19:28:12.093Z" },
-    { url = "https://files.pythonhosted.org/packages/47/6d/a3a239e938960df1a662b93d6230d4f3e9b4a22982d060fc38c42f45a56b/black-24.10.0-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:ddacb691cdcdf77b96f549cf9591701d8db36b2f19519373d60d31746068dbf2", size = 1760928, upload-time = "2024-10-07T19:24:15.233Z" },
-    { url = "https://files.pythonhosted.org/packages/dd/cf/af018e13b0eddfb434df4d9cd1b2b7892bab119f7a20123e93f6910982e8/black-24.10.0-cp313-cp313-win_amd64.whl", hash = "sha256:680359d932801c76d2e9c9068d05c6b107f2584b2a5b88831c83962eb9984c1b", size = 1436875, upload-time = "2024-10-07T19:24:42.762Z" },
-    { url = "https://files.pythonhosted.org/packages/8d/a7/4b27c50537ebca8bec139b872861f9d2bf501c5ec51fcf897cb924d9e264/black-24.10.0-py3-none-any.whl", hash = "sha256:3bb2b7a1f7b685f85b11fed1ef10f8a9148bceb49853e47a294a3dd963c1dd7d", size = 206898, upload-time = "2024-10-07T19:20:48.317Z" },
-]
-
 [[package]]
 name = "cachetools"
 version = "5.5.2"
@@ -1322,15 +1288,6 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/d8/30/9aec301e9772b098c1f5c0ca0279237c9766d94b97802e9888010c64b0ed/multidict-6.6.3-py3-none-any.whl", hash = "sha256:8db10f29c7541fc5da4defd8cd697e1ca429db743fa716325f236079b96f775a", size = 12313, upload-time = "2025-06-30T15:53:45.437Z" },
 ]
 
-[[package]]
-name = "mypy-extensions"
-version = "1.1.0"
-source = { registry = "https://pypi.org/simple" }
-sdist = { url = "https://files.pythonhosted.org/packages/a2/6e/371856a3fb9d31ca8dac321cda606860fa4548858c0cc45d9d1d4ca2628b/mypy_extensions-1.1.0.tar.gz", hash = "sha256:52e68efc3284861e772bbcd66823fde5ae21fd2fdb51c62a211403730b916558", size = 6343, upload-time = "2025-04-22T14:54:24.164Z" }
-wheels = [
-    { url = "https://files.pythonhosted.org/packages/79/7b/2c79738432f5c924bef5071f933bcc9efd0473bac3b4aa584a6f7c1c8df8/mypy_extensions-1.1.0-py3-none-any.whl", hash = "sha256:1be4cccdb0f2482337c4743e60421de3a356cd97508abadd57d47403e94f5505", size = 4963, upload-time = "2025-04-22T14:54:22.983Z" },
-]
-
 [[package]]
 name = "netcdf4"
 version = "1.7.2"
@@ -1953,25 +1910,6 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/c7/21/705964c7812476f378728bdf590ca4b771ec72385c533964653c68e86bdc/pygments-2.19.2-py3-none-any.whl", hash = "sha256:86540386c03d588bb81d44bc3928634ff26449851e99741617ecb9037ee5ec0b", size = 1225217, upload-time = "2025-06-21T13:39:07.939Z" },
 ]
 
-[[package]]
-name = "pyink"
-version = "24.10.1"
-source = { registry = "https://pypi.org/simple" }
-dependencies = [
-    { name = "black" },
-    { name = "click" },
-    { name = "mypy-extensions" },
-    { name = "packaging" },
-    { name = "pathspec" },
-    { name = "platformdirs" },
-    { name = "tomli", marker = "python_full_version < '3.11'" },
-    { name = "typing-extensions", marker = "python_full_version < '3.11'" },
-]
-sdist = { url = "https://files.pythonhosted.org/packages/d1/a1/e5e28626fca4266a94c2e1c9264fbf915b9e83e94f52e965190e48fd0cbf/pyink-24.10.1.tar.gz", hash = "sha256:5ec4339aa4953f796e88d90bcac3e3472161e4c36dbde203d80f5f76721ac718", size = 267230, upload-time = "2025-01-10T11:28:09.907Z" }
-wheels = [
-    { url = "https://files.pythonhosted.org/packages/1f/12/2f271b3601ae25731879f160d6b3941d80eb6b4f3e24be90289e33fb1dc4/pyink-24.10.1-py3-none-any.whl", hash = "sha256:6349bf6ab75e2ea39a5f0bc3dee7ede7f4af8529291472638026de5fd4af80d2", size = 137118, upload-time = "2025-01-10T11:28:06.138Z" },
-]
-
 [[package]]
 name = "pymdown-extensions"
 version = "10.21"
@@ -2150,6 +2088,31 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/64/8d/0133e4eb4beed9e425d9a98ed6e081a55d195481b7632472be1af08d2f6b/rsa-4.9.1-py3-none-any.whl", hash = "sha256:68635866661c6836b8d39430f97a996acbd61bfa49406748ea243539fe239762", size = 34696, upload-time = "2025-04-16T09:51:17.142Z" },
 ]
 
+[[package]]
+name = "ruff"
+version = "0.15.17"
+source = { registry = "https://pypi.org/simple" }
+sdist = { url = "https://files.pythonhosted.org/packages/8c/a9/3abdf488f1bf3d24c699415e454ed554a6350d5d89ce183be1ee0a3361ac/ruff-0.15.17.tar.gz", hash = "sha256:2ec446937fd16c8c4de2674a209cc5af64d9c6f17d21fbf1151054fa0bcf5219", size = 4743346, upload-time = "2026-06-11T17:54:47.663Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/db/4d/e11259f5da07cb6afb2d074c31bf09da9671993f7329d4f15d2fdc458301/ruff-0.15.17-py3-none-linux_armv6l.whl", hash = "sha256:d9feddb927fc68bd295f5eebc587a7e42cfaf9b65f60ca4a2386febff575da8f", size = 10856677, upload-time = "2026-06-11T17:54:49.533Z" },
+    { url = "https://files.pythonhosted.org/packages/29/3e/772d679e1a0dc058e58875bd2c0cb713a0530877b4a76fee3c7966df0d49/ruff-0.15.17-py3-none-macosx_10_12_x86_64.whl", hash = "sha256:25805a226d741c47d274a35ad5c10a7dde175fcddfa511d7cf3da0a21eb3eab7", size = 11223443, upload-time = "2026-06-11T17:55:00.573Z" },
+    { url = "https://files.pythonhosted.org/packages/68/58/bd41f7688b2fd5623012605130ed70e60aa7f2244baa3d5066bdd61530c8/ruff-0.15.17-py3-none-macosx_11_0_arm64.whl", hash = "sha256:f6ad73b14c2d18a3bf8ad7cb6974294d7f613a7898604826058e6ac64918ef4d", size = 10566458, upload-time = "2026-06-11T17:55:07.52Z" },
+    { url = "https://files.pythonhosted.org/packages/d8/5b/733371013fcf1ec339e477ece6ab42bfe10bdd9bba8ee88a9516aa56bfc0/ruff-0.15.17-py3-none-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:6ba0c1e4f95bcb3869d0d30cbd5917071ef2e28665abfec970cdab0492c713ed", size = 10914483, upload-time = "2026-06-11T17:55:05.501Z" },
+    { url = "https://files.pythonhosted.org/packages/bd/cc/6f24251cc0252f7239391ccb85833f320efad14ebe5b443943f37ced6332/ruff-0.15.17-py3-none-manylinux_2_17_armv7l.manylinux2014_armv7l.whl", hash = "sha256:81647960f10bff57d2e51cadd0c3950fe598400c852863a038720ef5b8cca91e", size = 10647497, upload-time = "2026-06-11T17:54:57.733Z" },
+    { url = "https://files.pythonhosted.org/packages/68/dd/0d10c17ce1a1624d6fc3156309c3f834fdb5dfaad026ec90c85684f3990e/ruff-0.15.17-py3-none-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:0e01a84ddbc8c16c23055ba3924476850f1bbc1917cebbb9376665a63e74260d", size = 11416967, upload-time = "2026-06-11T17:54:51.461Z" },
+    { url = "https://files.pythonhosted.org/packages/2f/91/556bfb156f6144f355e831c23db00b2fc4120f86b3ce81cc5f7fd2df51f3/ruff-0.15.17-py3-none-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:84fe9f653152f8f294f9f7e03bf3a453d8b4a27f7a59c78c8666167f2b17b96c", size = 12335770, upload-time = "2026-06-11T17:54:45.793Z" },
+    { url = "https://files.pythonhosted.org/packages/88/82/8b5999aa13355e926f06d9f42a32dcca862f623bf0363785ff89d607dffd/ruff-0.15.17-py3-none-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:8c0fe88a7676e7a05b73174d4d4a59cb2ac21ff8263583f87a81a6018475a978", size = 11575441, upload-time = "2026-06-11T17:54:32.661Z" },
+    { url = "https://files.pythonhosted.org/packages/11/93/f10377bb04109ca0e8cbc483ff1982c54b6d418210041776f93e8cdc7fa9/ruff-0.15.17-py3-none-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ecfc3c7878fff94633ab0348524e093f9ce3243080416dd7d14f8ba400174719", size = 11557614, upload-time = "2026-06-11T17:54:34.698Z" },
+    { url = "https://files.pythonhosted.org/packages/c7/a6/eeeae7f7d5493df41649ab3db92f086b2d0a30199e4efdf8e3dd7a033f24/ruff-0.15.17-py3-none-manylinux_2_31_riscv64.whl", hash = "sha256:b8461180b22420b1bdc289909410930761629fddf2a5aaf60fae1ab26cedc4c4", size = 11544450, upload-time = "2026-06-11T17:54:39.042Z" },
+    { url = "https://files.pythonhosted.org/packages/32/88/5991ce565129a24dd4a00db1254b3b5db2e53018cbe4018ea5a89738e727/ruff-0.15.17-py3-none-musllinux_1_2_aarch64.whl", hash = "sha256:6eccbe50a038b503e7140b441aa9c7fc8c1f36edf23ebef9f4165c2f28f568b7", size = 10892524, upload-time = "2026-06-11T17:55:09.432Z" },
+    { url = "https://files.pythonhosted.org/packages/f5/1d/0fdd248313425f55223968af04b0a42125466a8d88d21c1d99c6af0a51e8/ruff-0.15.17-py3-none-musllinux_1_2_armv7l.whl", hash = "sha256:382fc0521025f5a8ad447d8bdd523545d0d7646adb718eb1c2dac5065ec27c0f", size = 10659573, upload-time = "2026-06-11T17:54:36.824Z" },
+    { url = "https://files.pythonhosted.org/packages/9e/0e/072e8260deb9461062ce9311ced27a8e541229a6ffd483013dd37661e43e/ruff-0.15.17-py3-none-musllinux_1_2_i686.whl", hash = "sha256:456d41fcd1b2777ad63f09a6e7121d43f7b688bbc76a800c10f7f8fb1f912c3f", size = 11127818, upload-time = "2026-06-11T17:55:03.124Z" },
+    { url = "https://files.pythonhosted.org/packages/ab/b4/55060a34163121498014696b5f656db5b8c6963768f227dbf0d76b311073/ruff-0.15.17-py3-none-musllinux_1_2_x86_64.whl", hash = "sha256:b1a04bcc94ae6194e9db05d16ad31f298a7194bfbcb08258bbe589cee1d587b8", size = 11655901, upload-time = "2026-06-11T17:54:53.562Z" },
+    { url = "https://files.pythonhosted.org/packages/49/71/9b29d6b87cef468d697f43c6a91e3fae4a80185779d7d5a4ef27d173439f/ruff-0.15.17-py3-none-win32.whl", hash = "sha256:596065960ab1ff593f744220c9fe6580eda00a95003cffa9f4048bb5b1bf0392", size = 10925574, upload-time = "2026-06-11T17:54:55.723Z" },
+    { url = "https://files.pythonhosted.org/packages/3d/b2/8fc77f3723228836fa5d12497eb71c808f83782e10d058d2b15cfa14640b/ruff-0.15.17-py3-none-win_amd64.whl", hash = "sha256:6769e5fa1710b179b92e0bfa5a51735b35baea9013dadb06d5f44cbcf9547084", size = 12058788, upload-time = "2026-06-11T17:54:41.042Z" },
+    { url = "https://files.pythonhosted.org/packages/2d/c7/c53e8dbff9c9dc4b7928773421ae294a5d28fcb8dcda1a089579d3a7e510/ruff-0.15.17-py3-none-win_arm64.whl", hash = "sha256:f3be1fbb34bcdfd146240d8fb92a709d4c2c8191348580a3c044ec60fa0b4456", size = 11355275, upload-time = "2026-06-11T17:54:43.635Z" },
+]
+
 [[package]]
 name = "scipy"
 version = "1.15.3"
@@ -2624,7 +2587,7 @@ test = [
 dev = [
     { name = "maturin" },
     { name = "py-spy" },
-    { name = "pyink" },
+    { name = "ruff" },
     { name = "xarray-sql", extra = ["docs", "test"] },
 ]
 
@@ -2650,7 +2613,7 @@ provides-extras = ["dev", "docs", "test"]
 dev = [
     { name = "maturin", specifier = ">=1.9.1" },
     { name = "py-spy", specifier = ">=0.4.0" },
-    { name = "pyink", specifier = ">=24.10.1" },
+    { name = "ruff", specifier = ">=0.15.10" },
     { name = "xarray-sql", extras = ["docs"] },
     { name = "xarray-sql", extras = ["test"] },
 ]
diff --git a/xarray_sql/__init__.py b/xarray_sql/__init__.py
index 75b889b..d1e5984 100644
--- a/xarray_sql/__init__.py
+++ b/xarray_sql/__init__.py
@@ -1,7 +1,7 @@
 from . import cftime
+from .df import from_map
 from .reader import read_xarray, read_xarray_table
 from .sql import XarrayContext
-from .df import from_map
 
 __all__ = [
     "cftime",
diff --git a/xarray_sql/ds.py b/xarray_sql/ds.py
new file mode 100644
index 0000000..c33f996
--- /dev/null
+++ b/xarray_sql/ds.py
@@ -0,0 +1,838 @@
+"""Reconstruct xarray Datasets from SQL query results.
+
+The inverse of the forward Dataset-to-table pivot done by
+:func:`xarray_sql.df.pivot`. Internally defines an :class:`XarrayDataFrame`
+wrapper around the DataFusion ``DataFrame`` returned by
+:meth:`XarrayContext.sql`, with a :meth:`XarrayDataFrame.to_dataset`
+method that round-trips a query result back to ``xr.Dataset``.
+
+Reconstruction is controlled by the ``chunks`` argument to
+:meth:`XarrayDataFrame.to_dataset` -- the xarray idiom for tuning how a
+result is partitioned -- rather than by reflecting on the query plan:
+
+* **Eager** (``chunks=None``, or the default ``"inherit"`` when the
+  result keeps no multi-chunk source dimension): the plan executes
+  exactly once via ``execute_stream`` and the result is scattered into a
+  dense in-memory Dataset. This is the right default for reductions
+  (aggregations), whose results are small, and it never re-executes.
+* **Lazy / chunked** (``chunks`` is a mapping, ``"auto"``, or
+  ``"inherit"`` over a multi-chunk source dimension): data variables are
+  backed by :class:`SQLBackendArray` wrapped in
+  ``xarray.core.indexing.LazilyIndexedArray`` and chunked via xarray's
+  configured chunk manager (dask, cubed, ...). Each chunk maps onto the
+  source partitions and reads its coordinate range on access by
+  translating the indexer into a DataFusion ``filter`` expression, so only
+  the requested partitions are materialized as Arrow ``RecordBatch`` es
+  and scattered into numpy.
+
+``.compute()`` materializes the whole Dataset in memory.
+"""
+
+from __future__ import annotations
+
+import warnings
+from collections.abc import Mapping
+from typing import Any, Literal, cast
+
+import numpy as np
+import pandas as pd
+import pyarrow as pa
+import xarray as xr
+from datafusion import col, literal
+
+Sparsity = Literal["result", "template"]
+"""Output coordinate extent for a filtered round-trip.
+
+* ``"result"`` keeps only the dim values present in the query result, so
+  the output is sparse and equal to whatever rows came back.
+* ``"template"`` reindexes to the registered Dataset's full coord ranges
+  and fills absent cells with ``fill_value``.
+"""
+
+
+# ---------------------------------------------------------------------------
+# Private helpers
+# ---------------------------------------------------------------------------
+
+
+def _ds_var_dims(ds: xr.Dataset) -> list[str]:
+    """Return a Dataset's data-variable dim order.
+
+    The forward path validates that all data variables share the same dims
+    tuple, so the first var's dim order is canonical. Falls back to
+    ``ds.dims`` keys for empty Datasets. Always use this rather than
+    ``list(ds.dims)`` when round-tripping, since the latter is in
+    canonical name order and may not match the variable's axis order.
+    """
+    if ds.data_vars:
+        return list(next(iter(ds.data_vars.values())).dims)
+    return list(ds.dims)
+
+
+def _apply_template(ds: xr.Dataset, template: xr.Dataset) -> xr.Dataset:
+    """Recover metadata that the forward SQL pivot strips.
+
+    Adds back, where unambiguous:
+
+    * Data-variable ``attrs`` and ``encoding`` for vars present in
+      ``template`` (aggregation aliases like ``air_avg`` get nothing).
+      Dtype-bound encoding keys (``dtype``, ``_FillValue``,
+      ``missing_value``) are intentionally dropped: SQL may have
+      changed the column's dtype (e.g. ``int16`` -> ``float64`` after
+      ``AVG`` or a null-introducing filter), and reattaching the
+      source's packing would make a later ``ds.to_netcdf()`` write
+      corrupt values.
+    * Dim-coordinate dtype, where SQL upcasted (datetime is the
+      canonical case).
+    * Non-dim coordinates whose dims are all present in ``ds`` (scalar
+      coords attach as-is; vector coords use ``.sel``).
+    * Dataset-level ``attrs``.
+
+    Skipped coords are warned about once per call.
+    """
+    out = ds.copy()
+
+    # 1. Data-var attrs / encoding for vars present in the template.
+    #    Aggregation aliases absent from template intentionally inherit nothing.
+    for name in list(out.data_vars):
+        if name in template.data_vars:
+            out[name].attrs = dict(template[name].attrs)
+            # Drop dtype-bound encoding keys; SQL may have changed dtype.
+            enc = {
+                k: v
+                for k, v in template[name].encoding.items()
+                if k not in {"dtype", "_FillValue", "missing_value"}
+            }
+            out[name].encoding = enc
+
+    # 2. Restore dim-coordinate dtype when SQL changed it (e.g. datetime
+    #    upcast through pyarrow / pandas) and copy the source's dim-coord
+    #    attrs (``standard_name``, ``long_name``, ``units``, etc.).
+    for d in list(out.dims):
+        if d in template.coords:
+            tdt = template.coords[d].dtype
+            if out.coords[d].dtype != tdt:
+                try:
+                    out = out.assign_coords({d: out.coords[d].astype(tdt)})
+                except (ValueError, TypeError):
+                    pass  # incompatible cast; leave as-is
+            out[d].attrs = dict(template.coords[d].attrs)
+
+    # 3. Non-dim coordinates whose dims are all present in the result.
+    out_dims = set(out.dims)
+    skipped: list[str] = []
+    for cname, coord in template.coords.items():
+        if cname in template.dims:
+            continue  # dim coord; already in out
+        if not set(coord.dims) <= out_dims:
+            continue  # spans dims the result lacks
+        try:
+            if not coord.dims:
+                # Scalar coord (e.g. weather_dataset.reference_time).
+                out = out.assign_coords({cname: coord})
+            else:
+                sel = {d: out.coords[d] for d in coord.dims}
+                out = out.assign_coords({cname: coord.sel(sel)})
+        except (KeyError, ValueError, TypeError):
+            skipped.append(cname)
+
+    # 4. Dataset-level attrs.
+    out.attrs = dict(template.attrs)
+
+    if skipped:
+        warnings.warn(
+            f"Could not re-attach non-dim coordinates from template: {skipped}",
+            stacklevel=3,
+        )
+    return out
+
+
+def _scatter_batches_to_ndarray(
+    batches: list[pa.RecordBatch],
+    dimension_columns: list[str],
+    requested: dict[str, np.ndarray],
+    var_name: str,
+    out_shape: tuple[int, ...],
+    dtype: np.dtype,
+    drop_axes: list[int],
+) -> np.ndarray:
+    """Convert filtered Arrow ``RecordBatch`` rows into a dense N-D numpy array.
+
+    SQL query results arrive as flat rows; xarray expects N-D arrays.
+    This bridges the two: each row carries the dim-coord values that
+    identify its cell in the output cube plus the value to write there.
+    We look up the row's N-D position by binary-searching its coord
+    values within the caller's requested coord arrays
+    (``np.searchsorted``), then scatter-write the value at that index.
+
+    Missing combinations (sparse results from filtered queries) stay as
+    ``NaN`` for floating-point outputs by pre-filling the buffer; integer
+    outputs leave them as ``np.empty``-style undefined values.
+    """
+    # NaN fill for float outputs; default for int/datetime falls through
+    # to ``np.empty``-style undefined values (but every output cell is
+    # written below for non-sparse cases).
+    out = (
+        np.full(out_shape, np.nan, dtype=dtype)
+        if np.issubdtype(dtype, np.floating)
+        else np.empty(out_shape, dtype=dtype)
+    )
+
+    # ``requested[d]`` may be in any order (callers can iselect arbitrary
+    # positions, and template coords like air_temperature.lat are descending).
+    # ``np.searchsorted`` requires ascending input, so we sort each requested
+    # array once, search there, and remap back to the original positions.
+    sorted_idx = {d: np.argsort(requested[d]) for d in dimension_columns}
+    sorted_req = {d: requested[d][sorted_idx[d]] for d in dimension_columns}
+
+    for batch in batches:
+        if batch.num_rows == 0:
+            continue
+        schema_names = batch.schema.names
+        # Build per-dim position arrays for this batch (positions within
+        # the caller's requested coord order).
+        positions = []
+        for d in dimension_columns:
+            col_arr = batch.column(schema_names.index(d))
+            vals = col_arr.to_numpy(zero_copy_only=False)
+            pos_in_sorted = np.searchsorted(sorted_req[d], vals)
+            positions.append(sorted_idx[d][pos_in_sorted])
+        value_arr = batch.column(schema_names.index(var_name)).to_numpy(
+            zero_copy_only=False
+        )
+        out[tuple(positions)] = value_arr.astype(dtype, copy=False)
+
+    if drop_axes:
+        out = np.squeeze(out, axis=tuple(drop_axes))
+    return cast(np.ndarray, out)
+
+
+class SQLBackendArray(xr.backends.BackendArray):
+    """Read-only lazy N-D array view over a DataFusion DataFrame.
+
+    Bridges xarray's lazy-indexing interface
+    (:class:`xarray.backends.BackendArray`) to a DataFusion query result,
+    so an xarray ``Dataset`` can present a SQL query as if it were a
+    materialized N-D array without actually loading any data until the
+    caller asks for it. This is the workhorse that lets
+    :meth:`XarrayDataFrame.to_dataset` return a Dataset cheaply.
+
+    On each ``__getitem__`` call, the requested xarray indexer is
+    translated into a DataFusion filter expression (``df.filter(expr)``)
+    and a column projection (``df.select(*cols)``). The filtered
+    DataFrame is consumed via ``execute_stream`` as a sequence of Arrow
+    ``RecordBatch`` es and scattered into a preallocated numpy buffer,
+    so only the requested data is materialized.
+
+    Constraints and caveats:
+
+    - Read-only: there is no write path; the backend exists to surface
+      query results, not to round-trip writes into a SQL store.
+    - The underlying DataFusion ``DataFrame`` holds a reference to its
+      originating ``SessionContext``, which is not picklable. The class
+      therefore overrides ``__copy__`` and ``__deepcopy__`` to return
+      ``self`` -- this is safe because the backend is read-only.
+    - ``IndexingSupport.OUTER``: ``BasicIndexer`` and ``OuterIndexer``
+      are translated to filter predicates directly; ``VectorizedIndexer``
+      paths through xarray's adapter to outer-then-gather and so still
+      works, just less efficiently.
+
+    Raises:
+        ValueError, datafusion exceptions: propagated from the
+            underlying ``df.filter().select().execute_stream()`` chain
+            if a predicate refers to a missing column, the dtype of a
+            literal is incompatible, or the execution itself fails.
+        AssertionError: from ``np.searchsorted`` mis-alignment, which
+            indicates the result contains coordinate values not present
+            in the wrapper's pre-computed coord arrays -- usually a
+            symptom of a filtered query whose coord discovery missed a
+            value.
+
+    Constructed by :func:`_build_lazy_scan`; users should not instantiate
+    this class directly.
+    """
+
+    def __init__(
+        self,
+        inner_df: Any,
+        var_name: str,
+        dimension_columns: list[str],
+        coord_arrays: dict[str, np.ndarray],
+        shape: tuple[int, ...],
+        dtype: np.dtype,
+    ) -> None:
+        self._inner_df = inner_df
+        self._var_name = var_name
+        self._dimension_columns = list(dimension_columns)
+        self._coord_arrays = coord_arrays
+        self.shape = tuple(shape)
+        self.dtype = np.dtype(dtype)
+
+    def __getitem__(self, key: Any) -> np.ndarray:
+        return cast(
+            np.ndarray,
+            xr.core.indexing.explicit_indexing_adapter(
+                key,
+                self.shape,
+                xr.core.indexing.IndexingSupport.OUTER,
+                self._raw_getitem,
+            ),
+        )
+
+    def __copy__(self) -> "SQLBackendArray":
+        # The backend is read-only; the underlying DataFusion DataFrame
+        # holds a non-picklable SessionContext reference, so sharing the
+        # same backend across a copy is both safe and necessary.
+        return self
+
+    def __deepcopy__(self, memo: dict) -> "SQLBackendArray":
+        return self
+
+    # ------------------------------------------------------------------
+
+    def _raw_getitem(self, key: tuple) -> np.ndarray:
+        """Materialize the indexed region described by *key* via DataFusion + Arrow.
+
+        ``key`` is a tuple of ``int``/``slice``/1-D integer-array, one per
+        dim, in :attr:`_dimension_columns` order.
+        """
+        requested: dict[str, np.ndarray] = {}
+        # Dims whose indexer covers the full extent (slice(None) or
+        # equivalent). For these we omit the filter predicate entirely
+        # so DataFusion doesn't have to evaluate a tautology.
+        full_dims: set[str] = set()
+        drop_axes: list[int] = []
+        for axis, (dim, k) in enumerate(
+            zip(self._dimension_columns, key, strict=True)
+        ):
+            coord = self._coord_arrays[dim]
+            if isinstance(k, slice):
+                start = 0 if k.start is None else k.start
+                stop = len(coord) if k.stop is None else k.stop
+                step = 1 if k.step is None else k.step
+                requested[dim] = np.asarray(coord[start:stop:step])
+                if start == 0 and stop >= len(coord) and step == 1:
+                    full_dims.add(dim)
+            elif isinstance(k, (int, np.integer)):
+                requested[dim] = np.asarray([coord[int(k)]])
+                drop_axes.append(axis)
+            else:
+                arr = np.asarray(k)
+                requested[dim] = np.asarray(coord[arr])
+                if (
+                    len(arr) == len(coord)
+                    and (arr == np.arange(len(coord))).all()
+                ):
+                    full_dims.add(dim)
+
+        out_shape = tuple(len(requested[d]) for d in self._dimension_columns)
+        if any(n == 0 for n in out_shape):
+            empty = np.empty(out_shape, dtype=self.dtype)
+            squeezed = (
+                np.squeeze(empty, axis=tuple(drop_axes)) if drop_axes else empty
+            )
+            return cast(np.ndarray, squeezed)
+
+        # Build a single DataFusion filter expression as the AND of per-dim
+        # predicates. For a single requested value: equality. For multiple:
+        # OR-chain of equalities (DataFusion 52.0.0 does not expose a clean
+        # ``Expr.in_list`` from Python; OR-chained equalities constant-fold
+        # equivalently and stay typed).
+        predicates = []
+        for dim in self._dimension_columns:
+            if dim in full_dims:
+                continue
+            vals = requested[dim]
+            if len(vals) == 1:
+                predicates.append(col(f'"{dim}"') == literal(vals[0]))
+            else:
+                eq = col(f'"{dim}"') == literal(vals[0])
+                for v in vals[1:]:
+                    eq = eq | (col(f'"{dim}"') == literal(v))
+                predicates.append(eq)
+
+        filtered = self._inner_df
+        if predicates:
+            combined = predicates[0]
+            for p in predicates[1:]:
+                combined = combined & p
+            filtered = filtered.filter(combined)
+        projected = filtered.select(
+            *(col(f'"{c}"') for c in self._dimension_columns + [self._var_name])
+        )
+
+        # Consume the projected DataFrame as Arrow RecordBatches. The
+        # DataFusion wrapper exposes ``.to_pyarrow()`` to convert each
+        # batch into a true ``pyarrow.RecordBatch``.
+        batches = [b.to_pyarrow() for b in projected.execute_stream()]
+        return _scatter_batches_to_ndarray(
+            batches=batches,
+            dimension_columns=self._dimension_columns,
+            requested=requested,
+            var_name=self._var_name,
+            out_shape=out_shape,
+            dtype=self.dtype,
+            drop_axes=drop_axes,
+        )
+
+
+def _materialize(
+    inner_df: Any,
+    dimension_columns: list[str],
+    field_names: list[str],
+    field_types: dict[str, Any],
+) -> xr.Dataset:
+    """Execute the query once and build a dense in-memory Dataset.
+
+    Runs the plan exactly once via ``execute_stream()`` -- streaming the result
+    as Arrow ``RecordBatch`` es (``datafusion.RecordBatch.to_pyarrow()``) -- then
+    derives both the coordinates and every data variable from that single pass.
+    This is the eager path, used when no output chunking is requested. It never
+    re-executes, so an aggregation over a remote Zarr scan costs exactly one
+    scan, regardless of how many dimensions or variables the result has.
+    """
+    batches = [b.to_pyarrow() for b in inner_df.execute_stream()]
+
+    coord_arrays: dict[str, np.ndarray] = {}
+    for d in dimension_columns:
+        if not batches:
+            coord_arrays[d] = np.asarray([])
+            continue
+        vals = np.concatenate(
+            [
+                b.column(b.schema.names.index(d)).to_numpy(zero_copy_only=False)
+                for b in batches
+            ]
+        )
+        # Preserve the order coordinate values first appear in the result so an
+        # ORDER BY direction (e.g. ``ORDER BY level DESC``) carries through to
+        # the Dataset dimension instead of being force-sorted ascending.
+        # pd.unique keeps first-appearance order; the scatter below argsorts
+        # internally, so arbitrarily-ordered coordinates are placed correctly.
+        coord_arrays[d] = np.asarray(pd.unique(vals))
+    shape = tuple(len(coord_arrays[d]) for d in dimension_columns)
+
+    data_vars: dict[str, xr.Variable] = {}
+    for name in field_names:
+        if name in dimension_columns:
+            continue
+        np_dtype = np.dtype(field_types[name].to_pandas_dtype())
+        dense = _scatter_batches_to_ndarray(
+            batches=batches,
+            dimension_columns=dimension_columns,
+            requested=coord_arrays,
+            var_name=name,
+            out_shape=shape,
+            dtype=np_dtype,
+            drop_axes=[],
+        )
+        data_vars[name] = xr.Variable(dimension_columns, dense)
+
+    coords_arg = {d: coord_arrays[d] for d in dimension_columns}
+    return xr.Dataset(data_vars=data_vars, coords=coords_arg)
+
+
+def _build_lazy_scan(
+    inner_df: Any,
+    dimension_columns: list[str],
+    field_names: list[str],
+    field_types: dict[str, Any],
+) -> xr.Dataset:
+    """Build a lazy Dataset whose data vars are :class:`SQLBackendArray`.
+
+    Used when output chunking is requested: each data variable stays lazy and,
+    once wrapped by ``Dataset.chunk``, every chunk reads its coordinate range via
+    a pushdown filter on first access. Coordinates are discovered per dim via
+    ``inner_df.select(col(d)).distinct().sort(...).execute_stream()``; the table
+    provider projects to that single coordinate column and skips data variables,
+    so discovery reads coordinate values only (no data-variable I/O).
+    """
+    coord_arrays: dict[str, np.ndarray] = {}
+    for d in dimension_columns:
+        dim_only = (
+            inner_df.select(col(f'"{d}"')).distinct().sort(col(f'"{d}"').sort())
+        )
+        chunks = [b.to_pyarrow() for b in dim_only.execute_stream()]
+        if not chunks:
+            coord_arrays[d] = np.asarray([])
+            continue
+        coord_arrays[d] = np.concatenate(
+            [c.column(0).to_numpy(zero_copy_only=False) for c in chunks]
+        )
+    shape = tuple(len(coord_arrays[d]) for d in dimension_columns)
+
+    data_vars: dict[str, xr.Variable] = {}
+    for name in field_names:
+        if name in dimension_columns:
+            continue
+        np_dtype = field_types[name].to_pandas_dtype()
+        backend = SQLBackendArray(
+            inner_df=inner_df,
+            var_name=name,
+            dimension_columns=dimension_columns,
+            coord_arrays=coord_arrays,
+            shape=shape,
+            dtype=np_dtype,
+        )
+        lazy = xr.core.indexing.LazilyIndexedArray(backend)
+        data_vars[name] = xr.Variable(dimension_columns, lazy)
+
+    coords_arg = {d: coord_arrays[d] for d in dimension_columns}
+    return xr.Dataset(data_vars=data_vars, coords=coords_arg)
+
+
+def _auto_chunk_target_bytes() -> int:
+    """Byte target for ``chunks="auto"`` (the chunk manager's, else 128 MiB)."""
+    try:
+        import dask
+        from dask.utils import parse_bytes
+
+        return int(parse_bytes(dask.config.get("array.chunk-size")))
+    except Exception:
+        return 128 * 1024 * 1024
+
+
+def _auto_chunks(
+    template: xr.Dataset | None,
+    dimension_columns: list[str],
+    field_types: dict[str, Any],
+) -> dict[str, int] | None:
+    """Resolve ``chunks="auto"`` to a source-partition-aligned chunk spec.
+
+    Sizes chunks to roughly the chunk manager's byte target (dask's
+    ``array.chunk-size``, default 128 MiB) but snaps boundaries to whole source
+    partitions, so every chunk is a union of source partitions -- no chunk splits
+    a partition (which would make adjacent chunks re-read it). This is what makes
+    ``"auto"`` useful for finely partitioned sources (e.g. ERA5
+    ``chunks={"time": 1}``): it coarsens many tiny partitions into memory-sized,
+    aligned chunks. Returns ``None`` when there is no resolvable source grid to
+    align to, so the caller falls back to the chunk manager's own ``"auto"``.
+    """
+    if template is None:
+        return None
+    part = template.chunksizes  # dim -> tuple of source chunk lengths
+    chunked_dims = [
+        d for d in dimension_columns if d in part and len(part[d]) > 1
+    ]
+    if not chunked_dims:
+        return None
+
+    itemsizes = [
+        np.dtype(t.to_pandas_dtype()).itemsize
+        for name, t in field_types.items()
+        if name not in dimension_columns
+    ]
+    itemsize = max(itemsizes) if itemsizes else 8
+
+    # Bytes in one source-partition block: the nominal source chunk length per
+    # dimension (``part[d][0]``) multiplied across all dims, times itemsize.
+    block_bytes = itemsize
+    for d in dimension_columns:
+        if d in part:
+            block_bytes *= int(part[d][0])
+    # Number of source partitions to merge per chunk to approach the target.
+    merge = max(1, _auto_chunk_target_bytes() // max(block_bytes, 1))
+
+    # Absorb the coarsening into the most finely partitioned dimension; the rest
+    # keep their source chunk length. xarray caps an oversize chunk at the dim
+    # length, so an over-large merge simply yields a single chunk on that dim.
+    primary = max(chunked_dims, key=lambda d: len(part[d]))
+    return {
+        d: int(part[d][0]) * (merge if d == primary else 1)
+        for d in chunked_dims
+    }
+
+
+def _result_to_xarray(
+    inner_df: Any,
+    dimension_columns: list[str],
+    template: xr.Dataset | None,
+    sparsity: Sparsity,
+    fill_value: Any,
+    chunks: Mapping[str, int] | str | None,
+) -> xr.Dataset:
+    """Reconstruct an ``xr.Dataset`` from a SQL result.
+
+    ``chunks`` (already resolved by :meth:`XarrayDataFrame._resolve_chunks`)
+    selects the execution strategy:
+
+    * ``None`` -> eager: execute once and materialize a dense Dataset
+      (:func:`_materialize`). Correct for any query and the right default for
+      reductions, whose results are small.
+    * a mapping (or ``"auto"``) -> lazy/chunked: build :class:`SQLBackendArray`
+      data variables (:func:`_build_lazy_scan`) and wrap them with
+      ``Dataset.chunk`` so each chunk reads its coordinate range via filter
+      pushdown. The chunk grid maps onto the source partitions. Chunking goes
+      through xarray's configured chunk manager (dask, cubed, ...), so no
+      chunked-array backend is imported directly here.
+    """
+    if sparsity not in ("result", "template"):
+        raise ValueError(
+            f"sparsity must be 'result' or 'template', got {sparsity!r}"
+        )
+    if sparsity == "template" and template is None:
+        raise ValueError(
+            "sparsity='template' requires template= to be supplied"
+        )
+
+    schema = inner_df.schema()
+    field_names = [f.name for f in schema]
+    field_types = {f.name: f.type for f in schema}
+
+    if chunks is None:
+        ds = _materialize(inner_df, dimension_columns, field_names, field_types)
+    else:
+        ds = _build_lazy_scan(
+            inner_df, dimension_columns, field_names, field_types
+        )
+
+    if sparsity == "template":
+        assert template is not None
+        indexers = {
+            d: template.coords[d].values
+            for d in dimension_columns
+            if d in template.coords and d in template.dims
+        }
+        if indexers:
+            ds = ds.reindex(indexers, fill_value=fill_value)
+
+    if template is not None:
+        ds = _apply_template(ds, template)
+
+    if chunks is not None:
+        if chunks == "auto":
+            # Snap the byte-budgeted "auto" sizing to source partition
+            # boundaries; fall back to the chunk manager's own "auto" when there
+            # is no source grid to align to.
+            chunks = (
+                _auto_chunks(template, dimension_columns, field_types) or "auto"
+            )
+        # Wrap the lazy data variables in the configured chunk manager (dask by
+        # default). Each chunk reads its coordinate range via pushdown on access.
+        ds = ds.chunk(chunks)
+    return ds
+
+
+# ---------------------------------------------------------------------------
+# Public wrapper
+# ---------------------------------------------------------------------------
+
+
+class XarrayDataFrame:
+    """Wrapper around a DataFusion ``DataFrame`` with xarray-aware helpers.
+
+    Returned by :meth:`xarray_sql.XarrayContext.sql`. Forwards every
+    attribute it does not define itself to the wrapped DataFrame, so
+    ``.collect()``, ``.schema()``, ``.show()``, ``.count()`` all work
+    unchanged.
+
+    Carries a private snapshot of the context's registered Datasets so
+    :meth:`to_dataset` can default ``dims`` and recover metadata
+    dropped by the forward pivot.
+
+    Users should not construct this class directly; let
+    :meth:`XarrayContext.sql` produce it.
+    """
+
+    def __init__(
+        self,
+        inner: Any,
+        templates: dict[str, xr.Dataset] | None = None,
+    ) -> None:
+        """Construct a wrapper.
+
+        Args:
+            inner: The underlying ``datafusion.DataFrame`` returned by
+                :meth:`XarrayContext.sql`.
+            templates: Snapshot of the registered Datasets on the producing
+                context, keyed by the SQL identifier each was registered
+                under. Used by :meth:`to_dataset` to recover metadata that
+                the forward pivot strips. ``None`` means no metadata
+                recovery is possible from registrations alone; callers may
+                still pass ``template=`` to :meth:`to_dataset` explicitly.
+        """
+        object.__setattr__(self, "_inner", inner)
+        object.__setattr__(self, "_templates", dict(templates or {}))
+
+    def to_pandas(self) -> pd.DataFrame:
+        """Materialize the result as a ``pd.DataFrame`` (DataFusion API)."""
+        return self._inner.to_pandas()
+
+    def to_dataset(
+        self,
+        dims: list[str] | None = None,
+        template: xr.Dataset | str | None = None,
+        sparsity: Sparsity = "result",
+        fill_value: Any = np.nan,
+        chunks: Mapping[str, int] | str | None = "inherit",
+    ) -> xr.Dataset:
+        """Convert the result to an ``xr.Dataset``.
+
+        Args:
+            dims: Result columns to use as Dataset dimensions. When
+                ``None``, defaults to the dims of the registered Dataset
+                referenced by the SQL ``FROM`` clause (if exactly one
+                matches), or any single registered Dataset whose dims are
+                all present in the result columns.
+            template: Source to recover metadata (attrs, encoding, non-dim
+                coordinates, dim-coord dtype) from. Either an ``xr.Dataset``
+                used directly, or the name of a registered table (e.g.
+                ``"era5.surface"``) whose Dataset is looked up. When ``None``
+                and exactly one Dataset is registered, that one is used.
+            sparsity: ``"result"`` (default) keeps only dim values
+                present in the result. ``"template"`` reindexes to the
+                template's full coord ranges, filling absent cells with
+                ``fill_value``; requires a template.
+            fill_value: Used when ``sparsity="template"`` reindexes
+                to a wider extent. Defaults to ``np.nan``.
+            chunks: Output chunking, controlling laziness (an xarray idiom).
+
+                * ``"inherit"`` (default): reuse the source Dataset's chunk
+                  sizes, but only for dimensions that were genuinely split into
+                  multiple chunks in the input -- so the output chunk grid maps
+                  onto the source partitions. A reduction that drops the chunked
+                  dimension (e.g. a global aggregation) inherits nothing and so
+                  is materialized eagerly. Falls back to eager when no source
+                  Dataset is resolvable.
+                * ``None``: eager. Execute the query once and return a dense
+                  in-memory Dataset. Best for reductions (small results).
+                * a mapping (e.g. ``{"time": 100}``): chunk explicitly. Each
+                  chunk reads its coordinate range lazily via filter pushdown on
+                  access, through xarray's configured chunk manager (dask,
+                  cubed, ...).
+                * ``"auto"``: size chunks to the chunk manager's byte target but
+                  snap boundaries to whole source partitions, so each chunk is a
+                  union of source partitions. Useful for finely partitioned
+                  sources (e.g. ERA5 ``chunks={"time": 1}``), coarsening many
+                  tiny partitions into memory-sized, aligned chunks.
+
+        Returns:
+            An ``xr.Dataset`` with ``dims`` as dimensions and the
+            remaining result columns as data variables.
+
+        Raises:
+            ValueError: ``dims`` cannot be inferred, names a missing
+                column, or the result has duplicate dim tuples;
+                ``template`` names an unknown registered table; or
+                ``sparsity="template"`` is requested without a
+                resolvable template.
+        """
+        if not isinstance(template, xr.Dataset):
+            # ``template`` is a registered-table name or None; look it up.
+            template = self._resolve_template(template)
+        if dims is None:
+            dims = self._infer_dimension_columns(preferred_template=template)
+        resolved_chunks = self._resolve_chunks(chunks, template, dims)
+        return _result_to_xarray(
+            inner_df=self._inner,
+            dimension_columns=dims,
+            template=template,
+            sparsity=sparsity,
+            fill_value=fill_value,
+            chunks=resolved_chunks,
+        )
+
+    # ------------------------------------------------------------------
+    # Internals
+    # ------------------------------------------------------------------
+
+    @staticmethod
+    def _resolve_chunks(
+        chunks: Mapping[str, int] | str | None,
+        template: xr.Dataset | None,
+        dimension_columns: list[str],
+    ) -> Mapping[str, int] | str | None:
+        """Resolve the ``chunks`` argument to a concrete spec or ``None``.
+
+        ``None`` selects the eager path; anything else selects the lazy/chunked
+        path. ``"inherit"`` reuses the source Dataset's chunk sizes -- but only
+        for dimensions actually split into more than one chunk in the input
+        (a single full chunk is not "chunked"), so reductions that drop the
+        chunked dimension resolve to ``None`` (eager) automatically. Mappings
+        pass through unchanged; ``"auto"`` passes through here and is snapped to
+        source partition boundaries later (see :func:`_auto_chunks`).
+        """
+        if chunks is None:
+            return None
+        if chunks == "inherit":
+            if template is None:
+                return None
+            sizes = template.chunksizes  # dim -> tuple of chunk lengths
+            inherited = {
+                d: sizes[d][0]
+                for d in dimension_columns
+                if d in sizes and len(sizes[d]) > 1
+            }
+            return inherited or None
+        return chunks
+
+    def _resolve_template(self, name: str | None) -> xr.Dataset | None:
+        """Pick a template Dataset for metadata recovery by registered name.
+
+        Priority:
+          1. The named registered table (``name``).
+          2. If exactly one Dataset is registered on the context, use it.
+          3. None.
+        """
+        templates = self._templates
+        if name is not None:
+            if name not in templates:
+                raise ValueError(
+                    f"template={name!r} is not a registered table on this "
+                    f"context. Registered: {list(templates)}"
+                )
+            return templates[name]
+        if len(templates) == 1:
+            return next(iter(templates.values()))
+        return None
+
+    def _infer_dimension_columns(
+        self, preferred_template: xr.Dataset | None = None
+    ) -> list[str]:
+        """Pick a default ``dimension_columns`` from the registry, or raise.
+
+        Uses the data variable's dim order (via :func:`_ds_var_dims`) so
+        the round-trip preserves the original axis order.
+        """
+        result_cols = set(self._result_columns())
+        if (
+            preferred_template is not None
+            and set(preferred_template.dims) <= result_cols
+        ):
+            return _ds_var_dims(preferred_template)
+        if not self._templates:
+            raise ValueError(
+                "dims cannot be inferred (no registered "
+                "Dataset on this result); pass dims=[...] "
+                "explicitly."
+            )
+        candidates = [
+            _ds_var_dims(t)
+            for t in self._templates.values()
+            if set(t.dims) <= result_cols
+        ]
+        if len(candidates) == 1:
+            return candidates[0]
+        if not candidates:
+            raise ValueError(
+                "dims cannot be inferred: no registered "
+                "Dataset has all of its dims present in the result "
+                "columns. Pass dims=[...] explicitly."
+            )
+        raise ValueError(
+            "dims cannot be inferred unambiguously: multiple "
+            "registered Datasets are compatible with the result. Pass "
+            "dims=[...] explicitly."
+        )
+
+    def _result_columns(self) -> list[str]:
+        """Return the result's column names without materializing rows."""
+        return [field.name for field in self._inner.schema()]
+
+    def __getattr__(self, name: str) -> Any:
+        # Runs only when ``name`` is not found via normal lookup, so this
+        # safely forwards anything we have not overridden.
+        return getattr(self._inner, name)
+
+    def __repr__(self) -> str:
+        return repr(self._inner)
diff --git a/xarray_sql/sql.py b/xarray_sql/sql.py
index 403eaee..99efea0 100644
--- a/xarray_sql/sql.py
+++ b/xarray_sql/sql.py
@@ -5,12 +5,23 @@
 
 from . import cftime as cft
 from .df import Chunks
+from .ds import XarrayDataFrame
 from .reader import read_xarray_table
 
 
 class XarrayContext(SessionContext):
     """A datafusion `SessionContext` that also supports `xarray.Dataset`s."""
 
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        # Track registered xarray Datasets so XarrayDataFrame can recover
+        # defaults (dimension_columns) and metadata (var/dataset attrs,
+        # non-dim coords, dim-coord dtype) that the forward pivot drops.
+        # Keys are the fully-qualified table names users will reference
+        # in SQL (e.g. ``"air"`` for a uniform-dim Dataset, or
+        # ``"era5.surface"`` for one entry from a multi-dim-group split).
+        self._registered_datasets: dict[str, xr.Dataset] = {}
+
     def from_dataset(
         self,
         name: str,
@@ -79,6 +90,7 @@ def from_dataset(
         groups = _group_vars_by_dims(input_table)
 
         if len(groups) <= 1:
+            self._registered_datasets[name] = input_table
             return self._from_dataset(name, input_table, chunks)
 
         table_names = table_names or {}
@@ -89,9 +101,11 @@ def from_dataset(
             # Scalar variables group under empty dims, where "_".join(()) is
             # the empty string; fall back to a valid default table name.
             sub_name = table_names.get(dims, "_".join(dims) or "scalar")
-            self._from_dataset(
-                sub_name, input_table[var_names], chunks, schema=schema
-            )
+            sub_ds = input_table[var_names]
+            self._from_dataset(sub_name, sub_ds, chunks, schema=schema)
+            # Track the fully-qualified name so XarrayDataFrame metadata
+            # recovery can find this Dataset on round-trip.
+            self._registered_datasets[f"{name}.{sub_name}"] = sub_ds
 
         return self
 
@@ -123,6 +137,27 @@ def _maybe_register_cftime_udf(self, ds: xr.Dataset) -> None:
                     self.register_udf(cft.make_cftime_udf(units, cal))
                     break  # One UDF per context is enough.
 
+    def sql(self, query: str, *args, **kwargs) -> XarrayDataFrame:
+        """Run a SQL query, returning an :class:`XarrayDataFrame` wrapper.
+
+        Identical to ``datafusion.SessionContext.sql`` except the returned
+        object wraps the DataFusion DataFrame. The wrapper exposes
+        ``.to_pandas()`` (unchanged), forwards every other DataFusion
+        method via ``__getattr__``, and adds
+        ``.to_dataset(dimension_columns=[...])`` for round-tripping the
+        result back to an ``xr.Dataset``.
+
+        Args:
+            query: A SQL query string.
+            *args: Forwarded to ``SessionContext.sql``.
+            **kwargs: Forwarded to ``SessionContext.sql``.
+
+        Returns:
+            An :class:`XarrayDataFrame` wrapping the DataFusion DataFrame.
+        """
+        inner = super().sql(query, *args, **kwargs)
+        return XarrayDataFrame(inner, templates=self._registered_datasets)
+
 
 def _group_vars_by_dims(ds: xr.Dataset) -> dict[tuple[str, ...], list[str]]:
     """Group variables in the dataset based on shared dims.