astropy · thiagodcoq · May 7, 2026 · May 17, 2026 · May 20, 2026
diff --git a/astropy/coordinates/tests/test_representation_methods.py b/astropy/coordinates/tests/test_representation_methods.py
@@ -365,7 +365,7 @@ def test_move_axis(self):
 
     def test_roll_axis(self):
         # Goes via transpose so works without __array_function__ as well.
-        s0_10 = np.rollaxis(self.s0, 1)
+        s0_10 = np.moveaxis(self.s0, 1, 0)
         assert s0_10.shape == (self.s0.shape[1], self.s0.shape[0])
         assert np.all(representation_equal(self.s0.T, s0_10))
         assert np.may_share_memory(s0_10.lon, self.s0.lon)

@@ -2074,7 +2074,7 @@ def _prepare_inputs_model_set(self, params, inputs, model_set_axis_input, **kwar
                     new_input = _input.reshape(new_shape)
                 else:
                     pivot = _input.ndim - len(max_param_shape) - 1
-                    new_input = np.rollaxis(_input, model_set_axis_input, pivot + 1)
+                    new_input = np.moveaxis(_input, model_set_axis_input, pivot + 1)
             pivots.append(pivot)
             reshaped.append(new_input)
 
@@ -2279,7 +2279,7 @@ def _prepare_outputs_model_set(self, outputs, broadcasted_shapes, model_set_axis
         for idx, output in enumerate(outputs):
             pivot = pivots[idx]
             if pivot < output.ndim and pivot != model_set_axis:
-                outputs[idx] = np.rollaxis(output, pivot, model_set_axis)
+                outputs[idx] = np.moveaxis(output, pivot, model_set_axis)
         return tuple(outputs)
 
     def prepare_outputs(self, broadcasted_shapes, *outputs, **kwargs):

@@ -434,7 +434,7 @@ def _add_fitting_uncertainties(self, model, a, n_coeff, x, y, z=None, resids=Non
                         mask = y.mask[..., j].ravel()
                     xx = np.ma.array(x, mask=mask)
                     eval_y = model(xx, model_set_axis=False)
-                    eval_y = np.rollaxis(eval_y, model.model_set_axis)[j]
+                    eval_y = np.moveaxis(eval_y, model.model_set_axis, 0)[j]
                     RSS.append(
                         (1 / (xx.count() - n_coeff)) * np.sum((y[..., j] - eval_y) ** 2)
                     )
@@ -461,7 +461,7 @@ def _add_fitting_uncertainties(self, model, a, n_coeff, x, y, z=None, resids=Non
                     if model.model_set_axis == 1:
                         # model_set_axis passed when evaluating only refers to input shapes
                         # so output must be reshaped for model_set_axis=1.
-                        eval_z = np.rollaxis(eval_z, 1)
+                        eval_z = np.moveaxis(eval_z, 1, 0)
                     eval_z = eval_z[j]
                     RSS.append(
                         [(1 / (len(x) - n_coeff)) * np.sum((z[j] - eval_z) ** 2)]
@@ -694,7 +694,7 @@ def __call__(
                     # dimension along which models are stacked and transpose so
                     # the model axis is *last* (I think this resolves Erik's
                     # pending generalization from 80a6f25a):
-                    rhs = np.rollaxis(z, model_axis, z.ndim)
+                    rhs = np.moveaxis(z, model_axis, -1)
                     rhs = rhs.reshape(-1, rhs.shape[-1])
                 else:
                     # This "else" seems to handle the corner case where the
@@ -707,7 +707,7 @@ def __call__(
                     # Same for weights
                     if weights.ndim > 2:
                         # Separate 2D weights for each model:
-                        weights = np.rollaxis(weights, model_axis, weights.ndim)
+                        weights = np.moveaxis(weights, model_axis, -1)
                         weights = weights.reshape(-1, weights.shape[-1])
                     elif weights.ndim == z.ndim:
                         # Separate, flattened weights for each model:
@@ -799,7 +799,7 @@ def __call__(
             # Arrange the lhs as a stack of 2D matrices that we can iterate
             # over to get the correctly-orientated lhs for each model:
             if lhs.ndim > 2:
-                lhs_stack = np.rollaxis(lhs, -1, 0)
+                lhs_stack = np.moveaxis(lhs, -1, 0)
             else:
                 lhs_stack = np.broadcast_to(lhs, rhs.shape[-1:] + lhs.shape)
 
@@ -1049,11 +1049,11 @@ def __call__(self, model, x, y, z=None, weights=None, *, inplace=False, **kwargs
                         # Get views transposed appropriately for iteration
                         # over the set (handling data & mask separately due to
                         # NumPy issue #8506):
-                        data_T = np.rollaxis(filtered_data, model_set_axis, 0)
-                        mask_T = np.rollaxis(filtered_data.mask, model_set_axis, 0)
+                        data_T = np.moveaxis(filtered_data, model_set_axis, 0)
+                        mask_T = np.moveaxis(filtered_data.mask, model_set_axis, 0)
 
             if loop:
-                model_vals_T = np.rollaxis(model_vals, model_set_axis, 0)
+                model_vals_T = np.moveaxis(model_vals, model_set_axis, 0)
                 for row_data, row_mask, row_mod_vals in zip(
                     data_T, mask_T, model_vals_T
                 ):
@@ -2119,7 +2119,7 @@ def _convert_input(x, y, z=None, n_models=1, model_set_axis=0):
             # That is, each model should be represented by a column.  TODO:
             # Obviously this is a detail of np.linalg.lstsq and should be
             # handled specifically by any fitters that use it...
-            y = np.rollaxis(y, model_set_axis, y.ndim)
+            y = np.moveaxis(y, model_set_axis, -1)
             data_shape = y.shape[:-1]
         else:
             # Shape of z excluding model_set_axis

@@ -541,7 +541,7 @@ def fit_deriv(self, x, *params):
             v[1] = x
             for i in range(2, self.degree + 1):
                 v[i] = v[i - 1] * x2 - v[i - 2]
-        return np.rollaxis(v, 0, v.ndim)
+        return np.moveaxis(v, 0, -1)
 
     def prepare_inputs(self, x, **kwargs):
         inputs, broadcasted_shapes = super().prepare_inputs(x, **kwargs)
@@ -662,7 +662,7 @@ def fit_deriv(self, x, *params):
             v[1] = 2 * x
             for i in range(2, self.degree + 1):
                 v[i] = x2 * v[i - 1] - 2 * (i - 1) * v[i - 2]
-        return np.rollaxis(v, 0, v.ndim)
+        return np.moveaxis(v, 0, -1)
 
     def prepare_inputs(self, x, **kwargs):
         inputs, broadcasted_shapes = super().prepare_inputs(x, **kwargs)
@@ -841,10 +841,7 @@ def _hermderiv1d(self, x, deg):
             d[1] = x2
             for i in range(2, deg + 1):
                 d[i] = x2 * d[i - 1] - 2 * (i - 1) * d[i - 2]
-        return np.rollaxis(d, 0, d.ndim)
-
-
-class Legendre1D(_PolyDomainWindow1D):
+        return np.moveaxis(d, 0, -1)
     r"""
     Univariate Legendre series.
 
@@ -943,7 +940,7 @@ def fit_deriv(self, x, *params):
             v[1] = x
             for i in range(2, self.degree + 1):
                 v[i] = (v[i - 1] * x * (2 * i - 1) - v[i - 2] * (i - 1)) / i
-        return np.rollaxis(v, 0, v.ndim)
+        return np.moveaxis(v, 0, -1)
 
     @staticmethod
     def clenshaw(x, coeffs):
@@ -1061,7 +1058,7 @@ def fit_deriv(self, x, *params):
             v[1] = x
             for i in range(2, self.degree + 1):
                 v[i] = v[i - 1] * x
-        return np.rollaxis(v, 0, v.ndim)
+        return np.moveaxis(v, 0, -1)
 
     @staticmethod
     def horner(x, coeffs):
@@ -1497,7 +1494,7 @@ def _chebderiv1d(self, x, deg):
             d[1] = x
             for i in range(2, deg + 1):
                 d[i] = d[i - 1] * x2 - d[i - 2]
-        return np.rollaxis(d, 0, d.ndim)
+        return np.moveaxis(d, 0, -1)
 
 
 class Legendre2D(OrthoPolynomialBase):
@@ -1650,7 +1647,7 @@ def _legendderiv1d(self, x, deg):
             d[1] = x
             for i in range(2, deg + 1):
                 d[i] = (d[i - 1] * x * (2 * i - 1) - d[i - 2] * (i - 1)) / i
-        return np.rollaxis(d, 0, d.ndim)
+        return np.moveaxis(d, 0, -1)
 
 
 class _SIP1D(PolynomialBase):

@@ -885,7 +885,7 @@ def test_2d_set_axis_2_fitting_with_outlier_removal():
     )
 
     y, x = np.mgrid[0:5, 0:5]
-    z = np.rollaxis(np.array([x + y, 1 - 0.1 * x + 0.2 * y]), 0, 3)
+    z = np.moveaxis(np.array([x + y, 1 - 0.1 * x + 0.2 * y]), 0, -1)
     z[3, 3:5, 0] = 100.0  # outliers
 
     poly_set, filt_z = fitter(poly_set, x, y, z)

@@ -284,7 +284,7 @@ def test_linearlsqfitter(model_class):
 
     # Generate data for fitting 2 models and re-stack them along the last axis:
     y = np.array([2 * x + 1, x + 4])
-    y = np.rollaxis(y, 0, -1).T
+    y = np.moveaxis(y, 0, -1).T
 
     f = LinearLSQFitter()
     # This seems to fit the model_set correctly:
@@ -305,7 +305,7 @@ def test_model_set_axis_outputs():
     model_set = Polynomial2D(1, n_models=2, model_set_axis=2)
     y2, x2 = np.mgrid[:5, :5]
     # z = np.moveaxis([x2 + y2, 1 - 0.1 * x2 + 0.2 * y2]), 0, 2)
-    z = np.rollaxis(np.array([x2 + y2, 1 - 0.1 * x2 + 0.2 * y2]), 0, 3)
+    z = np.moveaxis(np.array([x2 + y2, 1 - 0.1 * x2 + 0.2 * y2]), 0, 3)
     model = fitter(model_set, x2, y2, z)
     res = model(x2, y2, model_set_axis=False)
     assert z.shape == res.shape
@@ -594,7 +594,7 @@ def test_linear_2d_separate_weights_axis_2(self):
             model,
             self.x2,
             self.y2,
-            np.rollaxis(self.z2, 0, 3),
+            np.moveaxis(self.z2, 0, 3),
             weights=self.w2[..., np.newaxis],
         )
         assert_allclose(model.c0_0, 1.0, atol=1e-12)

@@ -1153,9 +1153,9 @@ def test_two_model_nonzero_model_set_axis(self):
         # An example where the model set axis is the *last* axis of the
         # parameter arrays
         coeff = np.array([[[10, 20, 30], [30, 40, 50]], [[50, 60, 70], [70, 80, 90]]])
-        coeff = np.rollaxis(coeff, 0, 3)
+        coeff = np.moveaxis(coeff, 0, 3)
         e = np.array([[1, 2, 3], [3, 4, 5]])
-        e = np.rollaxis(e, 0, 2)
+        e = np.moveaxis(e, 0, 2)
         t = TParModel(coeff, e, n_models=2, model_set_axis=-1)
         assert len(t) == 2
         assert t.model_set_axis == -1

@@ -409,8 +409,8 @@ def test_sigma_clip_axis_tuple_3D():
     data = np.sin(0.78 * np.arange(27)).reshape(3, 3, 3)
     mask = np.zeros_like(data, dtype=np.bool_)
 
-    data_t = np.rollaxis(data, 1, 0)
-    mask_t = np.rollaxis(mask, 1, 0)
+    data_t = np.moveaxis(data, 1, 0)
+    mask_t = np.moveaxis(mask, 1, 0)
 
     # Loop over what was originally axis 1 and clip each plane directly:
     for data_plane, mask_plane in zip(data_t, mask_t):

@@ -434,7 +434,7 @@ def test_roll_axis(self, use_mask):
         # Goes via transpose so works without __array_function__ as well.
         self.create_data(use_mask)
 
-        t0_10 = np.rollaxis(self.t0, 1)
+        t0_10 = np.moveaxis(self.t0, 1, 0)
         assert t0_10.shape == (self.t0.shape[1], self.t0.shape[0])
         assert_time_all_equal(self.t0.T, t0_10)
         assert np.may_share_memory(t0_10.jd1, self.t0.jd1)

diff --git a/astropy/units/tests/test_quantity_non_ufuncs.py b/astropy/units/tests/test_quantity_non_ufuncs.py
@@ -120,9 +120,6 @@ def test_ravel(self):
     def test_moveaxis(self):
         self.check(np.moveaxis, 0, 1)
 
-    def test_rollaxis(self):
-        self.check(np.rollaxis, 0, 2)
-
     def test_swapaxes(self):
         self.check(np.swapaxes, 0, 1)
 

diff --git a/astropy/utils/masked/tests/test_function_helpers.py b/astropy/utils/masked/tests/test_function_helpers.py
@@ -96,8 +96,8 @@ def test_ravel(self):
     def test_moveaxis(self):
         self.check(np.moveaxis, 0, 1)
 
-    def test_rollaxis(self):
-        self.check(np.rollaxis, 0, 2)
+    def test_moveaxis_rollaxis_compat(self):
+        self.check(np.moveaxis, 0, 2)
 
     def test_swapaxes(self):
         self.check(np.swapaxes, 0, 1)