def integrate(da, over_dim, x=None, dx=None, method='trapz', cumulative=False, skipna=False):
    """ 
        Returns trapezoidal/rectangular integration along specified dimension 
        
        | Author: Dougie Squire
        | Date: 16/08/2018
        
        Parameters
        ----------
        da : xarray DataArray
            Array containing values to integrate
        over_dim : str
            Dimension to integrate
        x : xarray DataArray, optional
            Values to use for integrand. Must contain dimensions over_dim. If None, x is determined\
                    from the coords associated with over_dim
        dx : value, optional
            Integrand spacing used to compute the integral. If None, dx is determined from x
        method : str, optional
            Method of performing integral. Options are 'trapz' for trapezoidal integration, or 'rect'\
                    for rectangular integration
        cumulative : bool, optional
            If True, return the cumulative integral    
            
        Returns
        -------
        integral : xarray DataArray
            Array containing the integral along the specified dimension
            
        Examples
        --------
        >>> A = xr.DataArray(np.random.normal(size=(3,2)), coords=[('x', np.arange(3)), 
        ...                                                        ('y', np.arange(2))])
        >>> doppyo.utils.integrate(A, over_dim='x')
        <xarray.DataArray 'integral' (y: 2)>
        array([-0.20331 , -0.781251])
        Coordinates:
          * y        (y) int64 0 1

        See Also
        --------
        numpy.trapz
    """

    if x is None:
        x = da[over_dim]
    
    if len(x) == 1:
        if dx is None:
            raise ValueError('Must provide dx for integration along dimension with length 1')
        integral = (da * dx).drop(over_dim).squeeze()
    elif method == 'trapz':
        if dx is None:
            dx = x - x.shift(**{over_dim:1})
            dx = dx.fillna(0.0)

        if cumulative:
            integral = ((da.shift(**{over_dim:1}) + da) * dx / 2.0) \
                       .fillna(0.0) \
                       .cumsum(over_dim, skipna=skipna) 
        else:
            integral = ((da.shift(**{over_dim:1}) + da) * dx / 2.0) \
                       .fillna(0.0) \
                       .sum(over_dim, skipna=skipna) 
    elif method == 'rect':
        if dx is None:
            dx1 = x - x.shift(**{over_dim:1})
            dx2 = -(x - x.shift(**{over_dim:-1}))
            dx = dx1.combine_first(dx2)

        if cumulative:
            integral = (da * dx).cumsum(over_dim, skipna=skipna) 
        else:
            integral = (da * dx).sum(over_dim, skipna=skipna) 
    else:
        raise ValueError(f'{method} is not a recognised integration method')
    
    return integral