Peaks STAC Dataset: Getting Started¶
Pictured: Eric Gilbertson measuring the elevation of Mount Ararat (Ağrı Dağı),
the highest peak in Turkey, located in far eastern Anatolia near the
Armenian and Iranian borders, July 2025.
Welcome to the "getting started" tutorial for the Peaks STAC Dataset. This notebook is a user-friendly tour of a new, machine-learning optimized dataset built from Peakbagger GNSS surveys. By moving through this notebook, you will:
- Learn why the dataset matters and what STAC provides (TLDR; discoverability, provenance, versioning).
- Load the STAC collection and items.
- Read a GeoParquet for fast analysis.
- Build a simple ML-ready table.
- Train a beginner baseline machine learning model.
- Visualize the data in two distinct ways.
Let's get started!
1. Why this dataset matters¶
Peakbagger’s original format is great for humans, but it is not optimized for programmatic access, provenance tracking, or machine learning workflows. This project turns those surveys into a durable, cloud‑native geospatial dataset that is:
- Discoverable: STAC lets you index and search peaks like any other Earth data asset.
- Provenance‑aware: You can trace each elevation back to the survey method, datum, and source.
- Versionable: Heights can change over time while the summit identity remains stable.
- ML‑ready: You get fast, columnar GeoParquet outputs that fit cleanly into data pipelines.
That means you can build reliable downstream analytics, such as risk models or hydro/terrain studies, without re‑scraping or re‑parsing the website, and with the metadata you need to trust the measurements.
2. Setup¶
First, let's identify where the data live and point our script to the data.
from pathlib import Path
NOTEBOOK_DIR = Path.cwd()
PROCESSED_DIR = NOTEBOOK_DIR.parent / "data" / "processed"
STAC_DIR = NOTEBOOK_DIR.parent / "data" / "stac"
PROCESSED_DIR.exists(), STAC_DIR.exists()
The data have been prepared in three versions: GeoJSON, GeoParquet, and STAC. These formats allow users to work in familiar tools (e.g., QGIS), scale to large-area analysis with columnar storage, and integrate the data with modern, cloud-native geospatial and machine learning ecosystems.
- STAC is a metadata and indexing standard that enables coud-native discovery, access, and machine learning integration.
- GeoJSON is helpful for lightweight visualization, manual inspection.
- GeoParquet is columnar, compressed storage designed for big geospatial analytics, enabling fast filtering, joints, and processing.
3. Load the STAC collection¶
The SpatioTemporal Asset Catalog (STAC) provides a consistent, discoverable entry point for the dataset and all its assets. A spatiotemporal asset is any file that represents information about the earth captured in a certain space and time. The STAC Specification consists of four semi-independent specifications:
- STAC Item is the core atomic unit, representing a single asset as a GeoJSON feature plus datetime and links
- STAC Catalog is a simple, flexible JSON file of links that provides structure to organize and browser STAC ITems.
- STAC Collection is an extension of the STAC Catalog with additional information like extents, license, keywords, and providers that descrime the STAC ITems that fall within the Collection.
- STAC API provides a RESTful endpoint that enables search of STAC Items.
A STAC collection can be programmatically discovered, validated, and queried using standard tooling such as pystac.
import pprint
import pystac
# Load STAC collection
collection = pystac.Collection.from_file(str(STAC_DIR / "collection.json"))
collection.to_dict().keys()
# Basic metadata
print("Collection ID:", collection.id)
print("Description:", collection.description)
# Spatial extent
print(f"Spatial extent (bboxes): {collection.extent.spatial.bboxes}")
# Temporal extent
print(f"Temporal extent (intervals): {collection.extent.temporal.intervals}")
# Load items
items = list(collection.get_all_items())
print(f"Number of items: {len(items)}")
item = items[0]
print("Example item:")
print(f" Item ID: {item.id}")
print(f" Datetime: {item.datetime}")
print(f" Geometry type: {item.geometry['type']}")
print(f" Assets: {list(item.assets.keys())}")
One of the benefits of using STAC is that you can query each Item's metadata without needing to downlaod the actual data file. For example, we can easily identify all peaks that fall within the Northh Cascades".
target_range = "North Cascades"
north_cascades = [
x for x in items
if target_range in (x.properties.get("ranges"))
]
print(f"Peaks in {target_range}: {len(north_cascades)}")
# Short preview
for it in north_cascades[:20]:
props = it.properties
name = props.get("name", it.id)
peak_id = props.get("peak_id")
state = props.get("state_province")
country = props.get("country")
lat = props.get("latitude")
lon = props.get("longitude")
print(
f"- {name} "
f"(peak_id={peak_id}, {state}, {country}) "
f"@ ({lat}, {lon})"
)
4. Load the GeoParquet¶
GeoParquet provides a scalable, analytics-optimized dataset designed for efficient filtering, aggregation, and analysis. GeoParquet is a columnar storage format that extends Apache Parquet with standardized geospatial metadata, enabling geometry-aware operations while retaining the performance benefits of columnar data layouts.
In GeoParquet, each row represents a geographic feature, with attributes stored as columns and geometries encoded in a standardized binary format. Users can efficiently query and subset the dataset by spatial, categorical, or numeric attributes (e.g., country, range, elevation uncertainty) without loading the full dataset into memory.
GeoParquet is well suited for:
- Large-scale geospatial analytics and joins
- Cloud-native and distributed processing (e.g., DuckDB, Spark, BigQuery)
- Feature engineering and model training pipelines that benefit from fast, selective reads
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
parquet_path = PROCESSED_DIR / "peaks.parquet"
peaks = gpd.read_parquet(parquet_path)
print(f"Rows: {len(peaks)}")
print(f"Columns: {peaks.columns}")
peaks.head(5)
This data format makes it easy to run quick analyes. For instance, we may be interested in comparing true isolation across peaks. Isolation refers to the distance of a summit to the nearest higher land. It represents the radius of the area where the peak is the highest point. For instance, the isolation for Mt. Everest is undefined, since there is no nearest higher land, but every other summit can be asigned an isolation value.
# Example quick analyses
peaks.true_isolation_mi.describe()
The average isolation value in this dataset is 28.7 miles, with a maximum of 731.1 miles.
plt.figure(figsize=(12, 4))
peaks["true_isolation_mi"].dropna().plot(kind="hist", bins=30, color="lightgray")
plt.title("Isolation Distribution")
plt.xlabel("Isolation (miles)")
plt.ylabel("Count")
plt.tight_layout()
Let's examine which peak is uniquely isolated in this dataset.
peaks[peaks.true_isolation_mi > 700].name
Eight peaks in the dataset appear to be missing datum values, while 6 are missing a recorded elevation source method. This information allows us to easily correct for missing fields.
plt.figure(figsize=(12, 4))
peaks["elevation_m"].dropna().plot(kind="hist", bins=30, color="lightgray")
plt.title("Elevation Distribution")
plt.xlabel("Elevation (meters)")
plt.ylabel("Count")
plt.tight_layout()
According to the elevation_m field, a majority of peaks fall between approximately 2000-3000 meters high, with only a handful surpassing this value. Let's see which peak is tallest.
peaks[peaks.elevation_m > 5000].name
# Quick data-quality snapshot for important fields
required = [
"peak_id",
"name",
"country",
"feature_type",
"latitude",
"longitude",
"elevation_m",
"datum",
"source_method"
]
missing = peaks[required].isna().sum().sort_values(ascending=False)
missing
5. Visualization: Global Distribution¶
Using this geoparquet file, we can easily visualize patterns in the data. For instance, we may be interested in:
- Where are these peaks located, geographically?
- Where are the highest peaks in the dataset?
# Load a simple world basemap
world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
# Create a figure
fig, ax = plt.subplots(figsize=(10, 5))
# Plot the basemap
world.plot(
ax=ax,
color="lightgray",
edgecolor="white"
)
# Plot peak locations, colored by elevation
scatter = ax.scatter(
peaks["longitude"],
peaks["latitude"],
c=peaks["elevation_m"],
cmap="viridis",
s=8,
alpha=0.8
)
# Add a colorbar
plt.colorbar(scatter, ax=ax, label="Elevation (m)")
# Add labels and title
ax.set_title("Peak locations colored by elevation")
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
plt.tight_layout()
7. Simple Machine Learning Example¶
To establish a simple baseline, we can train a lightweight model that predicts peak elevation using only geographic location (lat/lon) and a small number of optional, commonly available attributes.
This baseline is intentially minimal: it helps us understand how much predictive power comes from basic spatial context alone, and provides a reference point for evaluating more sophisticated models later.
from sklearn.model_selection import train_test_split
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.metrics import mean_absolute_error
from sklearn.ensemble import RandomForestRegressor
# Select features and target
feature_cols = [
"latitude",
"longitude",
"country",
"clean_prominence_m",
"isolation_distance_km",
]
target_col = "elevation_m"
ml_df = peaks[feature_cols + [target_col]].copy()
X = ml_df[feature_cols]
y = ml_df[target_col]
# Set up train / test split
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)
# Preprocessing
numeric_features = [
"latitude",
"longitude",
"clean_prominence_m",
"isolation_distance_km",
]
categorical_features = ["country"]
preprocess = ColumnTransformer(
transformers=[
(
"numeric",
SimpleImputer(strategy="median"),
numeric_features,
),
(
"categorical",
Pipeline(
steps=[
("imputer", SimpleImputer(strategy="most_frequent")),
("onehot", OneHotEncoder(handle_unknown="ignore")),
]
),
categorical_features,
),
]
)
# Baseline model
model = RandomForestRegressor(
n_estimators=200,
random_state=42,
n_jobs=-1,
)
baseline_model = Pipeline(
steps=[
("preprocess", preprocess),
("model", model),
]
)
# Train and evaluate
baseline_model.fit(X_train, y_train)
y_pred = baseline_model.predict(X_test)
mae = mean_absolute_error(y_test, y_pred)
print(f"Baseline MAE: {mae:.2f} meters")
baseline_model
We train a simple baseline model using a random forest regressor to predict peak elevation from geographic location (latitude and longitude) and a small set of optional attributes (country, prominence, and isolation). The model includes basic preprocessing: median imputation for numeric features, and one-hot encoding for categorical variables. We evaluate our model using mean absolute error (MAE).
The resulting baseline achieves an MAE of ~368 meters, which is expectedly coarse for an elevation prediction task. While this performance is not sufficient for any applications, it provides a useful starting point and reference for incorporating more representative and predictive inputs, such as digital elevation models (DEMs), terrain derivatives, or remote sensing features, in subsequent models.
8. Conclusion¶
This notebook demonstrates how structuring geospatial data in machine learning-ready formats like STAC and GeoParquet makes it easy to discover, inspect, visualize, and model the data using standard tools. By separating data organization and metadata from modeling complexity, these formats lower the barrier to experimentation and enable more scalable, reproducible geospatial machine-learning workflows.