Warn on missing sphere_radius, set gradient/curl units, add vorticity example

- _compute_gradient: emit UserWarning when scale_by_radius=True but the grid
  lacks a sphere_radius attribute, instead of silently leaving the result on
  the unit sphere.
- Set units attrs on gradient and curl outputs so callers can distinguish
  physical (/m, 1/s) from unit-sphere (/rad) results.
- Add unit assertions to the existing radius-scaling tests.
- vector_calculus user guide: units callout under Gradient and a new
  solid-body-rotation vorticity example that validates curl against the
  analytic zeta = -2 Omega sin(phi).

Author: rajeeja

docs/user-guide/vector_calculus.ipynb

@@ -120,6 +120,12 @@
     "Gradients can be computed using the `UxDataArray.gradient()` method on a face-centered data variable."
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "66caacb5",
+   "metadata": {},
+   "source": "> **Units.** By default, `gradient()` and `curl()` divide by `uxgrid.sphere_radius` so derivatives carry physical units (e.g. `[data units]/m`, `1/s` for velocity curl). Pass `scale_by_radius=False` to keep results on the unit sphere (per radian). If the grid has no `sphere_radius` attribute, the call falls back to unit-sphere output and emits a `UserWarning`."
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -440,6 +446,70 @@
     "(vortex_vectors + curl_magnitude).opts(shared_axes=False)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "6e67a135",
+   "metadata": {},
+   "source": "### Example 3: Relative Vorticity from Solid-Body Rotation\n\nThe 2D curl on the sphere is **relative vorticity** $\\zeta = \\partial v/\\partial x - \\partial u/\\partial y$, a quantity meteorologists and oceanographers care about every day. With the default `scale_by_radius=True`, `u.curl(v)` returns $\\zeta$ in physical units of $s^{-1}$.\n\nA clean analytical check is **solid-body rotation about the polar axis** with angular speed $\\Omega$:\n\n$$\nu(\\varphi) = \\Omega R \\cos\\varphi, \\qquad v = 0\n$$\n\nFor this flow the relative vorticity on a sphere of radius $R$ is\n\n$$\n\\zeta = -\\frac{1}{R\\cos\\varphi}\\frac{\\partial(u\\cos\\varphi)}{\\partial \\varphi} = -2\\Omega \\sin\\varphi.\n$$\n\nWe construct the field on the existing MPAS subset and compare `u.curl(v)` against this closed form."
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3b5f3090",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "OMEGA = 7.292115e-5  # Earth's angular speed (rad/s)\n",
+    "R = uxds.uxgrid.sphere_radius\n",
+    "lat_rad = np.deg2rad(uxds.uxgrid.face_lat.values)\n",
+    "\n",
+    "u_sbr_data = OMEGA * R * np.cos(lat_rad)\n",
+    "v_sbr_data = np.zeros_like(u_sbr_data)\n",
+    "\n",
+    "u_sbr = ux.UxDataArray(\n",
+    "    u_sbr_data,\n",
+    "    dims=[\"n_face\"],\n",
+    "    uxgrid=uxds.uxgrid,\n",
+    "    name=\"u_sbr\",\n",
+    "    attrs={\"units\": \"m/s\"},\n",
+    ")\n",
+    "v_sbr = ux.UxDataArray(\n",
+    "    v_sbr_data,\n",
+    "    dims=[\"n_face\"],\n",
+    "    uxgrid=uxds.uxgrid,\n",
+    "    name=\"v_sbr\",\n",
+    "    attrs={\"units\": \"m/s\"},\n",
+    ")\n",
+    "\n",
+    "zeta = u_sbr.curl(v_sbr)\n",
+    "zeta_analytic = -2.0 * OMEGA * np.sin(lat_rad)\n",
+    "\n",
+    "finite = np.isfinite(zeta.values)\n",
+    "err = np.abs(zeta.values[finite] - zeta_analytic[finite])\n",
+    "\n",
+    "print(f\"curl units: {zeta.attrs['units']}\")\n",
+    "print(f\"sphere_radius (m): {R:.3e}\")\n",
+    "print(\n",
+    "    f\"computed   zeta range: [{zeta.values[finite].min():.3e}, {zeta.values[finite].max():.3e}]  1/s\"\n",
+    ")\n",
+    "print(\n",
+    "    f\"analytic   zeta range: [{zeta_analytic.min():.3e}, {zeta_analytic.max():.3e}]  1/s\"\n",
+    ")\n",
+    "print(f\"max |error|: {err.max():.3e}  1/s  (≈ {err.max() / (2 * OMEGA):.1%} of 2Omega)\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "434634a9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "zeta.plot(cmap=\"RdBu_r\", aspect=1).opts(\n",
+    "    title=\"Relative Vorticity (1/s) — Solid-Body Rotation\", colorbar=True\n",
+    ")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "23",

