Remote Sensing · Environmental Health · MDPI 2026

Satellites Can Predict
How Long You'll Live

Multimodal fusion of 11 satellite and agricultural data streams predicts US county-level life expectancy — no census, no surveys — purely from features observable from orbit across 3,108 counties over 20 years.

0.631±0.013
R² · 5-fold grouped CV
1.08yr
Mean absolute error
6.16×
Night vs daytime dominance
3,108
CONUS counties · 2000–2019
450
Features · 11 modalities
Methodology
Five-Stage Analysis Pipeline
1
Satellite Data Acquisition
All environmental features extracted via Google Earth Engine Python API from 11 globally available, continuously updated satellite and agricultural data streams. Spatial aggregation to 3,108 contiguous US county boundaries using ee.Image.reduceRegions() at 250 m scale. FAO livestock densities linearly interpolated between snapshot years (2005/2010/2015). Alaska, Hawaii, and US territories excluded — 860 rows removed from the 62,540-row raw dataset, leaving 61,680 CONUS county-year observations.
Google Earth Engine MODIS · Sentinel · Landsat FAO GLW3 250 m aggregation 61,680 obs · 2000–2019
2
Feature Engineering & Quality Control
425 base features from raw modalities joined with 25 derived cross-modal features (thermal–vegetation interactions, agricultural diversity indices, soil hydrological flags) to yield 450 total. Temporal drift validated via Mann-Kendall tests (NDVI τ=0.05, p=0.62; LST τ=0.08, p=0.43). Winsorisation at 0.1th/99.9th percentiles; median imputation for missing values. Two independent HPC runs on separate cluster nodes confirmed ΔR²<0.001 — full bit-identical reproduction.
450 total features 25 derived cross-modal Mann-Kendall drift test Winsorisation Reproduced × 2 HPC runs
3
Ablation: Single-Modality Benchmarking
11 separate Random Forest models — one per data stream, plus the 25-feature engineered set — each with modality-optimised hyperparameters. County-grouped 5-fold CV (GroupKFold by FIPS) ensures all 20 years of any given county appear exclusively in training or test. The engineered 25-feature set alone reaches R²=0.581, surpassing every individual raw modality despite using only 5.5% of the feature count.
11 single-modality models GroupKFold FIPS CV No spatial leakage Engineered-only R²=0.581
4
Full Fusion Model Training & Production CV
Production Random Forest trained on 193 pruned features (257 zero-impact features dropped at threshold 0.001) with the same county-grouped 5-fold CV. Hyperparameters: n_estimators=2000; max_depth=40; max_features=0.3; max_samples=0.9; min_samples_leaf=2. Multimodal fusion reduces unexplained variance by 32% over best single-modality model (LST: R²=0.442). Moran's I=0.08 (p=0.23) confirms zero residual spatial autocorrelation.
R² = 0.631 ± 0.013 MAE = 1.08 ± 0.02 yr 32% fusion gain vs LST Moran's I = 0.08 (p=0.23)
5
TreeSHAP Mechanistic Interpretation
TreeSHAP (interventional, marginal background n=500) computed on a stratified sample of 5,000 county-year observations. Parallelised across 5 batches with joblib — over 97 CPU-hours on a 64-core Juno node. LOWESS (bandwidth 0.20–0.40 adaptive) plus Savitzky-Golay smoothing applied to SHAP dependence curves; numerical differentiation on a 300-point grid identifies inflection points. Nighttime LST: Σ|SHAP|=0.860 yr vs daytime 0.140 yr — a 6.16× ratio. SHAP rankings corroborate RF impurity importance (Spearman ρ=0.89).
TreeSHAP interventional n=5,000 sample 97+ CPU-hours LOWESS threshold detection Night/Day = 6.16× Spearman ρ=0.89
Input Data
11 Satellite & Agricultural Data Streams
🌡️
MODIS LST
MOD11A2 v6.1 · 1 km · 8-day
14 features · Day & night channels · P10/P25/P50/P75/P90 · Mean · SD · Continuous backbone of the full 20-year analysis
🌿
MODIS Vegetation
MOD13A1 v6.1 · 500 m · 16-day
14 features · NDVI & EVI · 5 percentiles · Mean · SD
🌾
USDA Cropland
NASS CDL · 30 m · Annual
135 features · 30+ crop types · Forest classes · Development intensity levels · Largest single modality
🛰️
Landsat 8/9
Collection 2 Tier 1 SR · 30 m · 16-day
84 features · 6 SR bands · NDVI/EVI/NDMI/NDWI/SAVI/BSI · 7 statistics · Available from 2013
📡
Sentinel-1 SAR
GRD IW · 10 m · 6-day
42 features · VV/VH backscatter · GLCM textures · Cloud-penetrating · From 2015
🔭
Sentinel-2
Level-2A MSI · 10–20 m · 5-day
84 features · 6 SR bands · 6 spectral indices · QA60 cloud masking · From 2015
💧
JRC Surface Water
GSW v1.4 · 30 m · Annual
6 features · Permanent & seasonal water coverage · Water-free area fraction
⛰️
Copernicus DEM
GLO-30 · 30 m · Static
7 features · Mean elevation · SD (terrain ruggedness) · P10/P25/P50/P75/P90
🌱
ESA CCI Soil Moisture
Combined v6.1 · 0.25° · Daily→Annual
7 features · Volumetric index ×100 · Mean · SD · 5 percentiles · Continuous since 1978
🐄
FAO Livestock
GLW3 · ~10 km · 2005/10/15 epochs
32 features · 8 species (cattle, pig, chicken, goat, sheep, horse, duck, buffalo) · Density · Count · SD · SE · Linearly interpolated annual
⚗️
Derived Cross-Modal
Engineered during fusion preprocessing
25 features · NCE, TVI, Forest Buffer, Impervious Heat Index and more · R²=0.581 alone — beats all raw optical modalities → Features tab
Feature Engineering
25 Derived Cross-Modal Features

