From 094fbbd89618066560050a98dd780f5bb723cbb3 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Wed, 24 Jun 2026 14:24:43 -0700
Subject: [PATCH 1/3] Add .xrs.validate() contract compliance check (#3485)

---
 xrspatial/__init__.py                     |   1 +
 xrspatial/accessor.py                     |  46 ++
 xrspatial/tests/test_accessor.py          |   2 +
 xrspatial/tests/test_contract_validate.py | 314 +++++++++++++
 xrspatial/validate.py                     | 533 ++++++++++++++++++++++
 5 files changed, 896 insertions(+)
 create mode 100644 xrspatial/tests/test_contract_validate.py
 create mode 100644 xrspatial/validate.py

diff --git a/xrspatial/__init__.py b/xrspatial/__init__.py
index 506f0bb1b..3cf07911c 100644
--- a/xrspatial/__init__.py
+++ b/xrspatial/__init__.py
@@ -119,6 +119,7 @@
 from xrspatial.terrain_metrics import roughness  # noqa
 from xrspatial.terrain_metrics import tpi  # noqa
 from xrspatial.terrain_metrics import tri  # noqa
+from xrspatial.validate import validate  # noqa
 from xrspatial.hydro import twi  # noqa: unified wrapper
 from xrspatial.hydro import twi_d8  # noqa
 from xrspatial.polygon_clip import clip_polygon  # noqa
diff --git a/xrspatial/accessor.py b/xrspatial/accessor.py
index c7013e9cd..d3fe692fa 100644
--- a/xrspatial/accessor.py
+++ b/xrspatial/accessor.py
@@ -1505,6 +1505,30 @@ def multi_overlap(self, func, n_outputs, **kwargs):
         from .utils import multi_overlap
         return multi_overlap(self._obj, func, n_outputs, **kwargs)
 
+    # ---- Diagnostics ----
+
+    def validate(self, *, raise_on_error=False):
+        """Check this DataArray against the xarray-spatial input contract.
+
+        Returns a :class:`~xrspatial.validate.ValidationReport` listing
+        every contract violation as an error (a spatial op will fail) or
+        a warning (behavior degrades), each with a suggested fix. The
+        report is truthy when there are no error-level issues, so
+        ``if not da.xrs.validate(): ...`` reads naturally.
+
+        Parameters
+        ----------
+        raise_on_error : bool, default False
+            If True, raise :class:`~xrspatial.validate.XrsContractError`
+            when any error-level issue is found instead of only
+            recording it in the report.
+        """
+        from .validate import validate
+        report = validate(self._obj)
+        if raise_on_error:
+            report.raise_if_errors()
+        return report
+
 
 @xr.register_dataset_accessor("xrs")
 class XrsSpatialDatasetAccessor:
@@ -2116,6 +2140,28 @@ def rechunk_no_shuffle(self, **kwargs):
         from .utils import rechunk_no_shuffle
         return rechunk_no_shuffle(self._obj, **kwargs)
 
+    # ---- Diagnostics ----
+
+    def validate(self, *, raise_on_error=False):
+        """Check each data variable against the xarray-spatial contract.
+
+        Returns a
+        :class:`~xrspatial.validate.DatasetValidationReport` holding a
+        per-variable :class:`~xrspatial.validate.ValidationReport`. The
+        report is truthy when every variable is compliant.
+
+        Parameters
+        ----------
+        raise_on_error : bool, default False
+            If True, raise :class:`~xrspatial.validate.XrsContractError`
+            when any variable has an error-level issue.
+        """
+        from .validate import validate_dataset
+        report = validate_dataset(self._obj)
+        if raise_on_error:
+            report.raise_if_errors()
+        return report
+
 
 # ---------------------------------------------------------------------------
 # Surface standalone-function docstrings on accessor methods so that, e.g.,
diff --git a/xrspatial/tests/test_accessor.py b/xrspatial/tests/test_accessor.py
index 0c58db010..329ded8b4 100644
--- a/xrspatial/tests/test_accessor.py
+++ b/xrspatial/tests/test_accessor.py
@@ -100,6 +100,7 @@ def test_dataarray_accessor_has_expected_methods(elevation):
         'rechunk_no_shuffle',
         'fused_overlap',
         'multi_overlap',
+        'validate',
     ]
     for name in expected:
         assert name in names, f"Missing method: {name}"
@@ -118,6 +119,7 @@ def test_dataset_accessor_has_expected_methods():
         'proximity', 'allocation', 'direction', 'cost_distance',
         'ndvi', 'evi', 'arvi', 'savi', 'nbr', 'sipi',
         'rasterize',
+        'validate',
     ]
     for name in expected:
         assert name in names, f"Missing method: {name}"
