Vegetation index time series from Landsat and Sentinel-2 have great potential for following the dynamics of ecosystems and are the key to develop essential variables in the realm of biodiversity. Unfortunately, the removal of pixels covered mainly by clouds reduces the temporal resolution, producing irregularity in time series of satellite images. We propose a Bayesian approach based on a harmonic model, fitted on an annual base. To deal with data sparsity, we introduce hierarchical prior distribution that integrate information across the years. From the model, the mean and standard deviation of yearly Ecosystem Functional Attributes (i.e., mean, standard deviation, and peak’s day) plus the inter-year standard deviation are calculated. Accuracy is evaluated with a simulation that uses real cloud patterns found in the Peneda-Gêres National Park, Portugal. Sensitivity to the model’s abrupt change is evaluated against a record of multiple forest fires in the Bosco Difesa Grande Regional Park in Italy and in comparison with the BFAST software output. We evaluated the sensitivity in dealing with mixed patch of land cover by comparing yearly statistics from Landsat at 30m resolution, with a 2m resolution land cover of Murgia Alta National Park (Italy) using FAO Land Cover Classification System 2.