def _get_surface_index(layers, n_layers, surface_index):
    """Get index within each stack of the layer at the topographic surface.

    Parameters
    ----------
    layers : ndarray of shape `(n_layers, n_nodes)`
        Array of layer thicknesses. This array is the datastructure that
        contains all allocated layers, active or inactive.
    n_layers : int
        Number of active layers.
    surface_index : ndarray of shape `(n_nodes, )`

    Examples
    --------
    >>> import numpy as np
    >>> from landlab.layers.eventlayers import (_deposit_or_erode,
    ...                                         _get_surface_index)

    >>> layers = np.full((5, 3), 1.)
    >>> dz = np.array([-1., -2., -3.])

    Note here, if you are very confused by the use of ``_deposit_or_erode``
    we recommend you read the docstring associated with that function.

    >>> _deposit_or_erode(layers, 5, dz)
    >>> layers
    array([[ 1.,  1.,  1.],
           [ 1.,  1.,  1.],
           [ 1.,  1.,  0.],
           [ 1.,  0.,  0.],
           [ 0.,  0.,  0.]])

    >>> surface_index = np.empty(3, dtype=int)
    >>> _get_surface_index(layers, 5, surface_index)
    >>> surface_index
    array([3, 2, 1])
    """
    from .ext.eventlayers import get_surface_index
    layers = layers.reshape((layers.shape[0], -1))
    get_surface_index(layers, n_layers, surface_index)