diff --git a/xrspatial/tests/test_contract_validate.py b/xrspatial/tests/test_contract_validate.py
new file mode 100644
index 000000000..962b94f89
--- /dev/null
+++ b/xrspatial/tests/test_contract_validate.py
@@ -0,0 +1,314 @@
+"""Tests for the xarray-spatial input-contract check (`xrs.validate`)."""
+import numpy as np
+import pytest
+import xarray as xr
+
+import xrspatial  # noqa: F401  (registers the .xrs accessor)
+from xrspatial.validate import (
+    DatasetValidationReport,
+    ValidationReport,
+    XrsContractError,
+    validate,
+    validate_dataset,
+)
+
+from .general_checks import (
+    create_test_raster,
+    cuda_and_cupy_available,
+    dask_array_available,
+)
+
+
+def _compliant(name="dem"):
+    """A raster that satisfies the whole contract."""
+    data = np.ones((4, 4), dtype="float64")
+    return xr.DataArray(
+        data,
+        dims=["y", "x"],
+        coords={"x": np.arange(4.0), "y": np.arange(4.0)},
+        attrs={"crs": "EPSG:4326"},
+        name=name,
+    )
+
+
+def _checks(report):
+    return {i.check for i in report.issues}
+
+
+# ---------------------------------------------------------------------------
+# Compliant input
+# ---------------------------------------------------------------------------
+
+def test_compliant_raster_is_valid():
+    report = _compliant().xrs.validate()
+    assert isinstance(report, ValidationReport)
+    assert report.is_valid is True
+    assert bool(report) is True
+    assert report.issues == []
+
+
+# ---------------------------------------------------------------------------
+# Error-level checks
+# ---------------------------------------------------------------------------
+
+def test_non_dataarray_is_error():
+    report = validate([1, 2, 3])
+    assert not report
+    assert _checks(report) == {"type"}
+
+
+@pytest.mark.parametrize("shape", [(5,), (2, 2, 2, 2)])
+def test_wrong_ndim_is_error(shape):
+    da = xr.DataArray(np.zeros(shape, dtype="float64"))
+    report = da.xrs.validate()
+    assert not report
+    assert "ndim" in _checks(report)
+
+
+def test_string_dtype_is_error():
+    da = _compliant()
+    da = da.astype("<U5")
+    report = da.xrs.validate()
+    assert not report
+    assert "dtype" in _checks(report)
+
+
+def test_complex_dtype_is_error():
+    da = _compliant().astype("complex128")
+    report = da.xrs.validate()
+    assert not report
+    assert "dtype" in _checks(report)
+
+
+def test_non_numeric_coords_is_error():
+    da = _compliant()
+    da = da.assign_coords(x=np.array(["a", "b", "c", "d"]))
+    report = da.xrs.validate()
+    assert not report
+    assert "coords_numeric" in _checks(report)
+
+
+def test_non_monotonic_coords_is_error():
+    da = _compliant()
+    da = da.assign_coords(x=np.array([0.0, 2.0, 1.0, 3.0]))
+    report = da.xrs.validate()
+    assert not report
+    assert "monotonic" in _checks(report)
+
+
+def test_non_finite_cellsize_is_error():
+    # A non-finite coordinate makes the inferred cell size non-finite,
+    # which planar ops cannot use.
+    da = _compliant()
+    da = da.assign_coords(y=np.array([0.0, 1.0, 2.0, np.inf]))
+    report = da.xrs.validate()
+    assert not report
+    assert "cellsize" in _checks(report)
+
+
+# ---------------------------------------------------------------------------
+# Warning-level checks
+# ---------------------------------------------------------------------------
+
+def test_uneven_spacing_is_warning():
+    da = _compliant()
+    da = da.assign_coords(x=np.array([0.0, 1.0, 2.0, 5.0]))
+    report = da.xrs.validate()
+    assert report.is_valid is True  # warning only
+    assert "even_spacing" in _checks(report)
+
+
+def test_missing_coords_is_warning():
+    da = xr.DataArray(
+        np.zeros((4, 4), dtype="float64"), dims=["y", "x"],
+        attrs={"crs": "EPSG:4326"},
+    )
+    report = da.xrs.validate()
+    assert report.is_valid is True
+    assert "coords_present" in _checks(report)
+
+
+def test_unconventional_dim_names_is_warning():
+    da = xr.DataArray(
+        np.zeros((4, 4), dtype="float64"), dims=["row", "col"],
+        coords={"col": np.arange(4.0), "row": np.arange(4.0)},
+        attrs={"crs": "EPSG:4326"},
+    )
+    report = da.xrs.validate()
+    assert report.is_valid is True
+    assert "spatial_dims" in _checks(report)
+
+
+def test_projected_latlon_is_warning():
+    da = xr.DataArray(
+        np.zeros((3, 3), dtype="float64"), dims=["lat", "lon"],
+        coords={"lat": [0.0, 1.0, 2.0], "lon": [500.0, 501.0, 502.0]},
+        attrs={"crs": "EPSG:4326"},
+    )
+    report = da.xrs.validate()
+    assert report.is_valid is True
+    assert "geographic_range" in _checks(report)
+
+
+def test_missing_crs_is_warning():
+    da = _compliant()
+    da.attrs = {}
+    report = da.xrs.validate()
+    assert report.is_valid is True
+    assert "crs" in _checks(report)
+
+
+# ---------------------------------------------------------------------------
+# Severity semantics
+# ---------------------------------------------------------------------------
+
+def test_warnings_only_is_valid():
+    # missing coords + missing crs + odd dim names: all warnings.
+    da = xr.DataArray(np.zeros((4, 4), dtype="float64"), dims=["row", "col"])
+    report = da.xrs.validate()
+    assert report.warnings
+    assert not report.errors
+    assert report.is_valid is True
+
+
+def test_one_error_flips_validity():
+    da = xr.DataArray(np.zeros((4, 4), dtype="<U3"), dims=["row", "col"])
+    report = da.xrs.validate()
+    assert report.errors
+    assert report.is_valid is False
+
+
+# ---------------------------------------------------------------------------
+# raise_on_error
+# ---------------------------------------------------------------------------
+
+def test_raise_on_error_raises_and_names_every_error():
+    da = _compliant().astype("<U5").assign_coords(
+        x=np.array(["a", "b", "c", "d"])
+    )
+    report = da.xrs.validate()
+    assert {"dtype", "coords_numeric"} <= _checks(report)
+    with pytest.raises(XrsContractError) as exc:
+        da.xrs.validate(raise_on_error=True)
+    msg = str(exc.value)
+    # the aggregated message names both error-level issues by their text
+    assert "dtype is" in msg
+    assert "not numeric" in msg
+
+
+def test_raise_on_error_no_raise_for_warnings_only():
+    da = _compliant()
+    da.attrs = {}  # only a crs warning
+    # should not raise
+    report = da.xrs.validate(raise_on_error=True)
+    assert report.is_valid is True
+
+
+# ---------------------------------------------------------------------------
+# Backends — structural check is identical and never materializes data
+# ---------------------------------------------------------------------------
+
+def _backend_raster(backend, dtype="float64"):
+    data = np.ones((6, 6), dtype=dtype)
+    # create_test_raster supplies res + crs attrs and y/x coords.
+    return create_test_raster(data, backend=backend, name="dem")
+
+
+@pytest.mark.parametrize("backend", ["numpy"])
+def test_numpy_backend_compliant(backend):
+    assert _backend_raster(backend).xrs.validate().is_valid
+
+
+@dask_array_available
+def test_dask_backend_stays_lazy():
+    import dask.array as da
+    raster = _backend_raster("dask+numpy")
+    report = raster.xrs.validate()
+    assert report.is_valid
+    # validate must not have computed the dask graph.
+    assert isinstance(raster.data, da.Array)
+
+
+@dask_array_available
+def test_dask_backend_reports_same_as_numpy():
+    np_raster = _backend_raster("numpy", dtype="<U3")
+    dask_raster = _backend_raster("dask+numpy", dtype="<U3")
+    assert _checks(np_raster.xrs.validate()) == _checks(
+        dask_raster.xrs.validate()
+    )
+
+
+@cuda_and_cupy_available
+def test_cupy_backend_stays_on_device():
+    import cupy
+    raster = _backend_raster("cupy")
+    report = raster.xrs.validate()
+    assert report.is_valid
+    assert isinstance(raster.data, cupy.ndarray)
+
+
+@cuda_and_cupy_available
+def test_dask_cupy_backend_compliant():
+    raster = _backend_raster("dask+cupy")
+    assert raster.xrs.validate().is_valid
+
+
+# ---------------------------------------------------------------------------
+# Dataset accessor
+# ---------------------------------------------------------------------------
+
+def test_dataset_per_variable_report():
+    good = _compliant("dem")
+    bad = good.astype("<U3").rename("labels")
+    ds = xr.Dataset({"dem": good, "labels": bad})
+    report = ds.xrs.validate()
+    assert isinstance(report, DatasetValidationReport)
+    assert report.is_valid is False
+    assert report.reports["dem"].is_valid is True
+    assert report.reports["labels"].is_valid is False
+    assert "labels" in repr(report)
+
+
+def test_dataset_all_compliant_is_valid():
+    ds = xr.Dataset({"a": _compliant("a"), "b": _compliant("b")})
+    assert ds.xrs.validate().is_valid
+
+
+def test_dataset_raise_on_error():
+    ds = xr.Dataset(
+        {"dem": _compliant("dem"), "labels": _compliant("labels").astype("<U3")}
+    )
+    with pytest.raises(XrsContractError) as exc:
+        ds.xrs.validate(raise_on_error=True)
+    assert "labels" in str(exc.value)
+
+
+# ---------------------------------------------------------------------------
+# Repr
+# ---------------------------------------------------------------------------
+
+def test_text_repr_lists_messages_and_suggestions():
+    da = _compliant().astype("<U5")
+    report = da.xrs.validate()
+    text = repr(report)
+    assert "error" in text
+    assert "dtype" in text
+    assert "Fix:" in text
+
+
+def test_html_repr_renders_table():
+    da = _compliant().astype("<U5")
+    html = da.xrs.validate()._repr_html_()
+    assert "<table" in html
+    assert "ValidationReport" in html
+
+
+def test_compliant_repr_says_compliant():
+    assert "compliant" in repr(_compliant().xrs.validate())
+
+
+def test_validate_function_matches_accessor():
+    da = _compliant().astype("<U5")
+    assert _checks(validate(da)) == _checks(da.xrs.validate())
+    ds = xr.Dataset({"dem": _compliant()})
+    assert validate_dataset(ds).is_valid == ds.xrs.validate().is_valid
diff --git a/xrspatial/validate.py b/xrspatial/validate.py
new file mode 100644
index 000000000..bda1bfcd1
--- /dev/null
+++ b/xrspatial/validate.py
@@ -0,0 +1,533 @@
+"""Check a raster against the xarray-spatial input contract.
+
+xrspatial operations share an implicit contract about their input
+DataArray: it must be an ``xr.DataArray`` whose last two dims are the
+spatial axes (2-D, or 3-D with a leading band/time axis), with a real
+numeric dtype, numeric/monotonic/evenly-spaced coordinates, a positive
+finite cell size, lat/lon in geographic range, and a CRS for geodesic
+methods and GeoTIFF I/O. Those rules live in the per-op validators in
+``utils.py`` and surface one at a time when an op raises mid-pipeline.
+
+``validate`` checks all of them up front and returns a report of every
+violation with a suggested fix, rather than raising on the first one.
+Every check inspects structure only (dims, coords, dtype, attrs) and
+never reads the data buffer, so dask and cupy arrays validate without
+materializing and all four backends behave identically.
+"""
+
+from __future__ import annotations
+
+import html
+from dataclasses import dataclass
+
+import numpy as np
+import xarray as xr
+
+__all__ = [
+    "validate",
+    "validate_dataset",
+    "ValidationIssue",
+    "ValidationReport",
+    "DatasetValidationReport",
+    "XrsContractError",
+]
+
+# Names that mark a dimension as a recognised spatial axis. Mirrors the
+# lat/lon name sets used by the geodesic helpers in utils.py.
+_Y_DIM_NAMES = {"y", "lat", "latitude"}
+_X_DIM_NAMES = {"x", "lon", "longitude"}
+_LAT_NAMES = {"lat", "latitude"}
+_LON_NAMES = {"lon", "longitude"}
+
+
+class XrsContractError(ValueError):
+    """Raised by ``validate(raise_on_error=True)`` when errors are found."""
+
+
+@dataclass(frozen=True)
+class ValidationIssue:
+    """One way an array is out of compliance with the contract.
+
+    Parameters
+    ----------
+    severity : str
+        ``'error'`` (a spatial op will fail) or ``'warning'`` (behavior
+        degrades, e.g. distance-based results may be inaccurate).
+    check : str
+        Short identifier for the failed check, e.g. ``'dtype'``.
+    message : str
+        What is wrong.
+    suggestion : str
+        How to get back into compliance.
+    """
+
+    severity: str
+    check: str
+    message: str
+    suggestion: str
+
+    def __str__(self) -> str:
+        return f"[{self.severity}] {self.message} Fix: {self.suggestion}"
+
+
+class ValidationReport:
+    """Result of validating a single DataArray against the contract.
+
+    Truthy when the array has no error-level issues (warnings are
+    allowed), so ``if not da.xrs.validate(): ...`` reads naturally.
+    """
+
+    def __init__(self, issues, name=None):
+        self.issues = list(issues)
+        self.name = name
+
+    @property
+    def errors(self):
+        return [i for i in self.issues if i.severity == "error"]
+
+    @property
+    def warnings(self):
+        return [i for i in self.issues if i.severity == "warning"]
+
+    @property
+    def is_valid(self) -> bool:
+        """True when there are no error-level issues."""
+        return not self.errors
+
+    def __bool__(self) -> bool:
+        return self.is_valid
+
+    def raise_if_errors(self):
+        """Raise :class:`XrsContractError` listing every error-level issue."""
+        if not self.errors:
+            return
+        label = f" for {self.name!r}" if self.name else ""
+        lines = [f"raster{label} is not contract-compliant:"]
+        for issue in self.errors:
+            lines.append(f"  - {issue.message} Fix: {issue.suggestion}")
+        raise XrsContractError("\n".join(lines))
+
+    def __repr__(self) -> str:
+        label = f" for {self.name!r}" if self.name else ""
+        if not self.issues:
+            return f"ValidationReport{label}: compliant (0 issues)"
+        n_err = len(self.errors)
+        n_warn = len(self.warnings)
+        header = (
+            f"ValidationReport{label}: "
+            f"{n_err} error(s), {n_warn} warning(s)"
+        )
+        lines = [header]
+        for issue in self.issues:
+            lines.append(
+                f"  [{issue.severity}] {issue.message} "
+                f"Fix: {issue.suggestion}"
+            )
+        return "\n".join(lines)
+
+    def _repr_html_(self) -> str:
+        label = (
+            f" for <code>{html.escape(str(self.name))}</code>"
+            if self.name
+            else ""
+        )
+        if not self.issues:
+            return (
+                f"<div><b>ValidationReport</b>{label}: "
+                "compliant (0 issues)</div>"
+            )
+        rows = []
+        for issue in self.issues:
+            rows.append(
+                "<tr>"
+                f"<td><b>{html.escape(issue.severity)}</b></td>"
+                f"<td>{html.escape(issue.check)}</td>"
+                f"<td>{html.escape(issue.message)}</td>"
+                f"<td>{html.escape(issue.suggestion)}</td>"
+                "</tr>"
+            )
+        n_err = len(self.errors)
+        n_warn = len(self.warnings)
+        return (
+            f"<div><b>ValidationReport</b>{label}: "
+            f"{n_err} error(s), {n_warn} warning(s)"
+            "<table><tr><th>severity</th><th>check</th>"
+            "<th>issue</th><th>suggestion</th></tr>"
+            + "".join(rows)
+            + "</table></div>"
+        )
+
+
+class DatasetValidationReport:
+    """Per-variable validation result for a Dataset.
+
+    Holds an ordered ``{var_name: ValidationReport}`` map. Truthy when
+    every variable is compliant.
+    """
+
+    def __init__(self, reports):
+        # reports: dict mapping variable name -> ValidationReport
+        self.reports = dict(reports)
+
+    @property
+    def is_valid(self) -> bool:
+        return all(r.is_valid for r in self.reports.values())
+
+    def __bool__(self) -> bool:
+        return self.is_valid
+
+    def raise_if_errors(self):
+        """Raise if any variable has error-level issues."""
+        offending = {
+            name: r for name, r in self.reports.items() if not r.is_valid
+        }
+        if not offending:
+            return
+        lines = ["Dataset is not contract-compliant:"]
+        for name, r in offending.items():
+            for issue in r.errors:
+                lines.append(
+                    f"  - [{name}] {issue.message} Fix: {issue.suggestion}"
+                )
+        raise XrsContractError("\n".join(lines))
+
+    def __repr__(self) -> str:
+        if not self.reports:
+            return "DatasetValidationReport: no data variables"
+        status = "compliant" if self.is_valid else "non-compliant"
+        lines = [f"DatasetValidationReport: {status}"]
+        for name, r in self.reports.items():
+            if r.is_valid and not r.warnings:
+                lines.append(f"  {name}: compliant")
+            else:
+                lines.append(f"  {name}:")
+                for issue in r.issues:
+                    lines.append(
+                        f"    [{issue.severity}] {issue.message} "
+                        f"Fix: {issue.suggestion}"
+                    )
+        return "\n".join(lines)
+
+    def _repr_html_(self) -> str:
+        status = "compliant" if self.is_valid else "non-compliant"
+        parts = [f"<div><b>DatasetValidationReport</b>: {status}"]
+        for name, r in self.reports.items():
+            parts.append(f"<div><b>{html.escape(str(name))}</b></div>")
+            parts.append(r._repr_html_())
+        parts.append("</div>")
+        return "".join(parts)
+
+
+def _is_real_numeric(dtype) -> bool:
+    return np.issubdtype(dtype, np.number) and not np.issubdtype(
+        dtype, np.complexfloating
+    )
+
+
+def _discover_crs(agg):
+    """Return a CRS for *agg*, or None. Mirrors the polygonize convention.
+
+    Resolution order: ``attrs['crs']``, ``attrs['crs_wkt']``, then
+    ``agg.rio.crs`` if rioxarray is installed. The xrspatial.geotiff
+    "no georeference" marker forces None.
+    """
+    if agg.attrs.get("_xrspatial_no_georef"):
+        return None
+    crs = agg.attrs.get("crs")
+    if crs is not None:
+        return crs
+    crs_wkt = agg.attrs.get("crs_wkt")
+    if crs_wkt is not None:
+        return crs_wkt
+    try:
+        rio_crs = agg.rio.crs
+    except Exception:
+        return None
+    return rio_crs
+
+
+def _coord_values(coord):
+    """1-D coordinate values as a host numpy array, or None.
+
+    Coordinates are metadata, not the data buffer, so reading them does
+    not materialize a dask/cupy raster. A cupy-backed coordinate is
+    pulled to host with ``.get()``.
+    """
+    if coord is None or coord.ndim != 1:
+        return None
+    data = coord.data
+    get = getattr(data, "get", None)
+    if get is not None and not isinstance(data, np.ndarray):
+        # cupy ndarray -> host
+        try:
+            data = get()
+        except TypeError:
+            data = np.asarray(coord.values)
+    values = np.asarray(data)
+    if not np.issubdtype(values.dtype, np.number):
+        return None
+    return values
+
+
+def validate(agg) -> ValidationReport:
+    """Check a DataArray against the xarray-spatial input contract.
+
+    Parameters
+    ----------
+    agg : xarray.DataArray
+        Raster to check.
+
+    Returns
+    -------
+    ValidationReport
+        Lists every contract violation as an error (a spatial op will
+        fail) or a warning (behavior degrades), each with a suggested
+        fix. The report is truthy when there are no error-level issues.
+
+    Examples
+    --------
+    .. sourcecode:: python
+
+        >>> import numpy as np, xarray as xr
+        >>> import xrspatial  # registers the .xrs accessor
+        >>> bad = xr.DataArray(
+        ...     np.array([["a", "b"], ["c", "d"]]), dims=["y", "x"]
+        ... )
+        >>> report = bad.xrs.validate()
+        >>> bool(report)
+        False
+    """
+    issues = []
+
+    if not isinstance(agg, xr.DataArray):
+        issues.append(
+            ValidationIssue(
+                "error",
+                "type",
+                f"input is a {type(agg).__name__}, not an xarray.DataArray.",
+                "wrap the array with xr.DataArray(...).",
+            )
+        )
+        return ValidationReport(issues, name=None)
+
+    name = agg.name
+
+    # --- ndim ---
+    if agg.ndim not in (2, 3):
+        issues.append(
+            ValidationIssue(
+                "error",
+                "ndim",
+                f"array is {agg.ndim}D; spatial ops need 2D, or 3D with a "
+                "leading band/time axis.",
+                "reduce to 2D/3D, e.g. select a band/time slice with "
+                "da.isel(...).",
+            )
+        )
+
+    # --- dtype ---
+    if not _is_real_numeric(agg.dtype):
+        issues.append(
+            ValidationIssue(
+                "error",
+                "dtype",
+                f"dtype is {agg.dtype}; spatial ops need a real numeric "
+                "(integer or float) dtype.",
+                "cast with da.astype('float64').",
+            )
+        )
+
+    # Spatial-axis checks only make sense once there are >= 2 dims.
+    if agg.ndim >= 2:
+        ydim, xdim = agg.dims[-2], agg.dims[-1]
+
+        # --- spatial dim names ---
+        if (
+            str(ydim).lower() not in _Y_DIM_NAMES
+            or str(xdim).lower() not in _X_DIM_NAMES
+        ):
+            issues.append(
+                ValidationIssue(
+                    "warning",
+                    "spatial_dims",
+                    f"last two dims are ('{ydim}', '{xdim}'); xrspatial "
+                    "expects spatial axes named like ('y', 'x') or "
+                    "('lat', 'lon').",
+                    "rename with da.rename({'%s': 'y', '%s': 'x'})."
+                    % (ydim, xdim),
+                )
+            )
+
+        _check_spatial_coords(agg, ydim, xdim, issues)
+
+        # --- crs ---
+        if _discover_crs(agg) is None:
+            issues.append(
+                ValidationIssue(
+                    "warning",
+                    "crs",
+                    "no CRS found; geodesic methods and GeoTIFF I/O cannot "
+                    "georeference the raster.",
+                    "set one, e.g. da.attrs['crs'] = 'EPSG:4326'.",
+                )
+            )
+
+    return ValidationReport(issues, name=name)
+
+
+def _check_spatial_coords(agg, ydim, xdim, issues):
+    """Coordinate-dependent checks: presence, monotonic, spacing, cellsize."""
+    have_both = True
+    for dim in (ydim, xdim):
+        # A bare dimension is absent from ``coords`` even though
+        # ``coords.get(dim)`` synthesizes a default range index, so test
+        # membership explicitly.
+        if dim not in agg.coords:
+            have_both = False
+            issues.append(
+                ValidationIssue(
+                    "warning",
+                    "coords_present",
+                    f"dim '{dim}' has no coordinate; resolution-dependent "
+                    "ops (slope, hillshade, proximity) cannot infer cell "
+                    "size.",
+                    f"assign a coordinate, e.g. "
+                    f"da.assign_coords({dim}=np.arange(da.sizes['{dim}'])).",
+                )
+            )
+            continue
+        coord = agg.coords[dim]
+        if not np.issubdtype(coord.dtype, np.number):
+            have_both = False
+            issues.append(
+                ValidationIssue(
+                    "error",
+                    "coords_numeric",
+                    f"'{dim}' coordinate is {coord.dtype}, not numeric.",
+                    f"cast it, e.g. da['{dim}'] = "
+                    f"da['{dim}'].astype('float64').",
+                )
+            )
+            continue
+
+        values = _coord_values(coord)
+        if values is None or values.size < 2:
+            continue
+
+        diffs = np.diff(values)
+        if not np.all(np.isfinite(diffs)):
+            continue
+
+        # --- monotonic ---
+        if not (np.all(diffs > 0) or np.all(diffs < 0)):
+            have_both = False
+            issues.append(
+                ValidationIssue(
+                    "error",
+                    "monotonic",
+                    f"'{dim}' coordinate is not strictly monotonic.",
+                    f"sort with da.sortby('{dim}').",
+                )
+            )
+            continue
+
+        # --- even spacing --- (mirrors utils._warn_if_irregular_spacing)
+        if values.size >= 3:
+            step = abs(diffs[0])
+            if not np.allclose(np.abs(diffs), step, rtol=1e-5, atol=0):
+                issues.append(
+                    ValidationIssue(
+                        "warning",
+                        "even_spacing",
+                        f"'{dim}' coordinate is not evenly spaced; "
+                        "distance-based results may be inaccurate.",
+                        "resample to a regular grid, or set an explicit "
+                        "da.attrs['res'] = (dx, dy).",
+                    )
+                )
+
+    # --- cellsize --- needs both spatial coords usable.
+    if have_both:
+        _check_cellsize(agg, ydim, xdim, issues)
+
+    _check_geographic_range(agg, ydim, xdim, issues)
+
+
+def _check_cellsize(agg, ydim, xdim, issues):
+    """Cell size must be positive and finite on both axes."""
+    yv = _coord_values(agg.coords.get(ydim))
+    xv = _coord_values(agg.coords.get(xdim))
+    if yv is None or xv is None or yv.size < 2 or xv.size < 2:
+        return
+    cellsize_y = (yv[-1] - yv[0]) / (yv.size - 1)
+    cellsize_x = (xv[-1] - xv[0]) / (xv.size - 1)
+    for axis, cellsize in (("x", cellsize_x), ("y", cellsize_y)):
+        if not np.isfinite(cellsize) or cellsize == 0:
+            issues.append(
+                ValidationIssue(
+                    "error",
+                    "cellsize",
+                    f"{axis}-axis cell size is {cellsize}; planar ops need "
+                    "a positive, finite cell size.",
+                    "fix the coordinate spacing, or set "
+                    "da.attrs['res'] = (dx, dy).",
+                )
+            )
+
+
+def _check_geographic_range(agg, ydim, xdim, issues):
+    """When dims look like lat/lon, values must be in geographic range."""
+    is_lat = str(ydim).lower() in _LAT_NAMES
+    is_lon = str(xdim).lower() in _LON_NAMES
+    if not (is_lat or is_lon):
+        return
+    if is_lat:
+        lat = _coord_values(agg.coords.get(ydim))
+        if lat is not None and lat.size:
+            lo, hi = float(np.nanmin(lat)), float(np.nanmax(lat))
+            if lo < -90 or hi > 90 or (hi - lo) > 180:
+                issues.append(
+                    ValidationIssue(
+                        "warning",
+                        "geographic_range",
+                        f"'{ydim}' values [{lo}, {hi}] fall outside "
+                        "latitude range [-90, 90]; coordinates look "
+                        "projected, not geographic.",
+                        "use geographic (lat/lon) coordinates for geodesic "
+                        "methods, or rename the projected axis.",
+                    )
+                )
+    if is_lon:
+        lon = _coord_values(agg.coords.get(xdim))
+        if lon is not None and lon.size:
+            lo, hi = float(np.nanmin(lon)), float(np.nanmax(lon))
+            if lo < -180 or hi > 360 or (hi - lo) > 360:
+                issues.append(
+                    ValidationIssue(
+                        "warning",
+                        "geographic_range",
+                        f"'{xdim}' values [{lo}, {hi}] fall outside "
+                        "longitude range [-180, 360]; coordinates look "
+                        "projected, not geographic.",
+                        "use geographic (lat/lon) coordinates for geodesic "
+                        "methods, or rename the projected axis.",
+                    )
+                )
+
+
+def validate_dataset(ds) -> DatasetValidationReport:
+    """Check each data variable of a Dataset against the contract.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        Dataset whose data variables are rasters.
+
+    Returns
+    -------
+    DatasetValidationReport
+        A per-variable map of :class:`ValidationReport`. Truthy when
+        every variable is compliant.
+    """
+    reports = {name: validate(ds[name]) for name in ds.data_vars}
+    return DatasetValidationReport(reports)