test/core/test_vector_calculus.py

@@ -162,6 +162,9 @@ def test_gradient_scales_by_radius(self, gridpath, datasetpath):
             equal_nan=True,
         )
 
+        assert grad_scaled["zonal_gradient"].attrs["units"].endswith("/m")
+        assert grad_unit["zonal_gradient"].attrs["units"].endswith("/rad")
+
 
 class TestDivergenceQuadHex:
 
@@ -559,6 +562,9 @@ def test_curl_scales_by_radius(self, gridpath, datasetpath):
 
         nt.assert_allclose(curl_scaled, curl_unit / radius, equal_nan=True)
 
+        assert curl_scaled.attrs["units"]
+        assert curl_unit.attrs["units"].endswith("/rad")
+
     def test_curl_units_and_attributes(self, gridpath, datasetpath):
         """Test that curl preserves appropriate units and attributes"""
         uxds = ux.open_dataset(

uxarray/core/dataarray.py

@@ -1535,15 +1535,20 @@ def curl(
         # Compute curl = ∂v/∂x - ∂u/∂y
         curl_values = grad_v_zonal.values - grad_u_meridional.values
 
+        u_units = self.attrs.get("units", "")
+        has_sphere_radius = "sphere_radius" in self.uxgrid._ds.attrs
+        if scale_by_radius and has_sphere_radius:
+            curl_units = f"({u_units})/m" if u_units else "1/s"
+        else:
+            curl_units = f"({u_units})/rad" if u_units else "1/rad"
+
         # Create the result UxDataArray
         curl_da = UxDataArray(
             curl_values,
             dims=self.dims,
             attrs={
                 "long_name": f"Curl of ({self.name}, {other.name})",
-                "units": "1/s"
-                if "units" not in self.attrs
-                else f"({self.attrs.get('units', '1')})/m",
+                "units": curl_units,
                 "description": "Curl of vector field computed as ∂v/∂x - ∂u/∂y",
             },
             uxgrid=self.uxgrid,

uxarray/core/gradient.py

@@ -1,3 +1,5 @@
+import warnings
+
 import numpy as np
 from numba import njit, prange
 
@@ -186,20 +188,38 @@ def _compute_gradient(data, scale_by_radius=True):
             "Computing the gradient is only supported for face-centered data variables."
         )
 
+    has_sphere_radius = "sphere_radius" in uxgrid._ds.attrs
     if scale_by_radius:
-        radius = uxgrid.sphere_radius
-        grad_zonal = grad_zonal / radius
-        grad_meridional = grad_meridional / radius
+        if has_sphere_radius:
+            radius = uxgrid.sphere_radius
+            grad_zonal = grad_zonal / radius
+            grad_meridional = grad_meridional / radius
+        else:
+            warnings.warn(
+                "scale_by_radius=True but the grid has no 'sphere_radius' "
+                "attribute; result is left on the unit sphere. Set "
+                "uxgrid.sphere_radius or pass scale_by_radius=False.",
+                UserWarning,
+                stacklevel=2,
+            )
+
+    base_units = data.attrs.get("units", "")
+    if scale_by_radius and has_sphere_radius:
+        grad_units = f"{base_units}/m" if base_units else "1/m"
+    else:
+        grad_units = f"{base_units}/rad" if base_units else "1/rad"
 
     # Zonal
     grad_zonal_da = UxDataArray(
         data=grad_zonal, name="zonal_gradient", dims=data.dims, uxgrid=uxgrid
     )
+    grad_zonal_da.attrs["units"] = grad_units
 
     # Meridional
     grad_meridional_da = UxDataArray(
         data=grad_meridional, name="meridional_gradient", dims=data.dims, uxgrid=uxgrid
     )
+    grad_meridional_da.attrs["units"] = grad_units
 
     return grad_zonal_da, grad_meridional_da