Computed during combined-model preprocessing only — not used in single-modality ablation models. These capture synergistic interactions between sensors that no individual stream can encode alone. The full set achieves R²=0.581 alone, surpassing all 84-feature optical modalities and the 135-feature USDA CDL (R²=0.433), despite using 5.5% of the feature count.

#Feature NameConstructionMechanistic RationaleDomain

425 base features + 25 derived = 450 total.

Results
Key Findings Across Model & SHAP Analysis
Explained Variance · Production CV
0.631±0.013
R² on geographically held-out counties. Standard random CV would inflate by +0.04–0.07; this is a genuine lower bound. ≈78% of IHME census-model performance with zero census inputs.
Mean Absolute Error
1.08yr
63% of counties within ±1 year; 90% within ±2 years. Core brackets 74–79 yr (2,127 counties, 68% of sample) achieve MAE<1 year. Best county: 0.13 yr MAE.
Multimodal Fusion Gain
32%
Reduction in unexplained variance vs best single modality (MODIS LST: R²=0.442). Even SAR (0.111), Soil (0.111), and JRC water (0.066) each add unique information.
Spatial Autocorrelation · Residuals
0.08
Moran's I (p=0.23). Confirms genuine geographic generalisation — not spatial memorisation of neighbouring county LE patterns.
Engineered Features Alone
0.581
25 physics-grounded features beat every raw modality: USDA CDL (135 feat, R²=0.433), Landsat (84 feat, R²=0.173), Sentinel-2 (84 feat, R²=0.079).
Hardest Counties · Spectral Shadow
7.1yr MAE
9 Indigenous reservation counties (LE bracket 60–70 yr). IHS underfunding and historical trauma are invisible to any radiometric sensor. Excluding these 9 counties, overall MAE drops to 0.87 yr.
The Nighttime Thermal Paradox

The most unexpected finding: chronic minimum overnight cooling opportunity — not extreme daytime heat — is the dominant predictor of county-level longevity. NCE (Nighttime Cooling Efficiency) is the single most important feature in the entire 450-feature dataset.

🌙 Nighttime LST · 7 channels
0.860
cumulative mean |SHAP| · years
6.16×
dominance
☀️ Daytime LST · 6 channels
0.140
cumulative mean |SHAP| · years
SHAP Feature Hierarchy · Top 12

Ablation Summary
ModalityFeaturesR² (mean ± SD)MAE (yr)
Combined (193 pruned) 193 0.631 ± 0.013 1.08
Combined (full 450) 450 0.623 ± 0.013 1.09
Engineered ✦ 250.581 ± 0.0151.16
MODIS LST140.442 ± 0.0211.36
USDA CDL1350.433 ± 0.0111.38
FAO Livestock420.396 ± 0.0121.44
MODIS NDVI/EVI140.350 ± 0.0291.49
Sentinel-2840.079 ± 0.0041.85
JRC Water60.066 ± 0.0091.88

✦ 25 engineered features surpass all 84-feature optical modalities at 5.5% of the feature count.

Policy Implications
Four Satellite-Observable Policy Thresholds

SHAP derivative analysis (LOWESS + Savitzky-Golay + numerical differentiation on a 300-point grid) identifies sharp inflection points where the direction or rate of health impact changes qualitatively. Every threshold is directly measurable from existing freely available satellite archives.

