Shape-constrained regression analysis has applications in dose-response modelling, environmental risk assessment, disease screening and many other areas. Incorporating the shape constraints can improve estimation efficiency and avoid implausible results. We propose a novel method, focusing on monotone curve and surface estimation, which uses Gaussian process projections. Our inference is based on projecting posterior samples from the Gaussian process. We develop theory on continuity of the projection and rates of contraction. Our approach leads to simple computation with good performance in finite samples. The proposed projection method can also be applied to other constrained-function estimation problems, including those in multivariate settings. © 2014 Biometrika Trust.