Repository: astropy/astropy

Issue:
scipy and quantities
The example below works with released versions of scipy but fails with the development version. Correct me if I'm wrong but I think the example worked by accident. And this will need to be fixed in modeling.

```python
from scipy.interpolate import interpn
from astropy import units as u

points = (np.array([0, 1, 2, 4]) * u.pix,)
lookup_table = np.array([1, 2, 3, 4]) * u.one
inputs = np.array([3])

q = interpn(points, lookup_table, inputs*u.pix, method="nearest", bounds_error=False, fill_value=np.nan)
```
With released versions of scipy this returns a Quantity object, with the dev version it returns an array.
This only happens with the `nearest` interpolation method. It causes an [error](https://github.com/spacetelescope/gwcs/actions/runs/19424008864/job/55567202388) in gwcs.

Discussion:

Fix (patch):
diff --git a/astropy/modeling/tabular.py b/astropy/modeling/tabular.py
index 7203ee142cdd..15dce2aa72dc 100644
--- a/astropy/modeling/tabular.py
+++ b/astropy/modeling/tabular.py
@@ -247,7 +247,7 @@ def evaluate(self, *inputs):
 
         # return_units not respected when points has no units
         if isinstance(self.lookup_table, u.Quantity) and not isinstance(
-            self.points[0], u.Quantity
+            result, u.Quantity
         ):
             result = result * self.lookup_table.unit
 
diff --git a/docs/changes/modeling/18958.bugfix.rst b/docs/changes/modeling/18958.bugfix.rst
new file mode 100644
index 000000000000..babcd756541b
--- /dev/null
+++ b/docs/changes/modeling/18958.bugfix.rst
@@ -0,0 +1,1 @@
+Fixed a bug in ``modeling.tabular`` models when the ``lookup_table`` is a Quantity, where the result might lose its unit in some cases.


Test changes:
diff --git a/astropy/modeling/tests/test_models.py b/astropy/modeling/tests/test_models.py
index fc401f1b3e5b..40567f829c13 100644
--- a/astropy/modeling/tests/test_models.py
+++ b/astropy/modeling/tests/test_models.py
@@ -729,6 +729,17 @@ def test_tabular_interp_1d():
     assert_quantity_allclose(model(xnew.to(u.nm)), ans1)
     assert model.bounding_box == (0 * u.nm, 4 * u.nm)
 
+    # Test with no units on points.
+    model = LookupTable(points=points, lookup_table=values * u.nJy)
+    assert_quantity_allclose(model(xnew), ans1)
+    assert_quantity_allclose(model(xnew.to(u.nm)), ans1)
+    assert model.bounding_box == (0, 4)
+
+    model = LookupTable(points=points, lookup_table=values * u.nJy, method="nearest")
+    assert_quantity_allclose(model(xnew), np.array([1, 10, 10, 2, -3]) * u.nJy)
+    assert_quantity_allclose(model(xnew.to(u.nm)), np.array([1, 10, 10, 2, -3]) * u.nJy)
+    assert model.bounding_box == (0, 4)
+
     # Test fill value unit conversion and unitless input on table with unit
     model = LookupTable(
         [1, 2, 3],