For an arbitrary Markov operator P on a Lebesgue measure space ( X, m) we construct a projection S (called harmonic) from L ∞( X, m) onto the space of bounded P-harmonic functions H ∞( X, m, P). The projection S = S λ is obtained by applying a fixed measure-linear (medial) mean λ on ℤ + to the sequence of one-dimensional distributions of the Markov measure in the path space of the Markov chain { x n } associated with the operator P. If there are no non-constant bounded P-harmonic functions, S is a projection onto the space of constants. In the general situation when the space H ∞( X, m, P) is not necessarily trivial, the harmonic projection Sf can still be considered as a “space average” of f. We show that for any f ∈ L ∞( X, m) the harmonic projection Sf can be recovered from the averages of f (determined by the mean λ) along sample paths of the Markov chain { x n } by an integral Poisson formula in the same way as any bounded harmonic function is represented by the Poisson integral of its boundary values. In other words, for any f ∈ L ∞( X, m) its space average Sf is the Poisson integral of the time averages along sample paths.