In the near future, cosmology will enter the wide and deep galaxy survey area allowing high-precision studies of the large scale structure of the universe in three dimensions. To test cosmological models and determine their parameters accurately, it is natural to confront data with exact theoretical expectations expressed in the observational parameter space (angles and redshift). The data-driven galaxy number count fluctuations on redshift shells, can be used to build correlation functions $C(theta; z_1, z_2)$ on and between shells which can probe the baryonic acoustic oscillations, the distance-redshift distortions as well as gravitational lensing and other relativistic effects. Transforming the model to the data space usually requires the computation of the angular power spectrum $C_ell(z_1, z_2)$ but this appears as an artificial and inefficient step plagued by apodization issues. In this article we show that it is not necessary and present a compact expression for $C(theta; z_1, z_2)$ that includes directly the leading density and redshift space distortions terms from the full linear theory. It can be evaluated using a fast integration method based on Clenshaw-Curtis quadrature and Chebyshev polynomial series. This new method to compute the correlation functions without any Limber approximation, allows us to produce and discuss maps of the correlation function directly in the observable space and is a significant step towards disentangling the data from the tested models.