diff --git a/yt/_maintenance/numpy2_compat.py b/yt/_maintenance/numpy2_compat.py
index a5f13bc8121..5f53beb91bd 100644
--- a/yt/_maintenance/numpy2_compat.py
+++ b/yt/_maintenance/numpy2_compat.py
@@ -1,101 +1,12 @@
-# vendor functions that were moved from numpy 1.x to scipy
+# avoid deprecation warnings in numpy >= 2.0
 
-import functools
+from importlib.metadata import version
 
-import numpy as np
-from numpy.core import overrides
+from packaging.version import Version
 
-array_function_dispatch = functools.partial(
-    overrides.array_function_dispatch, module="yt"
-)
+NUMPY_VERSION = Version(version("numpy"))
 
-
-def _trapezoid_dispatcher(y, x=None, dx=None, axis=None):
-    return (y, x)
-
-
-# from numpy 1.25 (numpy.trapz), deprecated in numpy 2.0 the function is vendored to
-# avoid adding a runtime dependency on scipy.integrate.trapezoid the name of the
-# original numpy function is also avoided,
-# see https://github.com/scipy/scipy/issues/12924
-@array_function_dispatch(_trapezoid_dispatcher)
-def trapezoid(y, x=None, dx=1.0, axis=-1):
-    r"""
-    Integrate along the given axis using the composite trapezoidal rule.
-
-    If `x` is provided, the integration happens in sequence along its
-    elements - they are not sorted.
-
-    Integrate `y` (`x`) along each 1d slice on the given axis, compute
-    :math:`\int y(x) dx`.
-    When `x` is specified, this integrates along the parametric curve,
-    computing :math:`\int_t y(t) dt =
-    \int_t y(t) \left.\frac{dx}{dt}\right|_{x=x(t)} dt`.
-
-    Parameters
-    ----------
-    y : array_like
-        Input array to integrate.
-    x : array_like, optional
-        The sample points corresponding to the `y` values. If `x` is None,
-        the sample points are assumed to be evenly spaced `dx` apart. The
-        default is None.
-    dx : scalar, optional
-        The spacing between sample points when `x` is None. The default is 1.
-    axis : int, optional
-        The axis along which to integrate.
-
-    Returns
-    -------
-    trapezoid : float or ndarray
-        Definite integral of `y` = n-dimensional array as approximated along
-        a single axis by the trapezoidal rule. If `y` is a 1-dimensional array,
-        then the result is a float. If `n` is greater than 1, then the result
-        is an `n`-1 dimensional array.
-
-    See Also
-    --------
-    sum, cumsum
-
-    Notes
-    -----
-    Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
-    will be taken from `y` array, by default x-axis distances between
-    points will be 1.0, alternatively they can be provided with `x` array
-    or with `dx` scalar.  Return value will be equal to combined area under
-    the red lines.
-
-
-    References
-    ----------
-    .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
-
-    .. [2] Illustration image:
-           https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
-    """
-    y = np.asanyarray(y)
-    if x is None:
-        d = dx
-    else:
-        x = np.asanyarray(x)
-        if x.ndim == 1:
-            d = np.diff(x)
-            # reshape to correct shape
-            shape = [1] * y.ndim
-            shape[axis] = d.shape[0]
-            d = d.reshape(shape)
-        else:
-            d = np.diff(x, axis=axis)
-    nd = y.ndim
-    slice1 = [slice(None)] * nd
-    slice2 = [slice(None)] * nd
-    slice1[axis] = slice(1, None)
-    slice2[axis] = slice(None, -1)
-    try:
-        ret = (d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0).sum(axis)
-    except ValueError:
-        # Operations didn't work, cast to ndarray
-        d = np.asarray(d)
-        y = np.asarray(y)
-        ret = np.add.reduce(d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0, axis)
-    return ret
+if NUMPY_VERSION >= Version("2.0.0dev0"):
+    from numpy import trapezoid as trapezoid  # type: ignore [attr-defined]
+else:
+    from numpy import trapz as trapezoid  # noqa: F401