diff --git a/README.md b/README.md
index 5ddb40de..6868cc4c 100644
--- a/README.md
+++ b/README.md
@@ -125,12 +125,18 @@ result, and prediction interfaces.
 
 | Formula surface    | Supported terms                                                                                      |
 | ------------------ | ---------------------------------------------------------------------------------------------------- |
-| Univariate smooths | `s(..., bs='bs')`, `cr`, `cs`, `cc`, `cp`, `ds`, `gp`, `ps`, `tp`, `ts`                              |
+| Metric smooths     | `s(..., bs='bs')`, `cr`, `cs`, `cc`, `cp`, `ds`, `gp`, `ps`, `sos`, `tp`, `ts`                       |
 | Structured smooths | Markov random fields `mrf`, random effects `re`, factor smooths `fs`, sum-to-zero factor smooths `sz` |
-| Tensor products    | `te(...)` and `ti(...)` over supported numeric or MRF marginals                                      |
+| Tensor products    | `te(...)` and `ti(...)` over supported metric or MRF marginals                                       |
 | Parametric terms   | numeric and factor terms, supported interactions, intercept policies, and formula offsets            |
 | Shared smoothing   | supported `id=` groups, fixed/free smoothing parameters, `select=True`, and `pc=` on supported bases |
 
+For spherical splines, write
+`s(latitude, longitude, bs='sos')`: coordinates are degrees in latitude-first,
+longitude-second order. A spherical marginal inside `te()` or `ti()` needs
+`d=2`. Because the array API creates one smooth per input column, joint SOS
+terms use the formula interface.
+
 
 ### Shape-constrained functionality
 
diff --git a/docs/api/gam.rst b/docs/api/gam.rst
index 506ccd91..7f7e524a 100644
--- a/docs/api/gam.rst
+++ b/docs/api/gam.rst
@@ -30,6 +30,21 @@ High-level model
    model.summary()
    model.plot()
 
+Spherical smooths
+-----------------
+
+Use ``s(latitude, longitude, bs="sos")`` for an isotropic smooth on a
+sphere.  Coordinates are supplied in degrees, latitude first and longitude
+second.  The default ``m=0`` is the second-order Wendelberger spline; integer
+orders from ``-2`` through ``4`` select the upstream Duchon or Wahba kernel
+branches.  SOS margins in ``te`` and ``ti`` must be grouped with ``d=2``.
+The array API cannot express this joint two-coordinate term, so SOS models use
+the formula interface.  The ``m=-1`` null space is four-dimensional; combining
+that order with an ``fs`` factor smooth is rejected because the corresponding
+upstream penalty split is LAPACK-orientation dependent.  Upstream SOS smooths
+also do not define derivative matrices, and NAMpy does not expose the
+hemisphere-specific ``plot.gam`` schemes 0 and 1 through its generic plotter.
+
 Shape-constrained smooths
 -------------------------
 
diff --git a/docs/generate_notebooks.py b/docs/generate_notebooks.py
index 332025d9..e5e1700c 100644
--- a/docs/generate_notebooks.py
+++ b/docs/generate_notebooks.py
@@ -1515,6 +1515,7 @@ def gam_notebook() -> dict:
 | `cc` | cyclic cubic spline for periodic covariates |
 | `ps` | P-spline with difference penalties |
 | `gp` | low-rank Gaussian-process smooth with spherical, power-exponential, or Matérn covariance |
+| `sos` | isotropic spherical spline for latitude then longitude in degrees |
 | `mrf` | region effect coupled by a neighbor graph, polygon boundary, or supplied penalty |
 | `tp`, `ts` | thin-plate regression spline; `ts` adds shrinkage |
 | `te(...)` | scale-invariant tensor product including main-effect directions |
@@ -1540,6 +1541,9 @@ def gam_notebook() -> dict:
     "gaussian_process": GAM(
         formula="demand ~ s(temperature, humidity, bs='gp', k=20, m=[3, 1.0])"
     ),
