Terrestrial ecosystems exchange carbon fluxes (CFs) with the atmosphere. The three main CFs are gross primary production (GPP), terrestrial ecosystem respiration (TER), and net ecosystem exchange (NEE). GPP stands for the total amount of carbon that is fixed by plants through photosynthesis, whereas TER refers to carbon that is released into the atmosphere by respiration, and NEE is the difference between GPP and TER. As CFs are related to CO₂ assimilation and biomass production, they can be used to evaluate the state of ecosystems. This work aims to quantify long-term changes of CFs using EO data of terrestrial ecosystems at the regional scale in mainland Spain. Eight days' worth of CF images (and uncertainties), at a 1 km scale for the 2002–2023 period, were obtained from a daily global monitoring product elaborated on by the authors using machine learning tools. A data-driven approach based on a multi-output Gaussian process regression algorithm, blending MODIS products and in situ eddy covariance data, allowed us to jointly estimate GPP, TER, and NEE, preserving their physical relationship. A novel methodology based on nonlinear Bκ embeddings (NLEs) that can analyze time series at various temporal scales was applied. CF time series are nonlinearly embedded into an encompassing mathematical structure that depends on a continuous parameter k, which gives the NLE method a great flexibility and enormous potential to analyze the series at various temporal scales (just changing the k value). Its potential to map subtle, long-term changes has been demonstrated using NDVI data. NLE is able to quantify vegetation changes with very low uncertainty and high statistical significance by means of the slope Q of the trend (or inter-annual component of the time series). A map of statistically significant areas showing vegetation greening (positive Q) and browning (negative Q) is produced, showing the state of ecosystems in this highly vulnerable area.