A machine-learning pipeline that ranks undiscovered geothermal energy sites across the Western United States using real geophysical data.
Where in the Western US are the most promising undiscovered geothermal sites, and how confident are we?
The pipeline produces a calibrated probability score for every 0.1° grid cell (~11 km), SHAP explanations per site, and an interactive dashboard so a geologist can explore predictions and understand why a location scored high.
geoprospect/
├── app.py # Dash EDA dashboard (12-tab visual explorer)
├── requirements.txt
│
├── src/
│ ├── feature_engineering.py # Builds feature matrix (0.1° grid × 11 features)
│ ├── train.py # XGBoost + spatial CV + SHAP + calibration + MLflow
│ ├── train_cnn.py # 2D CNN on 5×5 geophysical patches + MLflow
│ └── download_known_sites.py # One-time: downloads known geothermal site labels
│
├── data/
│ ├── heat/ # SMU equilibrium HF + BHT CSVs
│ ├── gravity-geotiff/ # USGS Bouguer gravity XYZ
│ ├── Faults/ # USGS Quaternary faults SHP + magnetic XYZ
│ ├── known_sites/ # NREL geothermal site labels (CSV + GeoJSON)
│ └── processed/
│ ├── feature_matrix.parquet # Output of feature_engineering.py
│ ├── scored_grid.parquet # Output of train.py (XGBoost scores + SHAP)
│ └── scored_grid_cnn.parquet # Output of train_cnn.py (CNN scores)
│
├── models/
│ ├── xgb_model.pkl # XGBoost model + calibrator
│ └── cnn_model.pt # CNN weights + calibrator + normalization stats
│
└── mlruns/ # MLflow experiment tracking (local)
All data is free and publicly available.
| Dataset | Source | Records |
|---|---|---|
| Equilibrium Heat Flow | SMU via GDR/OpenEI | ~8,400 filtered (15,655 raw) |
| BHT Heat Flow | SMU via GDR/OpenEI | ~3,900 filtered (38,717 raw) |
| Bouguer Gravity | USGS mrdata | 233,661 pts in WUS |
| Magnetic Anomaly | USGS mrdata | 921,838 pts in WUS |
| Quaternary Faults | USGS Earthquake Hazards | 109,670 segments (51,127 Holocene) |
| Known Geothermal Sites | NREL (operating + developing + low-temp) | 1,553 WUS sites |
Geographic scope: Western US — lon −125 to −102, lat 31 to 49.
One row per 0.1° grid cell (~41,400 cells total). Aggregations use a 0.5° radius (~55 km).
| Feature | Description | Source |
|---|---|---|
hf_mean |
Mean equilibrium heat flow within 0.5° (mW/m²) | SMU HF |
hf_max |
Max heat flow within 0.5° | SMU HF |
hf_n |
Number of HF measurements within 0.5° | SMU HF |
bht_mean |
Mean BHT-derived heat flow within 0.5° | SMU BHT |
gravity_mean |
Nearest Bouguer gravity value (mGal) | USGS Gravity |
magnetics_mean |
Nearest magnetic anomaly (nT) | USGS Magnetics |
fault_length_total |
Total fault length within 0.5° (km) | USGS Faults |
fault_length_holocene |
Holocene-equivalent fault length within 0.5° (km) | USGS Faults |
dist_nearest_fault |
Distance to nearest fault (km) | USGS Faults |
dist_nearest_holocene |
Distance to nearest Holocene-equivalent fault (km) | USGS Faults |
Holocene definition: The USGS Quaternary fault database uses latest Quaternary (~last 15,000 yrs) and historic (~last 200 yrs) as the Holocene-equivalent age classes — not the string "Holocene". Feature engineering matches these explicitly.
Labeling: A cell is is_geothermal = 1 if a known site exists within 0.05° (~5.5 km). This tight radius (~985 positives, 2.4%) reflects the actual discovery problem — finding new sites — rather than regional classification.
dist_to_nearest_site is computed but never used as a model feature — only for filtering novel candidates post-inference.
Gradient-boosted trees trained on the full 11-feature tabular matrix.
- Imbalance handling:
scale_pos_weight = 41(ratio of negatives to positives) - Validation: 5-fold spatial
GroupKFoldwith 3°×3° geographic blocks - Calibration: Isotonic regression fit on out-of-fold predictions
- Explainability: SHAP values computed per cell via
TreeExplainer - Tracking: Logged to MLflow as
xgb-spatial-cv
OOF performance: ROC-AUC = 0.753 | PR-AUC = 0.057 (random baseline = 0.024)
Top features by gain: hf_max (0.178) → dist_nearest_fault (0.155) → dist_nearest_holocene (0.099) → hf_mean (0.091)
A small 2D CNN that reads 5×5 spatial patches (~55 km wide) of 4 geophysical channels and predicts geothermal probability from local spatial context rather than point features.
Architecture: 2 conv layers (3×3) → BatchNorm → ReLU → Global Average Pooling → FC → sigmoid. ~5,361 parameters.
4 channels: hf_mean, gravity_mean, magnetics_mean, fault_length_total — channels with spatial structure. Distance/count features are excluded (already aggregated, no spatial pattern to convolve over).
Design choices:
- Reflection padding at map edges (avoids "ocean" zero-padding artifact)
- Per-channel z-score normalization (heat flow ~30–200, magnetics ~−1,600 to +2,700)
- Same 3°×3° spatial GroupKFold as XGBoost — fair comparison
BCEWithLogitsLosswithpos_weight = 41for imbalance- Early stopping (patience = 8 epochs on val AP)
- Tracking: Logged to MLflow as
cnn-spatial-cv
OOF performance: ROC-AUC = 0.761 | PR-AUC = 0.067 (random baseline = 0.024)
A cell-level mean of calibrated XGBoost and CNN scores — no separate training required. Computed at dashboard load time.
Top-1% prospect agreement:
- Both models flag: 395 cells (22.8%)
- XGBoost only: 1,338 cells
- CNN only: 820 cells
The low overlap reflects genuinely different inductive biases — XGBoost learns from point feature combinations while CNN learns from spatial neighbourhood patterns. The ensemble surfaces cells that both models agree on, yielding higher-confidence candidates.
Standard random k-fold leaks because adjacent cells share geology. The WUS is divided into 3°×3° geographic tiles; GroupKFold(n_splits=5) ensures no tile appears in both train and validation. All reported metrics are OOF (out-of-fold) — honest estimates of performance on unseen geographies.
python -m venv .venv
source .venv/bin/activate # Windows: .venv\Scripts\activate
pip install -r requirements.txtpython src/feature_engineering.py
# Runtime: ~3-5 min
# Output: data/processed/feature_matrix.parquetpython src/train.py
# Runtime: ~1-2 min
# Output: data/processed/scored_grid.parquet, models/xgb_model.pklpython src/train_cnn.py
# Runtime: ~5-10 min (CPU) / ~1 min (GPU)
# Output: data/processed/scored_grid_cnn.parquet, models/cnn_model.ptpython app.py
# Open http://localhost:8050mlflow ui
# Open http://localhost:5000
# Experiment: geothermal-prospect-scoring
# Runs: xgb-spatial-cv, cnn-spatial-cv| Predictions | Faults |
|---|---|
![]() |
![]() |
| Heat Flow | Gravity |
|---|---|
![]() |
![]() |
| Feature Grid | Correlation |
|---|---|
![]() |
![]() |
| Tab | What it shows |
|---|---|
| Heat Flow | SMU equilibrium HF measurements, quality filter, histogram |
| BHT | Bottom-hole temperature derived heat flow |
| Gravity | Bouguer gravity anomaly grid |
| Magnetics | Magnetic anomaly grid |
| Faults | Quaternary fault traces by age category (Holocene highlighted) |
| Known Sites | NREL operating, developing, low-temp sites |
| Overlay | All layers toggled on a single map |
| KDE Analysis | Geothermal sites vs random WUS background — separability check |
| Coverage | Data density per 0.2° cell (heat flow sparse, gravity dense) |
| Correlation | Feature correlation matrix + gravity vs heat flow scatter |
| Feature Grid | Any feature visualized as a map with geothermal site overlay |
| Predictions | XGBoost / CNN / Ensemble scores with threshold slider, novel-only filter |
scored_grid.parquet (XGBoost):
| Column | Description |
|---|---|
score |
Calibrated geothermal probability |
score_raw |
Uncalibrated XGBoost output |
is_novel |
1 if > 50 km from any known site |
shap_* |
SHAP attribution per feature |
scored_grid_cnn.parquet (CNN):
| Column | Description |
|---|---|
cnn_score |
Calibrated CNN probability |
cnn_score_raw |
Uncalibrated sigmoid output |
is_novel |
1 if > 50 km from any known site |
Novel candidates = cells where is_novel == 1 AND is_geothermal == 0, ranked by score. These are the actionable predictions for follow-up exploration.
| Signal | Geothermal meaning |
|---|---|
High hf_max |
Direct evidence of elevated crustal heat flux |
Low dist_nearest_fault |
Proximity to fracture network → fluid pathway |
Low dist_nearest_holocene |
Active fault → permeable, likely fluid-bearing zone |
High hf_mean |
Sustained regional heat anomaly |
Negative gravity_mean |
Low-density, thermally expanded crust (Basin & Range) |
The Basin & Range province (Nevada, western Utah, Arizona) dominates high-scoring regions because it combines all five signals.
- Heat flow has 28% cell-level missingness (sparse measurement coverage). Imputed with column median before training.
- Gravity and magnetics features are nearest-neighbour lookups (dense grids), not kriged surfaces.
- The "negative" class contains unmapped geothermal sites — PU (Positive-Unlabeled) learning would be more rigorous than
scale_pos_weight. - CNN uses only 4 channels;
dist_nearest_faultand count features are excluded because they don't have spatial raster structure to convolve over. - Fold 2 CNN PR-AUC (0.049) is notably weaker than other folds — that geographic block has fewer positives (155) and likely less signal overlap with training regions.





