This work introduces the use of Gaussian processes (GPs) for the estimation and understanding of crop development and yield using multisensor satellite observations and meteo- rological data. The proposed methodology combines synergistic information on canopy greenness, biomass, soil, and plant water content from optical and microwave sensors with the atmospheric variables typically measured at meteorological stations. A com- posite covariance is used in the GP model to account for varying scales, nonstationary, and nonlinear processes. The GP model reports noticeable gains in terms of accuracy with respect to other machine learning approaches for the estimation of corn, wheat, and soybean yields consistently for four years of data across continental U.S. (CONUS). Sparse GPs allow obtaining fast and compact solutions up to a limit, where heavy sparsity compromises the credibility of confidence intervals. We further study the GP interpretability by sensitivity analysis, which reveals that remote sensing parameters accounting for soil moisture and greenness mainly drive the model predictions. GPs finally allow us to identify climate extremes and anomalies impacting crop productivity and their associated drivers.