BUG: Only allow 1d distance arrays

2d arrays would work, but in unpredictable and undocumented ways.

This at least makes gh-9401 give a better error message.

2 files changed, 7 insertions, 1 deletions

diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py
index 989030e03..91e7dc616 100644
--- a/numpy/lib/function_base.py
+++ b/numpy/lib/function_base.py
@@ -1697,8 +1697,10 @@ def gradient(f, *varargs, **kwargs):
         for i, distances in enumerate(dx):
             if distances.ndim == 0:
                 continue
+            elif distances.ndim != 1:
+                raise ValueError("distances must be either scalars or 1d")
             if len(distances) != f.shape[axes[i]]:
-                raise ValueError("distances must be either scalars or match "
+                raise ValueError("when 1d, distances must match "
                                  "the length of the corresponding dimension")
             diffx = np.diff(dx[i])
             # if distances are constant reduce to the scalar case
diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py
index ee2cc166a..d7d00758e 100644
--- a/numpy/lib/tests/test_function_base.py
+++ b/numpy/lib/tests/test_function_base.py
@@ -756,6 +756,10 @@ class TestGradient(TestCase):
         # 2D but axis specified
         gradient(f_2d, dx, axis=1)
 
+        # 2d coordinate arguments are not yet allowed
+        assert_raises_regex(ValueError, '.*scalars or 1d',
+            gradient, f_2d, np.stack([dx]*2, axis=-1), 1)
+
     def test_badargs(self):
         f_2d = np.arange(25).reshape(5, 5)
         x = np.cumsum(np.ones(5))