🌡️
≈ 7.3°C
Nighttime LST Benefit-to-Harm Threshold
Counties where nighttime LST 10th percentile exceeds ≈7.3°C (transitioning to harm above 10–12°C) are denied the overnight cooling window essential for cardiovascular and immune recovery. When minimum overnight temperatures close this physiological window, elevated resting heart rate, suppressed deep sleep, and impaired immune function accumulate nightly — a chronic burden invisible to daytime heat-wave indices. NCE is the single most important model feature (mean |SHAP|=0.309 yr), carrying 6.16× the attribution of all daytime LST channels combined. Detectable annually from MODIS MOD11A2.
↳ Cooling assistance prioritisation ↳ Nocturnal heat alert systems ↳ Urban tree canopy targeting
🌲
38% attenuation
Forest Cover >20% — Heat Penalty Buffer
Deciduous forest cover above 20% attenuates the heat-driven LE penalty by ≈38% in the hottest temperature quartile (Q4). The buffering is temperature-dose-dependent: negligible in cool quartiles Q1–Q2, rising to Δ=0.007 yr in Q4. Zero effect in cool counties — the buffering only activates under thermal stress above ≈20°C daytime LST. Detectable annually from USDA CDL (deciduous fraction) × MODIS LST.
↳ Urban forest conservation ↳ Nature-based climate adaptation ↳ Ecosystem co-benefit accounting
🐄
≈ 318 head/km²
Cattle Density — CAFO Oversight Inflection
A qualitative transition occurs at ≈318 head/km², with saturation at ≈1,899 head/km². Negative LE associations emerge at moderate densities, consistent with diffuse nitrogen loading and groundwater contamination. The current EPA CAFO threshold focuses on >1,000 animal units — empirical satellite data suggests earlier regulatory attention is warranted. Because this threshold is annually updatable from FAO GLW3 × USDA CDL, it provides a potential remote-sensing redline for triggering ground-level inspections.
↳ Tiered CAFO oversight ↳ Nitrogen loading monitoring ↳ Remote-sensing inspection triggers
🌱
SM ≈ 6–8
Soil Moisture Goldilocks Zone (ESA CCI ×100)
SHAP attribution reveals a nonlinear inverted-U: peak SHAP at SM≈5.8, optimal zone SM=6–8 corresponding to agricultural field capacity. Below this range, drought stress reduces LE. Above SM=8.5 (waterlogging threshold), SHAP turns strongly negative, encoding anaerobic soil conditions and flood risk. This single global measurement (daily observations since 1978) encodes the agricultural productivity–flood risk trade-off that neither NDVI nor LST can cleanly represent alone.
↳ Irrigation optimisation ↳ Flood risk early warning ↳ Agricultural drought monitoring
Final Summary
What This Work Establishes

This study demonstrates that satellite Earth observation alone — with zero census inputs — can predict life expectancy for every contiguous US county with R²=0.631 and MAE=1.08 years. Operating on 61,680 county-year observations from 2000–2019, the model achieves ≈78% of the performance of IHME's gold-standard sociodemographic model, which relies on income, education, and healthcare data that most of the world cannot provide. Two independent HPC runs confirmed ΔR²<0.001 — full bit-identical reproducibility.

The most important scientific finding is the Nighttime Thermal Paradox: chronic minimum overnight temperature — not extreme daytime heat — is the dominant predictor. Nighttime LST features carry 6.16× the cumulative SHAP attribution of all daytime channels (0.860 yr vs 0.140 yr). NCE alone has mean |SHAP|=0.309 yr — roughly 14× larger than the equivalent daytime feature (0.022 yr). This reorients the heat-health literature from an almost exclusive focus on peak daytime temperatures toward chronic nocturnal burden.

A striking secondary result: 25 engineered cross-modal features achieve R²=0.581 alone, surpassing every individual raw modality including the 135-feature USDA CDL (R²=0.433). Physics-grounded feature construction concentrates predictive signal more efficiently than raw spectral accumulation.

Three further policy-actionable thresholds: cattle density inflects at ≈318 head/km²; deciduous forest cover >20% attenuates the heat-driven LE penalty by ≈38% in the hottest quartile; and a soil moisture Goldilocks zone at SM=6–8 encodes the agricultural productivity–flood risk trade-off.

The model's systematic over-prediction in 9 Indigenous reservation counties (MAE≈7.1 yr) is not model failure — it is a diagnostic flag for where structural determinants invisible to any satellite sensor (chronic IHS underfunding, historical trauma) vastly outweigh environmental ones. A hybrid-surveillance design augmenting these counties with locally collected vital statistics is recommended for operational deployment.

Code: github.com/albertfaiz/Multimodal_geo_fusion_FM · Data: Zenodo doi:10.5281/zenodo.19229752 · Journal: MDPI Remote Sensing