+    "spherical_spline": GAM(
+        formula="demand ~ s(latitude, longitude, bs='sos', k=20, m=0)"
+    ),
     "markov_random_field": GAM(
         formula=(
             "demand ~ s(region, bs='mrf', "
diff --git a/docs/notebooks/01_gam.ipynb b/docs/notebooks/01_gam.ipynb
index 6280af49..71f9c04e 100644
--- a/docs/notebooks/01_gam.ipynb
+++ b/docs/notebooks/01_gam.ipynb
@@ -630,6 +630,7 @@
     "| `cc` | cyclic cubic spline for periodic covariates |\n",
     "| `ps` | P-spline with difference penalties |\n",
     "| `gp` | low-rank Gaussian-process smooth with spherical, power-exponential, or Mat\u00e9rn covariance |\n",
+    "| `sos` | isotropic spherical spline for latitude then longitude in degrees |\n",
     "| `mrf` | region effect coupled by a neighbor graph, polygon boundary, or supplied penalty |\n",
     "| `tp`, `ts` | thin-plate regression spline; `ts` adds shrinkage |\n",
     "| `te(...)` | scale-invariant tensor product including main-effect directions |\n",
@@ -682,6 +683,9 @@
     "    \"gaussian_process\": GAM(\n",
     "        formula=\"demand ~ s(temperature, humidity, bs='gp', k=20, m=[3, 1.0])\"\n",
     "    ),\n",
+    "    \"spherical_spline\": GAM(\n",
+    "        formula=\"demand ~ s(latitude, longitude, bs='sos', k=20, m=0)\"\n",
+    "    ),\n",
     "    \"markov_random_field\": GAM(\n",
     "        formula=(\n",
     "            \"demand ~ s(region, bs='mrf', \"\n",
diff --git a/nampy/gam/compiler/factory.py b/nampy/gam/compiler/factory.py
index 122f2479..e04cfaf2 100644
--- a/nampy/gam/compiler/factory.py
+++ b/nampy/gam/compiler/factory.py
@@ -21,6 +21,7 @@
 from ..smooths.univariate.ds import DuchonSplineTerm
 from ..smooths.univariate.gp import GaussianProcessTerm
 from ..smooths.univariate.ps import PSplineTerm1D
+from ..smooths.univariate.sos import SphericalSplineTerm
 from ..specs import LinearPredictorSpec, PenaltyGroupSpec, TermSpec
 from ..specs.smooth import (
     CubicRegressionSmoothSpec,
@@ -34,6 +35,7 @@
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
     ShapeConstrainedSmoothSpec,
+    SphericalSplineSmoothSpec,
     SumToZeroFactorSmoothSpec,
     TensorInteractionSmoothSpec,
     TensorProductSmoothSpec,
@@ -230,6 +232,25 @@ def instantiate_term(term_like: TermSpec | Any):
             metadata=metadata,
         )
 
+    if isinstance(smooth_spec, SphericalSplineSmoothSpec):
+        return SphericalSplineTerm(
+            feature=features,
+            k=smooth_spec.k,
+            m=smooth_spec.m,
+            label=label,
+            term_id=term_like.term_id,
+            smoothing_id=smoothing_id,
+            by=by,
+            sp=smooth_spec.sp,
+            select=smooth_spec.select,
+            fixed=smooth_spec.fx,
+            constraint_mode=smooth_spec.constraint_mode,
+            pc=smooth_spec.pc,
+            knots=smooth_spec.knots,
+            xt=smooth_spec.xt,
+            metadata=metadata,
+        )
+
     if isinstance(smooth_spec, MarkovRandomFieldSmoothSpec):
         return MarkovRandomFieldTerm(
             feature=features,
diff --git a/nampy/gam/diagnostics/plots.py b/nampy/gam/diagnostics/plots.py
index d858ed27..800fe625 100644
--- a/nampy/gam/diagnostics/plots.py
+++ b/nampy/gam/diagnostics/plots.py
@@ -190,6 +190,12 @@ def _prepare_smooth(model, tb, *, se, n, n2, n3, xlab, ylab, main, label,
     basis_name = str(getattr(tb, "basis_name", "")).lower()
     X_train = np.asarray(model.X_)
 
+    if basis_name == "sos":
+        raise NotImplementedError(
+            "plot() for bs='sos' requires mgcv's rotated hemisphere projection; "
+            "the generic rectangular 2D plot is not equivalent."
+        )
+
     if basis_name == "re":
         # plot.random.effect (plots.r:357-367): X is the identity; the plot
         # is a normal QQ plot of the estimated effects.
diff --git a/nampy/gam/smooths/__init__.py b/nampy/gam/smooths/__init__.py
index 65c16552..b77e568c 100644
--- a/nampy/gam/smooths/__init__.py
+++ b/nampy/gam/smooths/__init__.py
@@ -25,6 +25,7 @@
 from .univariate.ds import DuchonSplineTerm
 from .univariate.gp import GaussianProcessTerm
 from .univariate.ps import PSplineTerm1D
+from .univariate.sos import SphericalSplineTerm
 from .univariate.tp import ThinPlateSplineTerm
 
 # mgcv-facing smooth aliases keep formulas/tests readable without reintroducing
@@ -36,6 +37,7 @@
 cr = cs = cc = CubicSplineTerm
 ds = DuchonSplineTerm
 gp = GaussianProcessTerm
+sos = SphericalSplineTerm
 mrf = MarkovRandomFieldTerm
 cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
@@ -65,6 +67,7 @@
     "CubicSplineTerm",
     "DuchonSplineTerm",
     "GaussianProcessTerm",
+    "SphericalSplineTerm",
     "MarkovRandomFieldTerm",
     "DerivativeBSplineTerm1D",
     "PSplineTerm1D",
@@ -83,6 +86,7 @@
     "cc",
     "ds",
     "gp",
+    "sos",
     "mrf",
     "cp",
     "ps",
diff --git a/nampy/gam/smooths/categorical/fs.py b/nampy/gam/smooths/categorical/fs.py
index 6dd0a113..87036fce 100644
--- a/nampy/gam/smooths/categorical/fs.py
+++ b/nampy/gam/smooths/categorical/fs.py
@@ -20,6 +20,7 @@
 from ..univariate.ds import DuchonSplineTerm
 from ..univariate.gp import GaussianProcessTerm
 from ..univariate.ps import PSplineTerm1D
+from ..univariate.sos import SphericalSplineTerm
 from .categorical_utils import (
     as_object_1d,
     factor_indicator_matrix,
@@ -111,7 +112,7 @@ def _build_base_smooth_term(
     Build the per-level base smooth used inside fs/sz.
 
     Supported base smooth classes in the current codebase:
-    bs, cr, cs, cc, cp, ds, gp, mrf, ps, tp, ts
+    bs, cr, cs, cc, cp, ds, gp, mrf, ps, sos, tp, ts
     """
     base_bs = str(base_bs).lower()
     metric_features = list(metric_features)
@@ -122,10 +123,10 @@ def _build_base_smooth_term(
     if mode == "fs" and base_bs in {"cs", "ts"}:
         raise NotImplementedError(_fs_full_rank_base_error(base_bs))
 
-    if len(metric_features) > 1 and base_bs not in {"ds", "gp", "tp", "ts"}:
+    if len(metric_features) > 1 and base_bs not in {"ds", "gp", "sos", "tp", "ts"}:
         raise NotImplementedError(
             f"Current {mode} implementation supports multivariate base smooths only "
-            f"for bs in {{'ds','gp','tp','ts'}}, got base bs={base_bs!r}."
+            f"for bs in {{'ds','gp','sos','tp','ts'}}, got base bs={base_bs!r}."
         )
 
     if xt_rest is not None and base_bs not in {
@@ -135,11 +136,12 @@ def _build_base_smooth_term(
         "gp",
         "mrf",
         "ps",
+        "sos",
         "tp",
         "ts",
     }:
         raise NotImplementedError(
-            "Extra xt options are currently only supported for bs/cp/ds/gp/mrf/ps/tp/ts "
+            "Extra xt options are currently only supported for bs/cp/ds/gp/mrf/ps/sos/tp/ts "
             "base smooths, "
             f"got xt={xt_rest!r} with base bs={base_bs!r}."
         )
@@ -242,6 +244,24 @@ def _build_base_smooth_term(
             metadata=metadata,
         )
 
+    if base_bs == "sos":
+        return SphericalSplineTerm(
+            feature=metric_features,
+            k=k,
+            m=outer_m,
+            label=label,
+            smoothing_id=None,
+            by=by,
+            sp=None,
+            select=bool(select),
+            fixed=bool(fixed),
+            constraint_mode=str(constraint_mode),
+            pc=None,
+            knots=knots,
+            xt=xt_rest,
+            metadata=metadata,
+        )
+
     if base_bs == "mrf":
         return MarkovRandomFieldTerm(
             feature=metric_features[0],
@@ -280,7 +300,7 @@ def _build_base_smooth_term(
 
     raise NotImplementedError(
         f"Current {mode} implementation supports base bs in "
-        f"{{'bs','cr','cs','cc','cp','ds','gp','mrf','ps','tp','ts'}}, got {base_bs!r}."
+        f"{{'bs','cr','cs','cc','cp','ds','gp','mrf','ps','sos','tp','ts'}}, got {base_bs!r}."
     )
 
 
@@ -291,6 +311,8 @@ def _penalty_rank_from_base_term(base_term, basis_matrix, penalty_matrix) -> int
         return int(base_term._setup.rank)
     if isinstance(base_term, GaussianProcessTerm):
         return int(base_term._setup.rank)
+    if isinstance(base_term, SphericalSplineTerm):
+        return int(base_term._setup.rank)
     if isinstance(base_term, MarkovRandomFieldTerm):
         return int(base_term._setup.rank)
     if isinstance(base_term, PSplineTerm1D) and len(base_term.penalties) > 0:
@@ -697,6 +719,15 @@ def fit(self, X, feature_names):
             )
 
         self._base_term = base_term
+        if (
+            isinstance(base_term, SphericalSplineTerm)
+            and int(base_term._setup.null_space_dim) > 1
+        ):
+            raise NotImplementedError(
+                "bs='fs' with an SOS m=-1 base is not enabled: mgcv's four-way "
+                "repeated null eigenspace receives separate penalties whose "
+                "orientation is LAPACK-dependent. Use another SOS order."
+            )
         if (
             isinstance(base_term, MarkovRandomFieldTerm)
             and base_term._setup.used_low_rank
diff --git a/nampy/gam/smooths/tensor/marginals.py b/nampy/gam/smooths/tensor/marginals.py
index cec7fb08..8c523428 100644
--- a/nampy/gam/smooths/tensor/marginals.py
+++ b/nampy/gam/smooths/tensor/marginals.py
@@ -13,10 +13,11 @@
 from ..univariate.ds import DuchonSplineTerm
 from ..univariate.gp import GaussianProcessTerm
 from ..univariate.ps import PSplineTerm1D
+from ..univariate.sos import SphericalSplineTerm
 from ..univariate.tp import ThinPlateSplineTerm
 
 TENSOR_MARGINAL_BASES = frozenset(
-    {"bs", "cr", "cs", "cc", "cp", "ds", "gp", "mrf", "ps", "tp", "ts"}
+    {"bs", "cr", "cs", "cc", "cp", "ds", "gp", "mrf", "ps", "sos", "tp", "ts"}
 )
 
 
@@ -147,6 +148,27 @@ def make_tensor_marginal_term(
             metadata=metadata,
         )
 
+    if basis == "sos":
+        if len(marginal_features) != 2:
+            raise ValueError(
+                "Tensor marginal basis 'sos' requires a two-feature group "
+                "(latitude, longitude); supply d=2 for that marginal."
+            )
+        return SphericalSplineTerm(
+            feature=marginal_features,
+            k=k,
+            m=m,
+            xt=xt,
+            label=str(feature),
+            smoothing_id=None,
+            by=None,
+            select=False,
+            fixed=False,
+            constraint_mode=constraint_mode,
+            knots=knots,
+            metadata=metadata,
+        )
+
     if basis == "mrf":
         if len(marginal_features) != 1:
             raise ValueError("Tensor marginal basis 'mrf' only handles one feature.")
diff --git a/nampy/gam/smooths/univariate/__init__.py b/nampy/gam/smooths/univariate/__init__.py
index 290e0a69..636418bd 100644
--- a/nampy/gam/smooths/univariate/__init__.py
+++ b/nampy/gam/smooths/univariate/__init__.py
@@ -3,12 +3,14 @@
 from .ds import DuchonSplineTerm
 from .gp import GaussianProcessTerm
 from .ps import PSplineTerm1D
+from .sos import SphericalSplineTerm
 from .tp import ThinPlateSplineTerm
 
 bs = DerivativeBSplineTerm1D
 cr = cs = cc = CubicSplineTerm
 ds = DuchonSplineTerm
 gp = GaussianProcessTerm
+sos = SphericalSplineTerm
 cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
 
@@ -17,7 +19,8 @@
     "CubicSplineTerm",
     "DuchonSplineTerm",
     "GaussianProcessTerm",
+    "SphericalSplineTerm",
     "PSplineTerm1D",
     "ThinPlateSplineTerm",
 ]
-__all__ += ["bs", "cr", "cs", "cc", "cp", "ds", "gp", "ps", "tp", "ts"]
+__all__ += ["bs", "cr", "cs", "cc", "cp", "ds", "gp", "sos", "ps", "tp", "ts"]
diff --git a/nampy/gam/smooths/univariate/sos.py b/nampy/gam/smooths/univariate/sos.py
new file mode 100644
index 00000000..8c5afec1
--- /dev/null
+++ b/nampy/gam/smooths/univariate/sos.py
@@ -0,0 +1,247 @@
+"""Spherical spline smooth term (``bs='sos'``)."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from ...constraints.absorption import (
+    fit_single_penalty_with_constraint_policy,
+    fit_single_penalty_with_setup_basis,
+)
+from ...penalties.algebra import penalty_rescale_factor, scale_penalty
+from ...splines.univariate.sos import (
+    build_spherical_spline_setup,
+    predict_spherical_spline,
+)
+from ..registry import register_smooth
+from ..smooth_base import BaseSmoothTerm, _resolve_feature, columns_as_float_matrix
+
+
+@register_smooth("sos")
+class SphericalSplineTerm(BaseSmoothTerm):
+    term_type = "smooth"
+    basis_name = "sos"
+    supports_tensor_marginal = True
+
+    def __init__(
+        self,
+        feature,
+        k=-1,
+        m=None,
+        label=None,
+        term_id=None,
+        smoothing_id=None,
+        by=None,
+        sp=None,
+        select=False,
+        fixed=False,
+        constraint_mode="auto",
+        pc=None,
+        knots=None,
+        xt=None,
+        null_penalty_tol=1e-10,
+        metadata=None,
+    ):
+        features = list(feature) if not isinstance(feature, (str, int)) else [feature]
+        super().__init__(
+            feature=features,
+            label=label or f"s({', '.join(map(str, features))})",
+            term_id=term_id,
+            smoothing_id=smoothing_id,
+            by=by,
+            sp=sp,
+            metadata=metadata,
+        )
+        self.k = int(k)
+        self.m = m
+        self.select = bool(select)
+        self.fixed = bool(fixed)
+        self.constraint_mode = str(constraint_mode).lower()
+        self.pc = pc
+        self.knots = knots
+        self.xt = xt
+        self.null_penalty_tol = float(null_penalty_tol)
+        if len(features) != 2:
+            raise ValueError(
+                "Can only deal with a sphere: bs='sos' requires exactly two "
+                "features, latitude first and longitude second."
+            )
+        if self.constraint_mode not in {"auto", "factor_by", "always", "never"}:
+            raise ValueError(
+                "constraint_mode must be one of "
+                "{'auto', 'factor_by', 'always', 'never'}."
+            )
+
+        self._feature_indices = None
+        self._feature_names = None
+        self._by_state = None
+        self._basis_train = None
+        self._penalties = None
+        self._setup = None
+
+    @property
+    def expected_linked_penalty_count(self):
+        return None if self.select else 1
+
+    def fit(self, X, feature_names):
+        feature_indices = []
+        resolved_names = []
+        for feature in self.feature:
+            index, name = _resolve_feature(feature, feature_names)
+            feature_indices.append(index)
+            resolved_names.append(name)
+        values = columns_as_float_matrix(X, feature_indices)
+        self._set_by_state(X, feature_names)
+        self._feature_indices = feature_indices
+        self._feature_names = resolved_names
+        self._set_resolved_features(resolved_names)
+
+        shared_X = self._linked_id_setup_matrix(feature_names)
+        setup_values = (
+            columns_as_float_matrix(shared_X, feature_indices)
+            if shared_X is not None
+            else values
+        )
+        self._setup = build_spherical_spline_setup(
+            setup_values,
+            k=self.k,
+            m=self.m,
+            knots=self.knots,
+            xt=self.xt,
+        )
+        setup_base = np.asarray(self._setup.basis_train, dtype=np.float64)
+        base = (
+            setup_base
+            if shared_X is None
+            else np.asarray(
+                predict_spherical_spline(values, self._setup), dtype=np.float64
+            )
+        )
+        raw_penalty = np.asarray(self._setup.penalty, dtype=np.float64)
+        penalty = scale_penalty(setup_base, raw_penalty)
+        self._set_penalty_rescale_factors(
+            [penalty_rescale_factor(setup_base, raw_penalty)]
+        )
+
+        if self.pc is not None:
+            constrained_basis, constrained_penalties, transform, _ = (
+                self._apply_point_constraint(
+                    base,
+                    [penalty],
+                    self.pc,
+                    feature_names=self._feature_names,
+                    point_basis_fn=lambda points: predict_spherical_spline(
+                        points, self._setup
+                    )[0],
+                    fixed=self.fixed,
+                )
+            )
+            self._basis_train = np.asarray(constrained_basis, dtype=np.float64)
+            self._penalties = constrained_penalties
+            self._record_constraint_result("pc", transform, absorbed_by="runtime")
+            return self
+
+        auto_constrain = bool(self._by_state.is_constant)
+        if shared_X is None:
+            result = fit_single_penalty_with_constraint_policy(
+                base,
+                penalty,
+                self._by_state,
+                constraint_mode=self.constraint_mode,
+                fixed=self.fixed,
+                auto_constrain_when=auto_constrain,
+            )
+        else:
+            result = fit_single_penalty_with_setup_basis(
+                base,
+                setup_base,
+                penalty,
+                self._by_state,
+                constraint_mode=self.constraint_mode,
+                fixed=self.fixed,
+                auto_constrain_when=auto_constrain,
+            )
+        self._basis_train = result.basis_train
+        self._penalties = result.penalties
+        self._record_constraint_result(
+            result.constraint_kind,
+            result.constraint_transform,
+            absorbed_by=(
+                "runtime" if result.constraint_transform is not None else None
+            ),
+        )
+        return self
+
+    def get_penalty_definitions(self):
+        self._require_fitted()
+        if not self.penalties:
+            return []
+        metadata = {
+            "term_type": self.term_type,
+            "basis_name": self.basis_name,
+            "feature": list(self.feature),
+            "label": self.label,
+            "by": self.by,
+            "by_name": self._by_state.feature_name,
+            "by_is_constant": bool(self._by_state.is_constant),
+            "constraint_mode": self.constraint_mode,
+            "constraint_kind": self.constraint_kind,
+            "pc": self.pc,
+            "knots": self.knots,
+            "xt": self.xt,
+            "m": self.m,
+            "spherical_order": int(self._setup.order),
+            "original_null_space_dim": int(self._setup.null_space_dim),
+            "fixed": bool(self.fixed),
+        }
+        selection_metadata = {**metadata, "is_selection_penalty": True}
+        metadata = self._penalty_metadata_with_scale(metadata, penalty_index=0)
+        return self._build_penalty_block(
+            self.penalties[0],
+            rank=int(self._setup.rank),
+            smooth_metadata=metadata,
+            selection_metadata=selection_metadata,
+        )
+
+    def transform_new(self, X_new):
+        self._require_fitted()
+        values = columns_as_float_matrix(X_new, self._feature_indices)
+        basis = predict_spherical_spline(values, self._setup)
+        return self._apply_constraint_transform_and_by(basis, X_new)
+
+    def tensor_marginal_fit_matrices(
+        self, *, centered=False, apply_np=False, x_train=None
+    ):
+        del apply_np, x_train
+        self._require_fitted()
+        setup_base = np.asarray(self._setup.basis_train, dtype=np.float64)
+        setup_penalty = np.asarray(self._setup.penalty, dtype=np.float64)
+        if centered:
+            if (
+                self._linked_id_setup() is not None
+                and self.constraint_transform is not None
+            ):
+                transform = np.asarray(self.constraint_transform, dtype=np.float64)
+                scaled = scale_penalty(setup_base, setup_penalty)
+                return (
+                    np.asarray(setup_base @ transform, dtype=np.float64),
+                    np.asarray(transform.T @ scaled @ transform, dtype=np.float64),
+                    None,
+                )
+            return super().tensor_marginal_fit_matrices(centered=True)
+        return setup_base, setup_penalty, None
+
+    def tensor_marginal_predict_matrix(
+        self, X_new, *, centered=False, np_transform=None
+    ):
+        if centered:
+            basis = np.asarray(self.transform_new(X_new), dtype=np.float64)
+        else:
+            values = columns_as_float_matrix(X_new, self._feature_indices)
+            basis = predict_spherical_spline(values, self._setup)
+        if np_transform is not None:
+            basis = basis @ np.asarray(np_transform, dtype=np.float64)
+        return np.asarray(basis, dtype=np.float64)
+
+
+__all__ = ["SphericalSplineTerm"]
diff --git a/nampy/gam/specs/__init__.py b/nampy/gam/specs/__init__.py
index 0b849c2e..ae5855a0 100644
--- a/nampy/gam/specs/__init__.py
+++ b/nampy/gam/specs/__init__.py
@@ -17,6 +17,7 @@
     RandomEffectSmoothSpec,
     ShapeConstrainedSmoothSpec,
     SmoothSpec,
+    SphericalSplineSmoothSpec,
     SumToZeroFactorSmoothSpec,
     TensorInteractionSmoothSpec,
     TensorProductSmoothSpec,
@@ -37,6 +38,7 @@
     "DuchonSplineSmoothSpec",
     "FactorSmoothInteractionSpec",
     "GaussianProcessSmoothSpec",
+    "SphericalSplineSmoothSpec",
     "MarkovRandomFieldSmoothSpec",
     "PSplineSmoothSpec",
     "RandomEffectSmoothSpec",
diff --git a/nampy/gam/specs/build.py b/nampy/gam/specs/build.py
index 66d980f1..c6ff1a31 100644
--- a/nampy/gam/specs/build.py
+++ b/nampy/gam/specs/build.py
@@ -477,8 +477,8 @@ def _fs_by_without_factor_feature_base_spec(smooth_spec: FactorSmoothInteraction
 
     if base_bs == "ps":
         kwargs["m"] = None if xt_rest is None else xt_rest.get("m", None)
-    elif base_bs in {"ds", "gp", "mrf", "tp", "ts"}:
-        if base_bs in {"ds", "gp"}:
+    elif base_bs in {"ds", "gp", "mrf", "sos", "tp", "ts"}:
+        if base_bs in {"ds", "gp", "sos"}:
             kwargs["m"] = smooth_spec.m
         kwargs["xt"] = xt_rest
     elif xt_rest:
diff --git a/nampy/gam/specs/modeling.py b/nampy/gam/specs/modeling.py
index 60251958..c720aac3 100644
--- a/nampy/gam/specs/modeling.py
+++ b/nampy/gam/specs/modeling.py
@@ -77,6 +77,12 @@ def make_predictor_specs(model, feature_names, *, knots=None):
                     metadata={},
                 )
             )
+        elif basis == "sos":
+            raise NotImplementedError(
+                "The array API builds one smooth per feature and cannot express "
+                "the two-coordinate bs='sos' term; use a formula with both latitude "
+                "and longitude."
+            )
         elif basis in {"tp", "ts"}:
             main_terms.append(
                 TermSpec(
diff --git a/nampy/gam/specs/smooth.py b/nampy/gam/specs/smooth.py
index d06281bf..7b7aae87 100644
--- a/nampy/gam/specs/smooth.py
+++ b/nampy/gam/specs/smooth.py
@@ -76,6 +76,15 @@ class GaussianProcessSmoothSpec(BaseSmoothSpec):
     pc: Any = None
 
 
+@dataclass(frozen=True)
+class SphericalSplineSmoothSpec(BaseSmoothSpec):
+    bs: str = "sos"
+    m: Any = None
+    xt: Any = None
+    constraint_mode: str = "auto"
+    pc: Any = None
+
+
 @dataclass(frozen=True)
 class MarkovRandomFieldSmoothSpec(BaseSmoothSpec):
     bs: str = "mrf"
@@ -159,6 +168,7 @@ class TensorInteractionSmoothSpec(BaseSmoothSpec):
     DerivativeBSplineSmoothSpec,
     DuchonSplineSmoothSpec,
     GaussianProcessSmoothSpec,
+    SphericalSplineSmoothSpec,
     MarkovRandomFieldSmoothSpec,
     CubicShrinkageSmoothSpec,
     PSplineSmoothSpec,
diff --git a/nampy/gam/specs/smooth_build.py b/nampy/gam/specs/smooth_build.py
index 06685d58..b2e2dae7 100644
--- a/nampy/gam/specs/smooth_build.py
+++ b/nampy/gam/specs/smooth_build.py
@@ -21,6 +21,7 @@
     RandomEffectSmoothSpec,
     ShapeConstrainedSmoothSpec,
     SmoothSpec,
+    SphericalSplineSmoothSpec,
     SumToZeroFactorSmoothSpec,
     TensorInteractionSmoothSpec,
     TensorProductSmoothSpec,
@@ -146,6 +147,21 @@ def _build_s_gp(opts) -> GaussianProcessSmoothSpec:
     )
 
 
+def _build_s_sos(opts) -> SphericalSplineSmoothSpec:
+    return SphericalSplineSmoothSpec(
+        special="s",
+        k=opts["k"],
+        fx=opts["fx"],
+        select=opts["select"],
+        sp=opts["sp"],
+        knots=opts["knots"],
+        m=opts["m"],
+        xt=opts["xt"],
+        constraint_mode=opts["constraint_mode"],
+        pc=opts["pc"],
+    )
+
+
 def _build_s_mrf(opts) -> MarkovRandomFieldSmoothSpec:
     return MarkovRandomFieldSmoothSpec(
         special="s",
@@ -250,6 +266,7 @@ def _build_s_sz(opts) -> SumToZeroFactorSmoothSpec:
     "bs": _build_s_bs,
     "ds": _build_s_ds,
     "gp": _build_s_gp,
+    "sos": _build_s_sos,
     "mrf": _build_s_mrf,
     "tp": _build_s_tp,
     "ts": _build_s_ts,
@@ -351,6 +368,7 @@ def _is_vector_fx(fx) -> bool:
     "cs",
     "ds",
     "gp",
+    "sos",
     "ps",
     "tp",
     "ts",
@@ -568,7 +586,7 @@ def _default_k_for_smooth(kind, basis, features, default_k):
         # mgcv::te()/ti() default k to 5^d per marginal. The current
         # Python tensor surface supports one feature per marginal, so d = 1.
         return [5] * len(features)
-    if str(basis).lower() in {"bs", "ds", "gp", "mrf", "tp", "ts"}:
+    if str(basis).lower() in {"bs", "ds", "gp", "mrf", "sos", "tp", "ts"}:
         # mgcv/R/smooth.r::s() leaves k = -1. The basis constructor then
         # resolves its dimension-dependent default (TP/TS use M + 8/27/100;
         # DS uses M + 10/30/100; GP uses d + 1 + 10/30/100). A flat default here would be wrong in more
diff --git a/nampy/gam/splines/univariate/__init__.py b/nampy/gam/splines/univariate/__init__.py
index 650f5474..0ed01236 100644
--- a/nampy/gam/splines/univariate/__init__.py
+++ b/nampy/gam/splines/univariate/__init__.py
@@ -49,6 +49,14 @@
     pspline_knots,
     pspline_predict_matrix,
 )
+from .sos import (
+    SphericalSplineSetup,
+    build_spherical_spline_setup,
+    normalize_spherical_order,
+    predict_spherical_spline,
+    spherical_null_space_dimension,
+    spherical_spline_kernel,
+)
 from .tp import build_tprs_term_setup, predict_tprs_term
 
 __all__ = [
@@ -74,6 +82,12 @@
     "gp_polynomial_basis",
     "normalize_gp_definition",
     "predict_gaussian_process",
+    "SphericalSplineSetup",
+    "build_spherical_spline_setup",
+    "normalize_spherical_order",
+    "predict_spherical_spline",
+    "spherical_null_space_dimension",
+    "spherical_spline_kernel",
     "add_full_rank_shrinkage",
     "bspline_design_matrix",
     "cyclic_cubic_bd",
diff --git a/nampy/gam/splines/univariate/sos.py b/nampy/gam/splines/univariate/sos.py
new file mode 100644
index 00000000..f34d60c2
--- /dev/null
+++ b/nampy/gam/splines/univariate/sos.py
@@ -0,0 +1,327 @@
+"""Low-rank spherical splines matching mgcv ``bs='sos'``."""
+
+from __future__ import annotations
+
+import warnings
+from dataclasses import dataclass
+
+import numpy as np
+from scipy.special import spence
+
+from ...linalg.qr import r_linpack_qr_no_pivot
+from .ds import (
+    _duchon_setup_locations,
+    _normalize_duchon_knots,
+    _r_linpack_qty,
+)
+from .tp import _top_eigensystem
+
+
+def normalize_spherical_order(m) -> int:
+    """Mirror the order normalization in mgcv's SOS constructor."""
+    if m is None:
+        return 0
+    values = np.asarray(m, dtype=object).ravel()
+    if values.size != 1:
+        raise ValueError("For bs='sos', m must be a single numeric value.")
+    try:
+        value = float(values[0])
+    except (TypeError, ValueError) as exc:
+        raise ValueError("For bs='sos', m must be a single numeric value.") from exc
+    if np.isnan(value):
+        return 0
+    if not np.isfinite(value):
+        raise ValueError("For bs='sos', m must be finite.")
+    order = int(np.rint(value))
+    if order < -2:
+        order = -1
+    if order > 4:
+        order = 4
+    return order
+
+
+def spherical_null_space_dimension(order: int) -> int:
+    """Return the literal null-tail size used by mgcv 1.9-4."""
+    return 4 if int(order) == -1 else 1
+
+
+def _spherical_geometry(latitude, longitude, knot_latitude, knot_longitude):
+    latitude = np.deg2rad(np.asarray(latitude, dtype=np.float64).ravel())
+    longitude = np.deg2rad(np.asarray(longitude, dtype=np.float64).ravel())
+    knot_latitude = np.deg2rad(np.asarray(knot_latitude, dtype=np.float64).ravel())
+    knot_longitude = np.deg2rad(np.asarray(knot_longitude, dtype=np.float64).ravel())
+    cosine = np.sin(latitude)[:, None] * np.sin(knot_latitude)[None, :] + np.cos(
+        latitude
+    )[:, None] * np.cos(knot_latitude)[None, :] * np.cos(
+        longitude[:, None] - knot_longitude[None, :]
+    )
+    gamma = np.arccos(np.clip(cosine, -1.0, 1.0))
+    return latitude, longitude, knot_latitude, knot_longitude, gamma
+
+
+def spherical_spline_kernel(X, knots, order=0):
+    """Port mgcv's ``makeR`` and return its kernel and null-space tail."""
+    values = np.asarray(X, dtype=np.float64)
+    setup_knots = np.asarray(knots, dtype=np.float64)
+    if values.ndim != 2 or values.shape[1] != 2:
+        raise ValueError(
+            "Spherical spline data must have exactly two columns: latitude, longitude."
+        )
+    if setup_knots.ndim != 2 or setup_knots.shape[1] != 2:
+        raise ValueError(
+            "Spherical spline knots must have exactly two columns: latitude, longitude."
+        )
+    if not np.isfinite(values).all() or not np.isfinite(setup_knots).all():
+        raise ValueError("Spherical spline data and knots must be finite.")
+
+    order = int(order)
+    latitude, longitude, knot_latitude, knot_longitude, gamma = _spherical_geometry(
+        values[:, 0],
+        values[:, 1],
+        setup_knots[:, 0],
+        setup_knots[:, 1],
+    )
+
+    if order == -2:
+        distance = 2.0 * np.sin(gamma / 2.0)
+        distance = np.maximum(distance, np.finfo(np.float64).tiny * 10.0)
+        kernel = -distance
+        tail = np.ones((values.shape[0], 1), dtype=np.float64)
+        constraint_tail = np.ones((setup_knots.shape[0], 1), dtype=np.float64)
+    elif order == -1:
+        distance = 2.0 * np.sin(gamma / 2.0)
+        distance = np.maximum(distance, np.finfo(np.float64).tiny * 10.0)
+        with np.errstate(divide="ignore", invalid="ignore", under="ignore"):
+            kernel = distance * distance * np.log(distance) / (8.0 * np.pi)
+        kernel = np.nan_to_num(kernel, nan=0.0)
+
+        def _tail(lat, lon):
+            z = np.sin(lat)
+            x = np.cos(lat) * np.sin(lon)
+            y = np.cos(lat) * np.cos(lon)
+            return np.column_stack([np.ones(lat.size), x, y, z])
+
+        tail = _tail(latitude, longitude)
+        constraint_tail = _tail(knot_latitude, knot_longitude)
+    elif order == 0:
+        cosine = np.cos(gamma)
+        argument = np.clip((1.0 + cosine) / 2.0, 0.0, 1.0)
+        kernel = (1.0 - np.pi**2 / 6.0 + spence(1.0 - argument)) / (4.0 * np.pi)
+        tail = np.ones((values.shape[0], 1), dtype=np.float64)
+        constraint_tail = np.ones((setup_knots.shape[0], 1), dtype=np.float64)
+    elif order in {1, 2, 3, 4}:
+        z = 1.0 - np.cos(gamma)
+        z = np.maximum(z, np.finfo(np.float64).eps * 0.0001)
+        W = z / 2.0
+        C = np.sqrt(W)
+        A = np.log(1.0 + 1.0 / C)
+        C = 2.0 * C
+        if order == 1:
+            q = 2.0 * A * W - C + 1.0
+            kernel = (q - 0.5) / (2.0 * np.pi)
+        elif order == 2:
+            W2 = W * W
+            q = A * (6.0 * W2 - 2.0 * W) - 3.0 * C * W + 3.0 * W + 0.5
+            kernel = (q / 2.0 - 1.0 / 6.0) / (2.0 * np.pi)
+        elif order == 3:
+            W2 = W * W
+            W3 = W2 * W
+            q = (
+                A * (60.0 * W3 - 36.0 * W2)
+                + 30.0 * W2
+                + C * (8.0 * W - 30.0 * W2)
+                - 3.0 * W
+                + 1.0
+            ) / 3.0
+            kernel = (q / 6.0 - 1.0 / 24.0) / (2.0 * np.pi)
+        else:
+            W2 = W * W
+            W3 = W2 * W
+            W4 = W3 * W
+            q = (
+                A * (70.0 * W4 - 60.0 * W3 + 6.0 * W2)
+                + 35.0 * W3 * (1.0 - C)
+                + C * 55.0 * W2 / 3.0
+                - 12.5 * W2
+                - W / 3.0
+                + 0.25
+            )
+            kernel = (q / 24.0 - 1.0 / 120.0) / (2.0 * np.pi)
+        tail = np.ones((values.shape[0], 1), dtype=np.float64)
+        constraint_tail = np.ones((setup_knots.shape[0], 1), dtype=np.float64)
+    else:
+        raise ValueError("Spherical spline order must be in {-2,-1,0,1,2,3,4}.")
+
+    return (
+        np.asarray(kernel, dtype=np.float64),
+        np.asarray(tail, dtype=np.float64),
+        np.asarray(constraint_tail, dtype=np.float64),
+    )
+
+
+def _parse_spherical_xt(xt):
+    max_knots = 2000
+    seed = 1
+    if xt is None:
+        return max_knots, seed
+    if not isinstance(xt, dict):
+        raise NotImplementedError(
+            "For bs='sos', xt must be None or a dict with optional keys "
+            "{'max.knots', 'seed'}."
+        )
+    if xt.get("max.knots") is not None:
+        max_knots = int(xt["max.knots"])
+    if xt.get("seed") is not None:
+        seed = int(xt["seed"])
+    if max_knots < 1:
+        raise ValueError("For bs='sos', xt['max.knots'] must be positive.")
+    return max_knots, seed
+
+
+@dataclass
+class SphericalSplineSetup:
+    knots: np.ndarray
+    UZ: np.ndarray
+    order: int
+    null_space_dim: int
+    rank: int
+    bs_dim: int
+    basis_train: np.ndarray
+    penalty: np.ndarray
+    column_scale: np.ndarray
+    used_supplied_knots: bool
+    used_subsampling: bool
+
+
+def build_spherical_spline_setup(X, *, k=-1, m=None, knots=None, xt=None):
+    """Port ``smooth.construct.sos.smooth.spec`` and retain prediction state."""
+    values = np.asarray(X, dtype=np.float64)
+    if values.ndim != 2 or values.shape[1] != 2 or values.shape[0] == 0:
+        raise ValueError(
+            "Can only deal with a sphere: bs='sos' requires exactly latitude and longitude."
+        )
+    if not np.isfinite(values).all():
+        raise ValueError("Spherical spline data must be finite.")
+
+    n_obs = int(values.shape[0])
+    order = normalize_spherical_order(m)
+    null_space_dim = spherical_null_space_dimension(order)
+    bs_dim = 50 if int(k) < 0 else int(k)
+    if bs_dim < null_space_dim + 2:
+        raise ValueError(
+            f"For bs='sos' with m={order}, k must be at least {null_space_dim + 2}."
+        )
+
+    supplied = _normalize_duchon_knots(knots, 2)
+    if supplied is not None and supplied.shape[0] == 0:
+        supplied = None
+    if supplied is not None and supplied.shape[0] > n_obs:
+        warnings.warn(
+            "more knots than data in an sos term: knots ignored.",
+            stacklevel=2,
+        )
+        supplied = None
+    max_knots, seed = _parse_spherical_xt(xt)
+    setup_knots, used_subsampling = _duchon_setup_locations(
+        values,
+        np.zeros(2, dtype=np.float64),
+        supplied,
+        max_knots=max_knots,
+        seed=seed,
+    )
+    n_knots = int(setup_knots.shape[0])
+    if bs_dim > n_knots:
+        raise ValueError(
+            "Spherical spline requires at least as many unique knot locations "
+            "as basis coefficients."
+        )
+
+    radial, _, constraint_tail = spherical_spline_kernel(
+        setup_knots, setup_knots, order
+    )
+    if bs_dim < n_knots:
+        eigenvalues, eigenvectors = _top_eigensystem(
+            radial,
+            bs_dim,
+            tolerance_exponent=0.5,
+        )
+        diagonal_penalty = np.diag(eigenvalues)
+        constraint = (constraint_tail.T @ eigenvectors).T
+    else:
+        eigenvectors = np.eye(bs_dim, dtype=np.float64)
+        diagonal_penalty = radial
+        constraint = constraint_tail
+
+    packed_qr, qraux = r_linpack_qr_no_pivot(constraint)
+    first = _r_linpack_qty(packed_qr, qraux, diagonal_penalty)
+    reduced = _r_linpack_qty(
+        packed_qr,
+        qraux,
+        first[null_space_dim:, :].T,
+    )[null_space_dim:, :]
+    penalty = np.zeros((bs_dim, bs_dim), dtype=np.float64)
+    rank = int(bs_dim - null_space_dim)
+    penalty[:rank, :rank] = reduced
+    UZ = _r_linpack_qty(
+        packed_qr,
+        qraux,
+        eigenvectors.T,
+    )[null_space_dim:, :].T
+
+    setup = SphericalSplineSetup(
+        knots=np.asarray(setup_knots, dtype=np.float64),
+        UZ=np.asarray(UZ, dtype=np.float64),
+        order=int(order),
+        null_space_dim=int(null_space_dim),
+        rank=rank,
+        bs_dim=int(bs_dim),
+        basis_train=np.zeros((n_obs, bs_dim), dtype=np.float64),
+        penalty=np.asarray(penalty, dtype=np.float64),
+        column_scale=np.ones(bs_dim, dtype=np.float64),
+        used_supplied_knots=bool(supplied is not None),
+        used_subsampling=bool(used_subsampling),
+    )
+    basis = predict_spherical_spline(values, setup, apply_scale=False)
+    standard_deviation = np.std(basis, axis=0, ddof=1)
+    standard_deviation[standard_deviation == np.min(standard_deviation)] = 1.0
+    column_scale = 1.0 / standard_deviation
+    setup.column_scale = np.asarray(column_scale, dtype=np.float64)
+    setup.basis_train = np.asarray(basis * column_scale[None, :], dtype=np.float64)
+    setup.penalty = np.asarray(
+        column_scale[:, None] * penalty * column_scale[None, :],
+        dtype=np.float64,
+    )
+    return setup
+
+
+def predict_spherical_spline(X_new, setup: SphericalSplineSetup, *, apply_scale=True):
+    """Port ``Predict.matrix.sos.smooth`` including knot-sized chunks."""
+    values = np.asarray(X_new, dtype=np.float64)
+    if values.ndim != 2 or values.shape[1] != 2:
+        raise ValueError(
+            "Spherical spline prediction requires latitude and longitude columns."
+        )
+    if not np.isfinite(values).all():
+        raise ValueError("Spherical spline prediction data must be finite.")
+    n_obs = int(values.shape[0])
+    n_knots = int(setup.knots.shape[0])
+    out = np.empty((n_obs, setup.bs_dim), dtype=np.float64)
+    for start in range(0, n_obs, n_knots):
+        stop = min(start + n_knots, n_obs)
+        radial, tail, _ = spherical_spline_kernel(
+            values[start:stop, :], setup.knots, setup.order
+        )
+        out[start:stop, :] = np.column_stack([radial @ setup.UZ, tail])
+    if apply_scale:
+        out *= setup.column_scale[None, :]
+    return np.asarray(out, dtype=np.float64)
+
+
+__all__ = [
+    "SphericalSplineSetup",
+    "build_spherical_spline_setup",
+    "normalize_spherical_order",
+    "predict_spherical_spline",
+    "spherical_null_space_dimension",
+    "spherical_spline_kernel",
+]
diff --git a/pyproject.toml b/pyproject.toml
index 8ebdcccf..9b8a0a45 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -210,6 +210,7 @@ markers = [
     "smooth_ds: tests covering Duchon regression spline smooths",
     "smooth_gp: tests covering Gaussian-process smooths",
     "smooth_mrf: tests covering Markov-random-field smooths",
+    "smooth_sos: tests covering splines on the sphere",
     "smooth_ps: tests covering P-spline smooths",
     "smooth_tp: tests covering thin plate smooths",
     "smooth_ts: tests covering shrinkage thin plate smooths",
diff --git a/tests/SUBSYSTEM_COVERAGE.md b/tests/SUBSYSTEM_COVERAGE.md
index 1fbfbb88..d73e239d 100644
--- a/tests/SUBSYSTEM_COVERAGE.md
+++ b/tests/SUBSYSTEM_COVERAGE.md
@@ -15,7 +15,7 @@ local development notes rather than duplicated here.
 | Subsystem | Primary owner(s) | Primary tests | Notes |
 | --- | --- | --- | --- |
 | Formula/spec parsing | `nampy/gam/formula/`, `nampy/gam/specs/` | `tests/parity/test_mgcv_formula_parse_parity.py` | Direct formula parity vs `mgcv`. |
-| Smooth constructors / raw basis owners | `nampy/gam/smooths/`, `nampy/gam/splines/` | `tests/smooths/test_mgcv_raw_constructor_parity.py`, `tests/smooths/test_mgcv_smoothcon_parity.py`, `tests/parity/test_gam_spec_build_owner_contracts.py`, `tests/parity/test_mgcv_cp_combinations_parity.py`, `tests/parity/test_mgcv_bs_combinations_parity.py`, `tests/parity/test_mgcv_ds_combinations_parity.py`, `tests/parity/test_mgcv_gp_combinations_parity.py`, `tests/parity/test_mgcv_mrf_combinations_parity.py` | Basis, penalty, constructor, and smoothCon surfaces, including cyclic P-splines (`cp`), integrated-derivative B-splines (`bs`), multivariate Duchon splines (`ds`), five-family stationary/nonstationary Gaussian-process smooths (`gp`), and graph/polygon/direct-penalty Markov random fields (`mrf`) across prediction, selection, linked bases, tensor/factor-smooth combinations, and explicit malformed-upstream boundaries. |
+| Smooth constructors / raw basis owners | `nampy/gam/smooths/`, `nampy/gam/splines/` | `tests/smooths/test_mgcv_raw_constructor_parity.py`, `tests/smooths/test_mgcv_smoothcon_parity.py`, `tests/parity/test_gam_spec_build_owner_contracts.py`, `tests/parity/test_mgcv_cp_combinations_parity.py`, `tests/parity/test_mgcv_bs_combinations_parity.py`, `tests/parity/test_mgcv_ds_combinations_parity.py`, `tests/parity/test_mgcv_gp_combinations_parity.py`, `tests/parity/test_mgcv_mrf_combinations_parity.py`, `tests/parity/test_mgcv_sos_combinations_parity.py` | Basis, penalty, constructor, and smoothCon surfaces, including cyclic P-splines (`cp`), integrated-derivative B-splines (`bs`), multivariate Duchon splines (`ds`), five-family stationary/nonstationary Gaussian-process smooths (`gp`), graph/polygon/direct-penalty Markov random fields (`mrf`), and all seven spherical-spline kernels (`sos`) across prediction, selection, linked bases, tensor/factor-smooth combinations, and explicit malformed-upstream boundaries. |
 | `pc=` / linked `id=` routing | smooth metadata + linked basis owners | `tests/smooths/test_mgcv_pc_id_parity.py` | Localizes shared-smoothing and point-constraint issues, including `te`/`ti` point constraints. |
 | Design / pre-fit assembly | `nampy/gam/compiler/`, `nampy/gam/fit/penalized_system.py` | `tests/optimization/test_mgcv_gam_setup_assembly_parity.py`, `tests/optimization/test_mgcv_preoptimization_blocks_parity.py`, `tests/optimization/test_mgcv_preoptimization_reparam_parity.py` | Setup, blocks, and reparameterization parity, including one global shared-component block with overlapping linear-predictor indices. |
 | Term wrapping / by-variable / offset routing | predictor wrapping + compiled term owners | `tests/optimization/test_gam_term_wrapping_owner_contracts.py` | Localizes wrapped predictor blocks, offset routing, and general-family block ownership before broader prediction parity. |
@@ -47,7 +47,7 @@ closes the former seven-stage under-tested backlog.
 | Pipeline stage | Combination coverage now owned |
 | --- | --- |
 | 1. Formula parsing and canonical specs | Formula lists/shared components, transformed covariates, supported numeric/factor interactions, factor-by smooths, formula and fit offsets, intercept policies, weights, knots, `min_sp`, `drop_intercept`, and tensor-`m` warning/fallback behavior; the combined interaction recipe is rebuilt on newdata and compared with a committed `mgcv` reference fixture. |
-| 2. Runtime terms and low-level bases | Train/newdata pairing for every supported basis (`bs/cr/cs/cc/cp/ds/gp/mrf/ps/tp/ts/re/fs/sz/te/ti`), including univariate, multivariate Duchon/GP, categorical MRF, and structured boundary response/SE parity. Row-permutation contracts cover linked and identified-FS terms, PS/TP/TS bases, `te`/`ti`, SZ, and factor-by smooths. Non-unique constructor representations use column-space projectors and penalized response operators instead of arbitrary coefficient orientation. |
+| 2. Runtime terms and low-level bases | Train/newdata pairing for every supported basis (`bs/cr/cs/cc/cp/ds/gp/mrf/ps/sos/tp/ts/re/fs/sz/te/ti`), including univariate, multivariate Duchon/GP, categorical MRF, spherical SOS, and structured boundary response/SE parity. Row-permutation contracts cover linked and identified-FS terms, PS/TP/TS bases, `te`/`ti`, SZ, and factor-by smooths. Non-unique constructor representations use column-space projectors and penalized response operators instead of arbitrary coefficient orientation. |
 | 3. Constructed/wrapped terms | Numeric-by plus linked `id=`, factor-by plus linked `id=`, tensor-by, and mixed fixed/free/select penalty ownership, numerical coefficient-map composition, and fitted response/SE parity. |
 | 4. Predictor/model compilation | Multi-predictor layouts with unequal intercept and offset policies, overlapping shared-component coefficient indices, and three-term linked, fixed/free, select, and rank-deficient assembly; supported `gaulss` and `gammals` layouts are also fitted and compared with `mgcv`. |
 | 5. Side conditions and identifiability | Repeated, nested, reverse-formula-order, tensor/main-effect, three-way, both identified near-rank regimes, zero-width, no-intercept, ordered factor-by, linked, general-family, two-predictor, SZ, and exempt random/factor smooth cases. Deletion rank and behavior replace raw pivot-column identity where QR/eigen choices are non-unique. |
diff --git a/tests/_taxonomy_registry.py b/tests/_taxonomy_registry.py
index 0fef056d..0d787166 100644
--- a/tests/_taxonomy_registry.py
+++ b/tests/_taxonomy_registry.py
@@ -22,6 +22,7 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation:
     "ds": "smooth_ds",
     "gp": "smooth_gp",
     "mrf": "smooth_mrf",
+    "sos": "smooth_sos",
     "ps": "smooth_ps",
     "tp": "smooth_tp",
     "ts": "smooth_ts",
@@ -89,6 +90,7 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation:
     "test_mgcv_raw_constructor_parity.py": {"surface_smoothcon"},
     "test_mgcv_ds_combinations_parity.py": {"surface_snapshot"},
     "test_mgcv_mrf_combinations_parity.py": {"surface_snapshot"},
+    "test_mgcv_sos_combinations_parity.py": {"surface_snapshot"},
     "test_mgcv_score_hist_trace_parity.py": {"surface_trace"},
     "test_mgcv_optimization_lifecycle_parity.py": {"surface_trace"},
     "test_mgcv_linked_id_trace_parity.py": {"surface_trace"},
@@ -125,6 +127,7 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation:
     "test_mgcv_output_parity.py",
     "test_mgcv_ds_combinations_parity.py",
     "test_mgcv_mrf_combinations_parity.py",
+    "test_mgcv_sos_combinations_parity.py",
     "test_mgcv_score_hist_trace_parity.py",
     "test_mgcv_linked_id_trace_parity.py",
     "test_mgcv_score_gamma_parity.py",
@@ -156,6 +159,11 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation:
         "tests/parity/test_mgcv_mrf_combinations_parity.py",
         "tests/smooths/test_mgcv_raw_constructor_parity.py",
     ),
+    "smooth_sos": (
+        "tests/parity/test_mgcv_sos_combinations_parity.py",
+        "tests/smooths/test_mgcv_raw_constructor_parity.py",
+        "tests/smooths/test_mgcv_smoothcon_parity.py",
+    ),
     "smooth_cs": (
         "tests/smooths/test_mgcv_pc_id_parity.py",
         "tests/smooths/test_mgcv_raw_constructor_parity.py",
diff --git a/tests/diagnostics/test_mgcv_k_check_parity.py b/tests/diagnostics/test_mgcv_k_check_parity.py
index d1900f3d..99dcb01c 100644
--- a/tests/diagnostics/test_mgcv_k_check_parity.py
+++ b/tests/diagnostics/test_mgcv_k_check_parity.py
@@ -13,6 +13,7 @@
 from __future__ import annotations
 
 import numpy as np
+import pandas as pd
 import pytest
 
 from tests.mgcv_parity_utils import (
@@ -231,6 +232,15 @@ def _assert_k_check_parity(
 # --------------------------------------------------------------------------- #
 
 
+def _make_sos_kcheck_data(seed=603, n=180):
+    rng = np.random.default_rng(seed)
+    lo = rng.uniform(-180.0, 180.0, size=n)
+    la = np.rad2deg(np.arcsin(rng.uniform(-1.0, 1.0, size=n)))
+    y = np.sin(np.deg2rad(lo)) * np.cos(np.deg2rad(la - 12.0))
+    y += rng.normal(scale=0.12, size=n)
+    return pd.DataFrame({"y": y, "la": la, "lo": lo})
+
+
 class TestKCheckParity:
     """Compare k_check() output against mgcv::k.check() for each smooth type."""
 
@@ -277,6 +287,14 @@ class TestKCheckParity:
                 {"x0", "x1"},
                 1e-4,
             ),
+            (
+                _make_sos_kcheck_data,
+                'y ~ s(la, lo, bs="sos", k=12)',
+                "gaussian",
+                "REML",
+                {"la", "lo"},
+                1e-4,
+            ),
             (
                 lambda: _make_gaussian_data(seed=123, n=180),
                 'y ~ te(x0, x1, bs=["cr","cr"], k=[5,5])',
@@ -292,6 +310,7 @@ class TestKCheckParity:
             "gaussian_ps",
             "gaussian_cp",
             "gaussian_bs",
+            "gaussian_sos",
             "gaussian_te",
         ],
     )
diff --git a/tests/mgcv_invariant_policy.py b/tests/mgcv_invariant_policy.py
index 0643ee82..5b32e5c2 100644
--- a/tests/mgcv_invariant_policy.py
+++ b/tests/mgcv_invariant_policy.py
@@ -320,6 +320,16 @@ def _canonicalize_gp_raw_state(state):
     return state
 
 
+def _canonicalize_sos_raw_state(state):
+    extra = state["extra"]
+    extra.pop("used_supplied_knots", False)
+    extra.pop("used_subsampling", False)
+    state["S"] = [penalty_spectrum(S) for S in state["S"]]
+    state["X"] = matrix_self_gram(state["X"])
+    extra["UZ"] = stable_column_space_projector(extra["UZ"])
+    return state
+
+
 def _canonicalize_mrf_raw_state(state):
     extra = state["extra"]
     if extra.get("P") is not None:
@@ -389,6 +399,8 @@ def canonicalize_raw_representation_state(state: dict[str, Any]) -> dict[str, An
         return _canonicalize_duchon_raw_state(state)
     if class_name == "gp.smooth":
         return _canonicalize_gp_raw_state(state)
+    if class_name == "sos.smooth":
+        return _canonicalize_sos_raw_state(state)
     if class_name == "mrf.smooth":
         return _canonicalize_mrf_raw_state(state)
     if class_name == "fs.interaction":
diff --git a/tests/mgcv_parity_utils.py b/tests/mgcv_parity_utils.py
index dafc3353..b60530ff 100644
--- a/tests/mgcv_parity_utils.py
+++ b/tests/mgcv_parity_utils.py
@@ -1651,6 +1651,12 @@ def _run_mgcv_raw_constructor(
       shift = pack_vector(sm$shift, "numeric"),
       gp_defn = pack_vector(sm$gp.defn, "numeric")
     ),
+    "sos.smooth" = list(
+      knt = pack_vector(sm$knt, "numeric"),
+      UZ = pack_matrix(sm$UZ),
+      p_order = as.integer(sm$p.order),
+      xc_scale = pack_vector(sm$xc.scale, "numeric")
+    ),
     "mrf.smooth" = list(
       P = pack_matrix(sm$P),
       knots = pack_vector(sm$knots, "character"),
diff --git a/tests/parity/test_mgcv_sos_combinations_parity.py b/tests/parity/test_mgcv_sos_combinations_parity.py
new file mode 100644
index 00000000..0ebb8113
--- /dev/null
+++ b/tests/parity/test_mgcv_sos_combinations_parity.py
@@ -0,0 +1,253 @@
+"""Integrated parity coverage for spherical splines (``bs='sos'``)."""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from nampy.gam import GAM
+from nampy.gam.splines.univariate.sos import (
+    build_spherical_spline_setup,
+    predict_spherical_spline,
+)
+from tests.mgcv_parity_utils import (
+    _assert_basic_mgcv_parity,
+    _fit_nampy_model,
+    _fit_nampy_snapshot,
+    _run_mgcv_snapshot,
+)
+
+
+def _sos_data(seed=941, n=180):
+    rng = np.random.default_rng(seed)
+    lo = rng.uniform(-180.0, 180.0, size=n)
+    la = np.rad2deg(np.arcsin(rng.uniform(-1.0, 1.0, size=n)))
+    x = rng.uniform(-1.5, 1.5, size=n)
+    z = 0.8 + rng.uniform(-0.3, 0.5, size=n)
+    g = np.asarray(["a", "b", "c"], dtype=object)[np.arange(n) % 3]
+    y = (
+        0.2
+        + z * np.sin(np.deg2rad(lo)) * np.cos(np.deg2rad(la - 12.0))
+        + 0.2 * x
+        + 0.15 * (g == "b")
+        - 0.1 * (g == "c")
+        + rng.normal(scale=0.12, size=n)
+    )
+    return pd.DataFrame({"y": y, "la": la, "lo": lo, "x": x, "z": z, "g": g})
+
+
+def _assert_snapshot_fit(actual, expected, *, atol=1e-6):
+    for key in ("response", "link"):
+        np.testing.assert_allclose(
+            actual["predictions"][key],
+            expected["predictions"][key],
+            atol=atol,
+            rtol=atol,
+        )
+    np.testing.assert_allclose(
+        actual["fit"]["edf_total"],
+        expected["fit"]["edf_total"],
+        atol=atol,
+        rtol=atol,
+    )
+
+
+def test_sos_default_reml_matches_mgcv():
+    data = _sos_data(seed=942)
+    formula = 'y ~ s(la, lo, bs="sos", k=12)'
+    _assert_basic_mgcv_parity(
+        _fit_nampy_snapshot(data, formula, "gaussian", "REML"),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML"),
+        pred_atol=1e-6,
+        pred_rtol=1e-6,
+        sp_log_atol=2e-5,
+        criterion_atol=2e-6,
+    )
+
+
+def test_sos_duchon_null_space_select_matches_mgcv():
+    data = _sos_data(seed=943)
+    formula = 'y ~ s(la, lo, bs="sos", k=12, m=-1)'
+    model = _fit_nampy_model(data, formula, "gaussian", "REML", select=True)
+    runtime = [
+        term.predict_fn.__self__
+        for term in model.gam_result_.compiled_model.compiled_terms
+        if term.term_type == "smooth"
+    ][0]
+    assert runtime._setup.null_space_dim == 4
+    assert len(model.smoothing_params) == 2
+    _assert_basic_mgcv_parity(
+        model.parity_snapshot(X=data, include_covariances=True),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML", select=True),
+        pred_atol=1e-6,
+        pred_rtol=1e-6,
+        sp_log_atol=3e-5,
+        criterion_atol=3e-6,
+    )
+
+
+def test_sos_numeric_and_factor_by_fixed_sp_match_mgcv():
+    data = _sos_data(seed=944)
+    numeric = 'y ~ s(la, lo, by=z, bs="sos", k=12, m=2, sp=.7)'
+    _assert_snapshot_fit(
+        _fit_nampy_snapshot(data, numeric, "gaussian", "fixed"),
+        _run_mgcv_snapshot(data, numeric, "gaussian", "REML"),
+    )
+
+    factor = 'y ~ s(la, lo, by=g, bs="sos", k=12, m=-2, sp=.7)'
+    _assert_snapshot_fit(
+        _fit_nampy_snapshot(data, factor, "gaussian", "fixed"),
+        _run_mgcv_snapshot(data, factor, "gaussian", "REML"),
+        atol=2e-6,
+    )
+
+
+def test_sos_linked_id_and_point_constraint_match_mgcv():
+    data = _sos_data(seed=945, n=150)
+    data["la2"] = np.roll(data["la"].to_numpy(), 7)
+    data["lo2"] = np.roll(data["lo"].to_numpy(), 11)
+    linked = (
+        'y ~ s(la, lo, bs="sos", k=12, m=3, id="sphere", sp=.7)'
+        ' + s(la2, lo2, bs="sos", k=14, m=-1, id="sphere")'
+    )
+    _assert_snapshot_fit(
+        _fit_nampy_snapshot(data, linked, "gaussian", "fixed"),
+        _run_mgcv_snapshot(data, linked, "gaussian", "REML"),
+        atol=2e-6,
+    )
+
+    point = 'y ~ s(la, lo, bs="sos", k=12, pc=[0,0], sp=.7)'
+    _assert_snapshot_fit(
+        _fit_nampy_snapshot(data, point, "gaussian", "fixed"),
+        _run_mgcv_snapshot(data, point, "gaussian", "REML"),
+    )
+
+
+def test_sos_dynamic_null_penalties_link_and_fs_boundary():
+    data = _sos_data(seed=955, n=150)
+    data["la2"] = np.roll(data["la"].to_numpy(), 5)
+    data["lo2"] = np.roll(data["lo"].to_numpy(), 9)
+    linked = (
+        'y ~ s(la, lo, bs="sos", k=12, m=-1, id="sphere_null", sp=[.7,.9])'
+        ' + s(la2, lo2, bs="sos", k=14, m=3, id="sphere_null")'
+    )
+    _assert_snapshot_fit(
+        _fit_nampy_snapshot(data, linked, "gaussian", "fixed", select=True),
+        _run_mgcv_snapshot(data, linked, "gaussian", "REML", select=True),
+        atol=3e-6,
+    )
+
+    factor = 'y ~ s(g, la, lo, bs="fs", xt="sos", k=8, m=-1, sp=[.5,.6,.7,.8,.9])'
+    with pytest.raises(NotImplementedError, match="LAPACK-dependent"):
+        _fit_nampy_snapshot(data, factor, "gaussian", "fixed")
+
+
+@pytest.mark.parametrize(
+    "formula",
+    [
+        'y ~ te(la, lo, x, d=[2,1], bs=["sos","cr"], k=[8,4], sp=[.6,.8])',
+        'y ~ ti(la, lo, x, d=[2,1], bs=["sos","cr"], k=[8,4], sp=[.6,.8])',
+    ],
+    ids=["te", "ti"],
+)
+def test_sos_tensor_margin_fixed_sp_matches_mgcv(formula):
+    data = _sos_data(seed=946, n=150)
+    _assert_snapshot_fit(
+        _fit_nampy_snapshot(data, formula, "gaussian", "fixed"),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML"),
+        atol=3e-6,
+    )
+
+
+@pytest.mark.parametrize(
+    "formula",
+    [
+        'y ~ s(g, la, lo, bs="fs", xt="sos", k=8, sp=[.6,.8])',
+        'y ~ s(g, la, lo, bs="sz", xt="sos", k=8, sp=[.7,.7,.7])',
+    ],
+    ids=["fs", "sz"],
+)
+def test_sos_factor_smooth_base_fixed_sp_matches_mgcv(formula):
+    data = _sos_data(seed=947, n=150)
+    _assert_snapshot_fit(
+        _fit_nampy_snapshot(data, formula, "gaussian", "fixed"),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML"),
+        atol=3e-6,
+    )
+
+
+def test_sos_prediction_periodicity_poles_and_persistence(tmp_path):
+    data = _sos_data(seed=948, n=130)
+    formula = 'y ~ s(la, lo, bs="sos", k=12, m=-1, sp=.7)'
+    model = _fit_nampy_model(data, formula, "gaussian", "fixed")
+    newdata = pd.DataFrame(
+        {
+            "la": [-90.0, -45.0, 0.0, 45.0, 90.0],
+            "lo": [-170.0, -80.0, 10.0, 100.0, 170.0],
+        }
+    )
+    wrapped = newdata.copy()
+    wrapped["lo"] += 360.0
+    expected = model.predict(newdata, type="link", block_size=1)
+    np.testing.assert_allclose(
+        model.predict(wrapped, type="link", block_size=1), expected, atol=2e-13
+    )
+    poles = pd.DataFrame({"la": [-90.0] * 3 + [90.0] * 3, "lo": [-120, 0, 120] * 2})
+    pole_predictions = model.predict(poles, type="link")
+    np.testing.assert_allclose(pole_predictions[:3], pole_predictions[0], atol=2e-13)
+    np.testing.assert_allclose(pole_predictions[3:], pole_predictions[3], atol=2e-13)
+
+    path = tmp_path / "sos.pkl"
+    model.save_model(path)
+    restored = GAM.load_model(path)
+    np.testing.assert_allclose(restored.predict(newdata, type="link"), expected)
+
+
+def test_sos_order_knots_and_dimension_guards():
+    data = _sos_data(seed=949, n=60)
+    X = data[["la", "lo"]].to_numpy()
+    assert build_spherical_spline_setup(X, k=12, m=-3).order == -1
+    assert build_spherical_spline_setup(X, k=12, m=1.6).order == 2
+    assert build_spherical_spline_setup(X, k=12, m=5).order == 4
+    assert build_spherical_spline_setup(X, k=12, xt={"ignored": True}).bs_dim == 12
+    with pytest.warns(UserWarning, match="more knots than data in an sos term"):
+        ignored = build_spherical_spline_setup(
+            X, k=12, knots=[np.arange(61), np.arange(61)]
+        )
+    assert not ignored.used_supplied_knots
+    with pytest.raises(ValueError, match="at least 6"):
+        build_spherical_spline_setup(X, k=5, m=-1)
+    with pytest.raises(ValueError, match="at least as many unique knot locations"):
+        build_spherical_spline_setup(X, k=12, knots=[np.arange(10), np.arange(10)])
+    with pytest.raises(ValueError, match="single numeric value"):
+        build_spherical_spline_setup(X, k=12, m=[1, 2])
+
+    setup = build_spherical_spline_setup(X, k=12, m=0)
+    np.testing.assert_allclose(
+        predict_spherical_spline(X, setup), setup.basis_train, atol=2e-13
+    )
+    np.testing.assert_allclose(
+        setup.UZ.T @ np.ones(setup.knots.shape[0]), 0.0, atol=2e-12
+    )
+
+
+def test_sos_array_and_derivative_surfaces_are_explicitly_unsupported():
+    data = _sos_data(seed=950, n=80)
+    with pytest.raises(NotImplementedError, match="formula with both latitude"):
+        GAM(
+            family="gaussian",
+            basis="sos",
+            k=12,
+            optimize_smoothing=False,
+            smoothing_params=[0.7],
+        ).fit(X=data[["la"]], y=data["y"].to_numpy())
+
+    model = GAM(
+        formula='y ~ s(la, lo, bs="sos", k=12, sp=.7)',
+        optimize_smoothing=False,
+    ).fit(data=data)
+    with pytest.raises(NotImplementedError, match="only 1D smooths"):
+        model.derivative(data, smooth_number=1)
+    with pytest.raises(NotImplementedError, match="rotated hemisphere"):
+        model.plot()
diff --git a/tests/reference_fixtures/mgcv/010ff4c0a97943e32a1bc6840ad52cc732524a95d0295c777caf97e68d177ec7.json.gz b/tests/reference_fixtures/mgcv/010ff4c0a97943e32a1bc6840ad52cc732524a95d0295c777caf97e68d177ec7.json.gz
new file mode 100644
index 00000000..282e0ab2
Binary files /dev/null and b/tests/reference_fixtures/mgcv/010ff4c0a97943e32a1bc6840ad52cc732524a95d0295c777caf97e68d177ec7.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/052eae2874040cd8ec986cc7d051194159ad47ebeb129b03fa0f882778ca9c1e.json.gz b/tests/reference_fixtures/mgcv/052eae2874040cd8ec986cc7d051194159ad47ebeb129b03fa0f882778ca9c1e.json.gz
new file mode 100644
index 00000000..2c08375a
Binary files /dev/null and b/tests/reference_fixtures/mgcv/052eae2874040cd8ec986cc7d051194159ad47ebeb129b03fa0f882778ca9c1e.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/0bb9bc9f1923e289e226c6878abd5912ea7688388f95506ebc48442ded9c0c21.json.gz b/tests/reference_fixtures/mgcv/0bb9bc9f1923e289e226c6878abd5912ea7688388f95506ebc48442ded9c0c21.json.gz
new file mode 100644
index 00000000..31edd080
Binary files /dev/null and b/tests/reference_fixtures/mgcv/0bb9bc9f1923e289e226c6878abd5912ea7688388f95506ebc48442ded9c0c21.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/1291e6907564ad8fab20bd4a7622d054ba79a2a7830dda9f004578fa172a7a98.json.gz b/tests/reference_fixtures/mgcv/1291e6907564ad8fab20bd4a7622d054ba79a2a7830dda9f004578fa172a7a98.json.gz
new file mode 100644
index 00000000..1cc9b001
Binary files /dev/null and b/tests/reference_fixtures/mgcv/1291e6907564ad8fab20bd4a7622d054ba79a2a7830dda9f004578fa172a7a98.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/130e6356f5906cf2a274dba2d840d3af2e1634cf46407310cefe78d7de3f111b.json.gz b/tests/reference_fixtures/mgcv/130e6356f5906cf2a274dba2d840d3af2e1634cf46407310cefe78d7de3f111b.json.gz
new file mode 100644
index 00000000..5652a0e3
Binary files /dev/null and b/tests/reference_fixtures/mgcv/130e6356f5906cf2a274dba2d840d3af2e1634cf46407310cefe78d7de3f111b.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/28848b22a57366cb51cbc74c50f2e276a0db39075a517fb3d6b74540be0e656b.json.gz b/tests/reference_fixtures/mgcv/28848b22a57366cb51cbc74c50f2e276a0db39075a517fb3d6b74540be0e656b.json.gz
new file mode 100644
index 00000000..21749bf9
Binary files /dev/null and b/tests/reference_fixtures/mgcv/28848b22a57366cb51cbc74c50f2e276a0db39075a517fb3d6b74540be0e656b.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/29e27e156572fb3cbe03fc94fd013cecb2b08c6bdf9407a69f95b32f222d1b41.json.gz b/tests/reference_fixtures/mgcv/29e27e156572fb3cbe03fc94fd013cecb2b08c6bdf9407a69f95b32f222d1b41.json.gz
new file mode 100644
index 00000000..11914194
Binary files /dev/null and b/tests/reference_fixtures/mgcv/29e27e156572fb3cbe03fc94fd013cecb2b08c6bdf9407a69f95b32f222d1b41.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/2ad8b7db3a1313ab7ffb9ce5721aa29c22842148dd3d7341836405e20aebf597.json.gz b/tests/reference_fixtures/mgcv/2ad8b7db3a1313ab7ffb9ce5721aa29c22842148dd3d7341836405e20aebf597.json.gz
new file mode 100644
index 00000000..2f121ff8
Binary files /dev/null and b/tests/reference_fixtures/mgcv/2ad8b7db3a1313ab7ffb9ce5721aa29c22842148dd3d7341836405e20aebf597.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/32d878c7fbd1ecd42ab8ff6a637fb6499506d6220a4b03a0b48d4065b54d3f1d.json.gz b/tests/reference_fixtures/mgcv/32d878c7fbd1ecd42ab8ff6a637fb6499506d6220a4b03a0b48d4065b54d3f1d.json.gz
new file mode 100644
index 00000000..4f3ffdc7
Binary files /dev/null and b/tests/reference_fixtures/mgcv/32d878c7fbd1ecd42ab8ff6a637fb6499506d6220a4b03a0b48d4065b54d3f1d.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/3f6080bc733831c0fb434e4ecab8b0d908ce24bd8e98ad53ba23eacb37c397b5.json.gz b/tests/reference_fixtures/mgcv/3f6080bc733831c0fb434e4ecab8b0d908ce24bd8e98ad53ba23eacb37c397b5.json.gz
new file mode 100644
index 00000000..415bcee1
Binary files /dev/null and b/tests/reference_fixtures/mgcv/3f6080bc733831c0fb434e4ecab8b0d908ce24bd8e98ad53ba23eacb37c397b5.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/428a284d995d4dd74f5cf260e083822c44d0eff5ca68f66b4d0e0016bc66cc6f.json.gz b/tests/reference_fixtures/mgcv/428a284d995d4dd74f5cf260e083822c44d0eff5ca68f66b4d0e0016bc66cc6f.json.gz
new file mode 100644
index 00000000..2d9b0c7d
Binary files /dev/null and b/tests/reference_fixtures/mgcv/428a284d995d4dd74f5cf260e083822c44d0eff5ca68f66b4d0e0016bc66cc6f.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/4941bb9846e5f3ce6fbc520e8faa021e339ec25cf8940126998057e692c91dc0.json.gz b/tests/reference_fixtures/mgcv/4941bb9846e5f3ce6fbc520e8faa021e339ec25cf8940126998057e692c91dc0.json.gz
new file mode 100644
index 00000000..d5ffcc45
Binary files /dev/null and b/tests/reference_fixtures/mgcv/4941bb9846e5f3ce6fbc520e8faa021e339ec25cf8940126998057e692c91dc0.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/622d0b9578ac311d6688702e51e61f95aeff4622a7c1bf249362d6781638a717.json.gz b/tests/reference_fixtures/mgcv/622d0b9578ac311d6688702e51e61f95aeff4622a7c1bf249362d6781638a717.json.gz
new file mode 100644
index 00000000..6200643b
Binary files /dev/null and b/tests/reference_fixtures/mgcv/622d0b9578ac311d6688702e51e61f95aeff4622a7c1bf249362d6781638a717.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/636bd9a820e0dbe7799ae16cb3c8b4aaea7998ccc7e5076f22c4b23dfc98140e.json.gz b/tests/reference_fixtures/mgcv/636bd9a820e0dbe7799ae16cb3c8b4aaea7998ccc7e5076f22c4b23dfc98140e.json.gz
new file mode 100644
index 00000000..6213592e
Binary files /dev/null and b/tests/reference_fixtures/mgcv/636bd9a820e0dbe7799ae16cb3c8b4aaea7998ccc7e5076f22c4b23dfc98140e.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/69b2840ca487dca4c4d31907fdd773891a704c32c7f723b568bba78c07f336df.json.gz b/tests/reference_fixtures/mgcv/69b2840ca487dca4c4d31907fdd773891a704c32c7f723b568bba78c07f336df.json.gz
new file mode 100644
index 00000000..fe2be4b5
Binary files /dev/null and b/tests/reference_fixtures/mgcv/69b2840ca487dca4c4d31907fdd773891a704c32c7f723b568bba78c07f336df.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/6bdfda1ab76f68aa8a3fea340a11d2371a8ba7bc72ca9ee02f7304d9875c3d5f.json.gz b/tests/reference_fixtures/mgcv/6bdfda1ab76f68aa8a3fea340a11d2371a8ba7bc72ca9ee02f7304d9875c3d5f.json.gz
new file mode 100644
index 00000000..095bdc6c
Binary files /dev/null and b/tests/reference_fixtures/mgcv/6bdfda1ab76f68aa8a3fea340a11d2371a8ba7bc72ca9ee02f7304d9875c3d5f.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/79d204be16e17f784762534a3fac4a980dd48acedd7b96a74a6a7d759a1e3ad0.json.gz b/tests/reference_fixtures/mgcv/79d204be16e17f784762534a3fac4a980dd48acedd7b96a74a6a7d759a1e3ad0.json.gz
new file mode 100644
index 00000000..1f29b2cf
Binary files /dev/null and b/tests/reference_fixtures/mgcv/79d204be16e17f784762534a3fac4a980dd48acedd7b96a74a6a7d759a1e3ad0.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/83fc300c9d5c40e21bc19226ed1825e3c898c72bcd687d70816eb83beafc5de4.json.gz b/tests/reference_fixtures/mgcv/83fc300c9d5c40e21bc19226ed1825e3c898c72bcd687d70816eb83beafc5de4.json.gz
new file mode 100644
index 00000000..df68b357
Binary files /dev/null and b/tests/reference_fixtures/mgcv/83fc300c9d5c40e21bc19226ed1825e3c898c72bcd687d70816eb83beafc5de4.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/89fc90c7327dd7cd0d5b16fad6d504edb58985ab11dedbcd7f810d847abf6749.json.gz b/tests/reference_fixtures/mgcv/89fc90c7327dd7cd0d5b16fad6d504edb58985ab11dedbcd7f810d847abf6749.json.gz
new file mode 100644
index 00000000..56ca7084
Binary files /dev/null and b/tests/reference_fixtures/mgcv/89fc90c7327dd7cd0d5b16fad6d504edb58985ab11dedbcd7f810d847abf6749.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/9319e3210aed8657ed65c30fe0cc31695dcf474a78b98a6b1ed1c23b634d1acd.json.gz b/tests/reference_fixtures/mgcv/9319e3210aed8657ed65c30fe0cc31695dcf474a78b98a6b1ed1c23b634d1acd.json.gz
new file mode 100644
index 00000000..3418cbd5
Binary files /dev/null and b/tests/reference_fixtures/mgcv/9319e3210aed8657ed65c30fe0cc31695dcf474a78b98a6b1ed1c23b634d1acd.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/9ad52f7b90dde69108755c78a0544a99b3faf21de7e595aa3dff5c53bbb8e3f5.json.gz b/tests/reference_fixtures/mgcv/9ad52f7b90dde69108755c78a0544a99b3faf21de7e595aa3dff5c53bbb8e3f5.json.gz
new file mode 100644
index 00000000..d1e70429
Binary files /dev/null and b/tests/reference_fixtures/mgcv/9ad52f7b90dde69108755c78a0544a99b3faf21de7e595aa3dff5c53bbb8e3f5.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/9fcf5a735b3f2c36f0ff5db1a6132df2e59b4a25287282ec3fc13ed791265570.json.gz b/tests/reference_fixtures/mgcv/9fcf5a735b3f2c36f0ff5db1a6132df2e59b4a25287282ec3fc13ed791265570.json.gz
new file mode 100644
index 00000000..c9f8ba13
Binary files /dev/null and b/tests/reference_fixtures/mgcv/9fcf5a735b3f2c36f0ff5db1a6132df2e59b4a25287282ec3fc13ed791265570.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/a787afb1de97b549a3d1de25a438f8fe6617544a4aba29508a410fda8fa6a936.json.gz b/tests/reference_fixtures/mgcv/a787afb1de97b549a3d1de25a438f8fe6617544a4aba29508a410fda8fa6a936.json.gz
new file mode 100644
index 00000000..69065a70
Binary files /dev/null and b/tests/reference_fixtures/mgcv/a787afb1de97b549a3d1de25a438f8fe6617544a4aba29508a410fda8fa6a936.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/b5c8654cbd979b79d365c9a2c4723735e942b59736295a0c624fd08d08de32c1.json.gz b/tests/reference_fixtures/mgcv/b5c8654cbd979b79d365c9a2c4723735e942b59736295a0c624fd08d08de32c1.json.gz
new file mode 100644
index 00000000..63a645e0
Binary files /dev/null and b/tests/reference_fixtures/mgcv/b5c8654cbd979b79d365c9a2c4723735e942b59736295a0c624fd08d08de32c1.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/be77e69a95933662f66cb2f23b1b2124fff85cc44ad6bc5b318e60d68d294cd0.json.gz b/tests/reference_fixtures/mgcv/be77e69a95933662f66cb2f23b1b2124fff85cc44ad6bc5b318e60d68d294cd0.json.gz
new file mode 100644
index 00000000..70f14662
Binary files /dev/null and b/tests/reference_fixtures/mgcv/be77e69a95933662f66cb2f23b1b2124fff85cc44ad6bc5b318e60d68d294cd0.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/d1306dbd3e0a16d73fc6c3dd28778909cf1a8fd73bfd21b1e65d1b4a80a7cc15.json.gz b/tests/reference_fixtures/mgcv/d1306dbd3e0a16d73fc6c3dd28778909cf1a8fd73bfd21b1e65d1b4a80a7cc15.json.gz
new file mode 100644
index 00000000..50b12a5f
Binary files /dev/null and b/tests/reference_fixtures/mgcv/d1306dbd3e0a16d73fc6c3dd28778909cf1a8fd73bfd21b1e65d1b4a80a7cc15.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/e8ab6048c7725e5ed4239080579513677a508d1320bf04bddc0d894c3d13082f.json.gz b/tests/reference_fixtures/mgcv/e8ab6048c7725e5ed4239080579513677a508d1320bf04bddc0d894c3d13082f.json.gz
new file mode 100644
index 00000000..1862d496
Binary files /dev/null and b/tests/reference_fixtures/mgcv/e8ab6048c7725e5ed4239080579513677a508d1320bf04bddc0d894c3d13082f.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/f1e5a90f8c8a9dd7cef5c43e69fb888e127acd7266afd788b06523d6ea3f81fe.json.gz b/tests/reference_fixtures/mgcv/f1e5a90f8c8a9dd7cef5c43e69fb888e127acd7266afd788b06523d6ea3f81fe.json.gz
new file mode 100644
index 00000000..5b4cd70e
Binary files /dev/null and b/tests/reference_fixtures/mgcv/f1e5a90f8c8a9dd7cef5c43e69fb888e127acd7266afd788b06523d6ea3f81fe.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/f30d10dc8481a7999c1c27940fb29cd52ccea6d3498a60bd7d393ab9f1f8e152.json.gz b/tests/reference_fixtures/mgcv/f30d10dc8481a7999c1c27940fb29cd52ccea6d3498a60bd7d393ab9f1f8e152.json.gz
new file mode 100644
index 00000000..591ae183
Binary files /dev/null and b/tests/reference_fixtures/mgcv/f30d10dc8481a7999c1c27940fb29cd52ccea6d3498a60bd7d393ab9f1f8e152.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/fec8f5fb7f7907b1d0547bf4a37fa77b90d9425240363cbfdee37d09a9f9810e.json.gz b/tests/reference_fixtures/mgcv/fec8f5fb7f7907b1d0547bf4a37fa77b90d9425240363cbfdee37d09a9f9810e.json.gz
new file mode 100644
index 00000000..57015b77
Binary files /dev/null and b/tests/reference_fixtures/mgcv/fec8f5fb7f7907b1d0547bf4a37fa77b90d9425240363cbfdee37d09a9f9810e.json.gz differ
diff --git a/tests/smooths/test_mgcv_raw_constructor_parity.py b/tests/smooths/test_mgcv_raw_constructor_parity.py
index 262de4db..ab6cf0f8 100644
--- a/tests/smooths/test_mgcv_raw_constructor_parity.py
+++ b/tests/smooths/test_mgcv_raw_constructor_parity.py
@@ -32,6 +32,7 @@
 from nampy.gam.smooths.univariate.ds import DuchonSplineTerm
 from nampy.gam.smooths.univariate.gp import GaussianProcessTerm
 from nampy.gam.smooths.univariate.ps import PSplineTerm1D
+from nampy.gam.smooths.univariate.sos import SphericalSplineTerm
 from nampy.gam.smooths.univariate.tp import ThinPlateSplineTerm
 from nampy.gam.specs.build import build_formula_model
 from nampy.gam.splines.basis.cr import cr_exact_null_basis_from_knots
@@ -150,6 +151,17 @@ def _make_numeric_mrf_factor_data(seed=906, n=100):
     return pd.DataFrame({"y": y, "region": region, "group": group})
 
 
+def _make_spherical_data(seed=911, n=100):
+    rng = np.random.default_rng(seed)
+    longitude = rng.uniform(-180.0, 180.0, size=n)
+    latitude = np.rad2deg(np.arcsin(rng.uniform(-1.0, 1.0, size=n)))
+    y = (
+        np.sin(np.deg2rad(longitude)) * np.cos(np.deg2rad(latitude - 15.0))
+        + rng.normal(scale=0.1, size=n)
+    )
+    return pd.DataFrame({"y": y, "la": latitude, "lo": longitude})
+
+
 def _mrf_knots_with_unobserved(_data):
     levels = ["a", "b", "c", "d", "e"]
     return {"region": pd.Categorical(levels, categories=levels)}
@@ -827,6 +839,49 @@ def _build_gp_case_matrix():
     ]
 
 
+def _build_sos_case_matrix():
+    cases = [
+        _case(
+            f"sos_order_{order}",
+            _factory(_make_spherical_data, seed=920 + order, n=90),
+            f'y ~ s(la, lo, bs="sos", k=12, m={order})',
+            atol=2e-7,
+        )
+        for order in (-2, -1, 0, 1, 2, 3, 4)
+    ]
+    cases.extend(
+        [
+            _case(
+                "sos_default_k",
+                _factory(_make_spherical_data, seed=930, n=90),
+                'y ~ s(la, lo, bs="sos")',
+                atol=5e-7,
+            ),
+            _case(
+                "sos_supplied_truncated",
+                _factory(_make_spherical_data, seed=931, n=90),
+                'y ~ s(la, lo, bs="sos", k=12, m=2)',
+                atol=2e-7,
+                knots_factory=_observed_row_knots(["la", "lo"], 20),
+            ),
+            _case(
+                "sos_supplied_full",
+                _factory(_make_spherical_data, seed=932, n=90),
+                'y ~ s(la, lo, bs="sos", k=12, m=-1)',
+                atol=2e-7,
+                knots_factory=_observed_row_knots(["la", "lo"], 12),
+            ),
+            _case(
+                "sos_max_knots_xt",
+                _factory(_make_spherical_data, seed=933, n=100),
+                'y ~ s(la, lo, bs="sos", k=12, xt={"max.knots":25,"seed":7})',
+                atol=2e-7,
+            ),
+        ]
+    )
+    return cases
+
+
 def _build_re_case_matrix():
     penalty_multi = {
         "S": [
@@ -1131,6 +1186,7 @@ def _build_tensor_case_matrix():
     *_build_tprs_case_matrix(),
     *_build_duchon_case_matrix(),
     *_build_gp_case_matrix(),
+    *_build_sos_case_matrix(),
     *_build_mrf_case_matrix(),
     *_build_re_case_matrix(),
     *_build_factor_smooth_case_matrix(),
@@ -1370,6 +1426,25 @@ def _serialize_gp_raw(term):
     )
 
 
+def _serialize_sos_raw(term):
+    setup = term._setup
+    return _common_raw_state(
+        "sos.smooth",
+        np.asarray(setup.basis_train, dtype=np.float64),
+        [np.asarray(setup.penalty, dtype=np.float64)],
+        rank=int(setup.rank),
+        null_space_dim=int(setup.null_space_dim),
+        extra={
+            "knt": np.concatenate([setup.knots[:, 0], setup.knots[:, 1]]),
+            "UZ": np.asarray(setup.UZ, dtype=np.float64),
+            "p_order": int(setup.order),
+            "xc_scale": np.asarray(setup.column_scale, dtype=np.float64),
+            "used_supplied_knots": bool(setup.used_supplied_knots),
+            "used_subsampling": bool(setup.used_subsampling),
+        },
+    )
+
+
 def _serialize_mrf_raw(term):
     setup = term._setup
     return _common_raw_state(
@@ -1616,6 +1691,8 @@ def _serialize_term_raw(term, X):
         return _serialize_duchon_raw(term)
     if isinstance(term, GaussianProcessTerm):
         return _serialize_gp_raw(term)
+    if isinstance(term, SphericalSplineTerm):
+        return _serialize_sos_raw(term)
     if isinstance(term, MarkovRandomFieldTerm):
         return _serialize_mrf_raw(term)
     if isinstance(term, RandomEffectTerm):
diff --git a/tests/smooths/test_mgcv_smoothcon_parity.py b/tests/smooths/test_mgcv_smoothcon_parity.py
index c34b3ef0..dfbb0b0f 100644
--- a/tests/smooths/test_mgcv_smoothcon_parity.py
+++ b/tests/smooths/test_mgcv_smoothcon_parity.py
@@ -1072,6 +1072,85 @@ def test_ds_2d_custom_order_smoothcon_basis_and_penalty_match_mgcv(self):
         )
 
 
+class TestSphericalSplineSmooth:
+    """Spherical-spline smoothCon parity against mgcv 1.9-4."""
+
+    @staticmethod
+    def _make_data(seed=951, n=130):
+        rng = np.random.default_rng(seed)
+        lo = rng.uniform(-180.0, 180.0, size=n)
+        la = np.rad2deg(np.arcsin(rng.uniform(-1.0, 1.0, size=n)))
+        y = np.sin(np.deg2rad(lo)) * np.cos(np.deg2rad(la - 10.0))
+        return pd.DataFrame({"y": y, "la": la, "lo": lo})
+
+    @staticmethod
+    def _assert_basis_and_penalties(data, formula, expression, *, n_penalties=1):
+        design = _compile_formula_design(data, formula)
+        expected_x = _run_mgcv_smoothcon_matrix(data, expression)
+        actual_x = np.asarray(design.design_matrix, dtype=np.float64)
+        target_x = np.asarray(expected_x["X"], dtype=np.float64)
+        np.testing.assert_allclose(
+            actual_x @ np.linalg.pinv(actual_x),
+            target_x @ np.linalg.pinv(target_x),
+            atol=2e-8,
+            rtol=2e-8,
+        )
+        actual_s = [
+            np.asarray(block.matrix, dtype=np.float64)
+            for block in design.compiled_penalties
+        ]
+        assert len(actual_s) == n_penalties
+        if n_penalties == 0:
+            return
+        expected_s = _run_mgcv_smoothcon_penalties(
+            data, expression, absorb_cons=True, scale_penalty=True
+        )
+        penalty_payload = expected_s["S"]
+        if isinstance(penalty_payload, dict):
+            penalty_payload = list(penalty_payload.values())
+        target_s = [np.asarray(S, dtype=np.float64) for S in penalty_payload]
+        assert len(target_s) == n_penalties
+        np.testing.assert_allclose(
+            penalized_response_operator(actual_x, actual_s),
+            penalized_response_operator(target_x, target_s),
+            atol=2e-8,
+            rtol=2e-8,
+        )
+
+    def test_sos_default_smoothcon_basis_and_penalty_match_mgcv(self):
+        data = self._make_data(seed=951)
+        self._assert_basis_and_penalties(
+            data,
+            'y ~ s(la, lo, bs="sos", k=12, sp=.7)',
+            's(la, lo, bs="sos", k=12, sp=.7)',
+        )
+
+    def test_sos_duchon_tail_smoothcon_basis_and_penalty_match_mgcv(self):
+        data = self._make_data(seed=952)
+        self._assert_basis_and_penalties(
+            data,
+            'y ~ s(la, lo, bs="sos", k=12, m=-1, sp=.7)',
+            's(la, lo, bs="sos", k=12, m=-1, sp=.7)',
+        )
+
+    def test_sos_pc_smoothcon_basis_and_penalty_match_mgcv(self):
+        data = self._make_data(seed=953)
+        self._assert_basis_and_penalties(
+            data,
+            'y ~ s(la, lo, bs="sos", k=12, pc=[0,0], sp=.7)',
+            's(la, lo, bs="sos", k=12, pc=c(0,0), sp=.7)',
+        )
+
+    def test_sos_fixed_smoothcon_basis_matches_mgcv_without_penalty(self):
+        data = self._make_data(seed=954)
+        self._assert_basis_and_penalties(
+            data,
+            'y ~ s(la, lo, bs="sos", k=12, fx=True)',
+            's(la, lo, bs="sos", k=12, fx=TRUE)',
+            n_penalties=0,
+        )
+
+
 class TestPSplineSmooth(_SharedTestPSplineSmooth):
     """P-spline (bs='ps') standalone parity against mgcv."""