From 9c0f38816f63136307533ba271768cbc9934eb88 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Wed, 24 Jun 2026 14:25:32 -0700
Subject: [PATCH 2/3] Document validate in utilities reference and README
 matrix (#3485)

---
 README.md                           |  1 +
 docs/source/reference/utilities.rst | 10 ++++++++++
 2 files changed, 11 insertions(+)

diff --git a/README.md b/README.md
index 1183870ea..0f9b0cb91 100644
--- a/README.md
+++ b/README.md
@@ -266,6 +266,7 @@ Built-in Numba JIT and CUDA projection kernels bypass pyproj for per-pixel coord
 | [rechunk_no_shuffle](xrspatial/utils.py) | Rechunk dask arrays using whole-chunk multiples (no shuffle) | Custom | 🔼 | 🔼 | 🔼 | 🔼 |
 | [fused_overlap](xrspatial/utils.py) | Fuse sequential map_overlap calls into a single pass | Custom | 🔼 | 🔼 | 🔼 | 🔼 |
 | [multi_overlap](xrspatial/utils.py) | Run multi-output kernel in a single overlap pass | Custom | 🔼 | 🔼 | 🔼 | 🔼 |
+| [validate](xrspatial/validate.py) | Check a raster against the xarray-spatial input contract (`.xrs.validate()`) | Custom | ✅ | 🔼 | 🔼 | 🔼 |
 
 -----------
 
diff --git a/docs/source/reference/utilities.rst b/docs/source/reference/utilities.rst
index 526c7a836..4aaa5331d 100644
--- a/docs/source/reference/utilities.rst
+++ b/docs/source/reference/utilities.rst
@@ -75,3 +75,13 @@ Diagnostics
     :toctree: _autosummary
 
     xrspatial.diagnostics.diagnose
+
+Validation
+==========
+Check a raster against the xarray-spatial input contract. Also available
+on the accessor as ``da.xrs.validate()`` and ``ds.xrs.validate()``.
+
+.. autosummary::
+    :toctree: _autosummary
+
+    xrspatial.validate.validate

From 8413f8bc6af3bc879e609a8a0367ec01a4918435 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Wed, 24 Jun 2026 14:27:57 -0700
Subject: [PATCH 3/3] Address review: guard geographic check on bare dims,
 align even-spacing step (#3485)

---
 xrspatial/tests/test_contract_validate.py | 13 +++++++++++++
 xrspatial/validate.py                     | 20 +++++++++++++-------
 2 files changed, 26 insertions(+), 7 deletions(-)

diff --git a/xrspatial/tests/test_contract_validate.py b/xrspatial/tests/test_contract_validate.py
index 962b94f89..77dd57cc6 100644
--- a/xrspatial/tests/test_contract_validate.py
+++ b/xrspatial/tests/test_contract_validate.py
@@ -150,6 +150,19 @@ def test_projected_latlon_is_warning():
     assert "geographic_range" in _checks(report)
 
 
+def test_bare_latlon_dim_does_not_trigger_geographic_warning():
+    # A dim named 'lat' with no real coordinate must not be checked
+    # against the geographic range using xarray's synthesized indices.
+    da = xr.DataArray(
+        np.zeros((100, 4), dtype="float64"), dims=["lat", "lon"],
+        attrs={"crs": "EPSG:4326"},
+    )
+    report = da.xrs.validate()
+    assert "geographic_range" not in _checks(report)
+    # the real problem (no coords) is still reported
+    assert "coords_present" in _checks(report)
+
+
 def test_missing_crs_is_warning():
     da = _compliant()
     da.attrs = {}
diff --git a/xrspatial/validate.py b/xrspatial/validate.py
index bda1bfcd1..fd40a3ade 100644
--- a/xrspatial/validate.py
+++ b/xrspatial/validate.py
@@ -431,9 +431,10 @@ def _check_spatial_coords(agg, ydim, xdim, issues):
             )
             continue
 
-        # --- even spacing --- (mirrors utils._warn_if_irregular_spacing)
+        # --- even spacing --- (mirrors utils._warn_if_irregular_spacing,
+        # which compares each step to the averaged span/(n-1) resolution)
         if values.size >= 3:
-            step = abs(diffs[0])
+            step = abs(values[-1] - values[0]) / (values.size - 1)
             if not np.allclose(np.abs(diffs), step, rtol=1e-5, atol=0):
                 issues.append(
                     ValidationIssue(
@@ -476,13 +477,18 @@ def _check_cellsize(agg, ydim, xdim, issues):
 
 
 def _check_geographic_range(agg, ydim, xdim, issues):
-    """When dims look like lat/lon, values must be in geographic range."""
-    is_lat = str(ydim).lower() in _LAT_NAMES
-    is_lon = str(xdim).lower() in _LON_NAMES
+    """When dims look like lat/lon, values must be in geographic range.
+
+    Only real coordinates are checked: a bare dimension synthesizes a
+    default integer index whose range would otherwise trip a false
+    "looks projected" warning.
+    """
+    is_lat = str(ydim).lower() in _LAT_NAMES and ydim in agg.coords
+    is_lon = str(xdim).lower() in _LON_NAMES and xdim in agg.coords
     if not (is_lat or is_lon):
         return
     if is_lat:
-        lat = _coord_values(agg.coords.get(ydim))
+        lat = _coord_values(agg.coords[ydim])
         if lat is not None and lat.size:
             lo, hi = float(np.nanmin(lat)), float(np.nanmax(lat))
             if lo < -90 or hi > 90 or (hi - lo) > 180:
@@ -498,7 +504,7 @@ def _check_geographic_range(agg, ydim, xdim, issues):
                     )
                 )
     if is_lon:
-        lon = _coord_values(agg.coords.get(xdim))
+        lon = _coord_values(agg.coords[xdim])
         if lon is not None and lon.size:
             lo, hi = float(np.nanmin(lon)), float(np.nanmax(lon))
             if lo < -180 or hi > 360 or (hi - lo) > 360: