📆 Project Period | December 2025 – January 2026 and September – December, 2026 |
👤 CIN Visiting Researcher |
Project Summary
- Applied an unsupervised Multivariate Gaussian Hidden Markov Model (K = 10 states) to a 17-year monthly time series (June 1999 – June 2016) of remote-sensing and reanalysis data covering three Mediterranean karst areas (Apulia/Italy, Classical Karst/Slovenia–Italy, Notranjska/Slovenia).
- Inferred latent ecological states from 18 features per pixel-month (climate variables + satellite-derived FAPAR and albedo) without any field-labeled training data.
- Showed that the hidden states align with independent classifications (MEDALUS desertification categories, Köppen climate zones, local meteorological station records), confirming ecological interpretability from satellite data alone.
- Quantified state transition dynamics: the transition matrix reveals which states are stable (S9: mean persistence ~826 months) and which are transient (S3: ~4 months).
- Projected the state distribution forward to December 2025 using the learned transition matrix; the system converges to its stationary distribution within ~2 years, with the dominant degraded state S2 collapsing from 46 % to ~11 %.
- The framework is fully transferable to any region with equivalent open datasets (ERA5, MODIS, CMIP6-derived indices) and scales to the Mediterranean basin at coarser resolution.
Development Tools
The core of the development is written in Python 3.11. The HMM is fully unsupervised: no labeled data is required at any stage. It is implemented using hmmlearn (Multivariate Gaussian, EM training, MPM decoding), with NumPy/SciPy for preprocessing and matrix algebra, and scikit-learn for standardization. UMAP-learn is used for state visualization. Geospatial data handling relies on xarray, rioxarray, and GDAL; spatial masking on GeoPandas/Shapely.
Development Outputs
- Code
- HMM Analysis Pipeline (private repository, to be made public upon article acceptance):
Full pipeline from raw ERA5/MODIS download to HMM training, state inference, transition matrix analysis, projection, and figure generation.
Language: Python. Key modules: dataset assembly, BIC model selection, EM training, MPM inference, validation against MEDALUS/Köppen.
- Data
- Processed monthly pixel-level dataset: 18,290 pixels × 205 months × 18 features (z-scored baseline + anomaly), covering June 1999 – June 2016 across AOIs 1–3.
- Trained HMM model parameters (transition matrix A, emission means μ and covariances Σ, initial distribution π) in serialized format.
- Publications
- Gregori F., Naylor P., Ravbar N., De Giglio M., Dubbini M.
- Conference poster:
“Hidden Markov Models for Desertification Assessment in Mediterranean Karst Ecosystems” (in preparation — target journal: Applied Geography)
Presented at the 5th ECMWF–ESA Machine Learning Workshop, 13–17 April 2026.
Project Description
- Background and Motivation
Karst landscapes cover approximately 40 % of the Mediterranean basin. They are characterized by thin soils, soluble carbonate bedrock, and highly interconnected hydrological systems, making them among the most vulnerable environments to desertification and land degradation. Standard frameworks such as MEDALUS classify the current desertification status of a pixel but cannot characterize trajectories: they reveal where a system stands, not where it is heading.
Capturing ecosystem dynamics requires a model capable of handling multivariate time series with temporal autocorrelation, producing directly interpretable parameters, and operating without labeled ground truth, because in situ ecological surveys are scarce, expensive, and spatially uneven across Mediterranean karst.
- Approach
We applied a Multivariate Gaussian Hidden Markov Model (HMM) to a spatially gridded dataset of 18,290 pixels across three Areas of Interest (AOIs):
AOI | Location | Area (km²) | Climate |
1 | Apulia, Italy | ~11,400 | Csa (hot-summer Mediterranean) |
2 | Classical Karst, Slovenia–Italy | ~520 | Cfb transitional |
3 | Notranjska Karst, Slovenia | ~1,030 | Cfb/Dfb |
Input data (~1 km, monthly, June 1999 – June 2016):
Variable | Description | Source |
Tasmax / Tasmin | Max / Min air temperature | ERA5-Land |
Pr | Precipitation | ERA5-Land |
RSDS | Shortwave radiation | ERA5-Land |
VPD | Vapor-pressure deficit | ERA5-Land |
CMI | Climatic Moisture Index | Derived |
SPEI-12 | 12-month standardized drought index | Derived |
FAPAR | Fraction absorbed PAR | MODIS MOD15A2H |
Albedo | Surface broadband albedo | MODIS MCD43A3 |
Each variable was decomposed into a 12-month rolling baseline (seasonal cycle) and a monthly anomaly (departure from the baseline), yielding 18 features per pixel-month, and then standardized to z-scores. This design ensures that the HMM operates on climate-normalized signals, capturing structural differences between states rather than the mean seasonal forcing.
The model number of states (K = 10) was selected via the Bayesian Information Criterion (BIC). Training used the Expectation-Maximization (EM) algorithm; state assignment was performed with the Most Probable Marginal (MPM) decoder, which assigns each pixel-month independently to the state with the highest marginal posterior probability.
- Results
Ecological gradient. The 10 hidden states span a continuous gradient from hyper-arid degraded conditions (S6, CMI ≈ −662 mm yr⁻¹) to mesic forest states (S4, CMI ≈ +915 mm yr⁻¹). The ordering is stable and consistent across all three AOIs, confirming that the model has recovered a single latent gradient rather than region-specific clusters.
Transition dynamics. The transition matrix is diagonally dominant (all states are self-persistent). Still, the off-diagonal structure is ecologically meaningful: the main fluxes occur between adjacent states on the gradient (S6↔S5, S4↔S3, S2→S7), not between distant extremes, consistent with smooth ecological transitions rather than abrupt regime shifts.
Validation. Although the HMM was trained without any labels: - Degraded states (S1, S2) overlap with MEDALUS “affected/fragile” classes in >85 % of pixels in AOI 1. - S1 and S2 fall within Köppen Csa in 99 % and 93 % of pixels respectively; resilient states S0 and S4 fall within Cfb/Dfb in >95 % of pixels. - Meteorological station records in AOI 2 (Godnje, Komen, Lipica, Novelo, Opatje Selo, Razdrto, Senožeče) are consistent with the emission profiles of the states assigned to those locations.
Projection. Iterating the learned transition matrix from the June 2016 state distribution forward to December 2025 reveals rapid convergence to the stationary distribution π* within approximately 2 years. The dominant degraded state S2, which occupied 46 % of pixels in 2016, converges to ~10.6 %; state S7 rises from 3 % to ~18 %. This rapid convergence indicates that year-specific predictive power is limited beyond the long-run equilibrium encoded in π*.
- Collaboration with ESA Φ-lab
The collaboration with Peter Naylor (ESA Φ-lab) contributed to the methodological design of the HMM pipeline, the interpretation of the ecological gradient, and the framing of the results for the Earth Observation community. The ESA Φ-lab context enabled access to Earth observation expertise and provided a satellite-science perspective on the validation strategy.