The millihertz gravitational-wave frequency band is expected to contain a rich symphony of signals with sources ranging from galactic white dwarf binaries to extreme mass ratio inspirals. Many of these gravitational-wave signals will not be individually resolvable. Instead, they will incoherently add to produce stochastic gravitational-wave confusion noise whose frequency content will be governed by the dynamics of the sources. The angular structure of the power of the confusion noise will be modulated by the distribution of the sources across the sky. Measurement of this structure can yield important information about the distribution of sources on galactic and extra-galactic scales, their astrophysics and their evolution over cosmic timescales. Moreover, since the confusion noise is part of the noise budget of LISA, mapping it will also be essential for studying resolvable signals. In this paper, we present a Bayesian algorithm to probe the angular distribution of the stochastic gravitational-wave confusion noise with LISA using a spherical harmonic basis. We develop a technique based on Clebsch-Gordan coefficients to mathematically constrain the spherical harmonics to yield a non-negative distribution, making them optimal for expanding the gravitational-wave power and amenable to Bayesian inference. We demonstrate these techniques using a series of simulations and analyses, including recovery of simulated distributed and localized sources of gravitational-wave power. We also apply this method to map the gravitational-wave foreground from galactic white-dwarfs using a simplified model of the galactic white dwarf distribution.