def ndintegrate(data, unit_conv, limits, unit='ppm', noise_limits=None):
    """
    Integrate one nD data array within limits given in units. Data points must
    be equally spaced. Can only integrate one region per function call.

    The integration error due to baseline noise is calculated as:

    .. math::
        error = \sigma_{vol} = \sigma \sqrt{\prod_i^{d} n_{i}},

    if the noise_limits are set.

    Where:
    sigma is the standard deviation of the baseline noise. n is the number
    of bins in the integration range for each d dimensions.

    See integrate for more information.

    Parameters
    ----------
    data: array like
        1d array of intensity
    unit_convs: [`fileiobase.unit_conversion`, ] list
        list of unit_conversion object associated with each dim of data.
    limits: array like
        With shape (2,) or (d, 2). 1D Array with lower and upper integration
        limits for 1D . Or array with d rows of lower and upper integration
        limits for each dimension.
    noise_limits: Optional[array like]
        With shape(2, ). Array with lower and upper limits to section of data
        with only noise. A larger range will likely yield a more accurate
        estimate. It is unwise to use the very end of the spectrum for the
        noise range.

    Returns
    -------
    array
        [value, ...] integration values

    if noise_limits is given:

    array
        [[value, error], ...] where error a one sigma estimate of the error
        only from the spectrum noise
    """

    # determine the dimensionality of the data.
    d = np.ndim(data)

    try:
        iter(unit_conv)
    except __HOLE__:
        unit_conv = [unit_conv, ]

    if d != len(unit_conv):
        mesg = 'A unit_conversion object is needed for each dimension.'
        raise ValueError(mesg)

    limits = np.array(limits)
    if limits.ndim == 1:
        limits = np.expand_dims(limits, axis=0)

    if limits.shape[0] != d and limits.shape[1] != 2:
        mesg = 'A lower and upper limit is needed for each dimension.'
        raise ValueError(mesg)

    inds = [(uc(x[0], unit), uc(x[1], unit))
            for (uc, x) in zip(unit_conv, limits)]
    inds = [sorted(x) for x in inds]

    # the integrate_nd needs to be scaled by the bin width in ppm
    ppm_scales = [x.ppm_scale() for x in unit_conv]
    dx = np.prod(np.array([abs(x[1]-x[0]) for x in ppm_scales]))

    slice_sum = (data[[slice(x[0], x[1])for x in np.flipud(inds)]]).sum()

    value = slice_sum * dx

    if noise_limits is None:
        return value

    else:
        noise_limits = np.array(noise_limits)
        if noise_limits.size == 2:
            noise_limits = np.expand_dims(noise_limits, axis=0)

        if noise_limits.shape[0] != d and noise_limits.shape[1] != 2:
            mesg = 'If given, a noise limit is needed for each dimension.'
            raise ValueError(mesg)

        noise_inds = [(uc(x[0], unit), uc(x[1], unit))
                      for (uc, x) in zip(unit_conv, noise_limits)]
        noise_inds = [sorted(x) for x in noise_inds]

        # see docstring of integrate
        nm = np.prod(np.array([abs(x[1]-x[0]) for x in noise_inds]))
        std = np.std(data[[slice(x[0], x[1])for x in np.flipud(noise_inds)]])

        error = std * np.sqrt(nm)
        return np.hstack((value, error))
dataset/ETHPy150Open jjhelmus/nmrglue/nmrglue/analysis/integration.py/ndintegrate