In this article, we propose an approach based on Gaussian Processes (GP) for large scale land cover pixel-based classification with Sentinel-2 satellite image time-series (SITS). We used a sparse approximation of the posterior combined with variational inference to learn the GP's parameters. We applied stochastic gradient descent and GPU computing to optimize our GP models on massive data sets. The proposed GP model can be trained with hundreds of thousands of samples, compared to few thousands for traditional GP methods. Moreover, we included the spatial information by adding the geographic coordinates into the GP's covariance function to efficiently exploit the spatiospectro-temporal structure of the SITS. We ran experiments with Sentinel-2 SITS of the full year 2018 over an area of 200 000 km 2 (about 2 billion pixels) in the south of France, which is representative of an operational setting. Adding the spatial information significantly improved the results in terms of classification accuracy. With spatial information, GP models have an overall accuracy of 79.8. They are more than three points above Random Forest (the method used for current operational systems) and more than one point above a multi-layer perceptron. Compared to a Transformer-based model (which provides state of the art results in the literature, but are not applied in operational systems), GP models are only one point below.