def convolve_map_hpx_gauss(m, sigmas, imin=0, imax=None, wmap=None):
    """
    Perform an energy-dependent convolution on a sequence of 2-D spatial maps.

    Parameters
    ----------

    m : `HpxMap`
       2-D map containing a sequence of 1-D HEALPix maps.  First
       dimension should be energy.

    sigmas : `~numpy.ndarray`
       1-D map containing a sequence gaussian widths for smoothing

    imin : int
       Minimum index in energy dimension.

    imax : int
       Maximum index in energy dimension.

    wmap :  `~numpy.ndarray`
       2-D map containing a sequence of 1-D HEALPix maps of weights.  First
       dimension should be energy. This map should have the same dimension as m.

    """
    islice = slice(imin, imax)

    o = np.zeros(m.data.shape)

    nside = m.geom.nside
    nest = m.geom.nest

    # Loop over energy
    for i, ms in enumerate(m.data[islice, ...]):
        sigma = sigmas[islice][i]
        # Need to be in RING scheme
        if nest:
            ms = hp.pixelfunc.reorder(ms, n2r=True)

        o[islice, ...][i] = hp.sphtfunc.smoothing(ms, sigma=sigma)
        if nest:
            o[islice, ...][i] = hp.pixelfunc.reorder(
                o[islice, ...][i], r2n=True)
        if wmap is not None:
            o[islice, ...][i] *= wmap.data[islice, ...][i]

    return HpxNDMap(m.geom, o)