Download this notebook

Step 4: Add Metadata and Analyze Spatial Patterns

This notebook turns metric results into spatial summaries. You will use PGA and FAS as two separate examples, so you can see how the same spatial tests can highlight different model-performance patterns for different metrics.

Imports

These functions prepare metric fields and calculate spatial summaries.

from spatial_vtk.config.notebook import notebook_timer, register_svtk_cell_timer

with notebook_timer():
    import os

    os.environ.setdefault("LOKY_MAX_CPU_COUNT", "1")

    from pathlib import Path
    import pandas as pd
    from IPython.display import Markdown, display

    from spatial_vtk.config import SpatialVTKConfig
    from spatial_vtk.config.labels import metric_display_name
    from spatial_vtk.io import read_config_table, write_output_tables
    from spatial_vtk.spatial.calculate import (
        bootstrap_contrast_table,
        build_distance_bin_summary,
        build_metric_field,
        build_station_feature_table,
        center_field_by_event,
        compute_global_morans_i,
        compute_pca_spatial_modes,
        moran_result_to_frame,
        normalize_metrics_table,
        run_residual_feature_clustering,
        summarize_station_bias,
    )
    from spatial_vtk.spatial.map import plot_residual_grid, plot_station_bias_map
    from spatial_vtk.spatial.map.pca import plot_pca_summary
    from spatial_vtk.spatial.plot import plot_distance_correlation_by_metric, plot_geology_contrast
    register_svtk_cell_timer()

Run time: 1.42 s

Configuration

Load the config and read the spatial-statistics settings from the tutorial run scenario.

# Use the repository root so paths match the public source checkout.
from spatial_vtk.tutorials import tutorial_root
repo_root = tutorial_root()
# Metadata table paths are project-relative.
import os
os.chdir(repo_root)
config_path = repo_root / "data/examples/configuration/example_spatial_vtk_config.yaml"

# Load the tutorial run scenario and make it the active config for later package calls.
cfg = SpatialVTKConfig.from_file(config_path, run_scenario="tutorial").activate()

# Use the native ln metric results calculated in Step 3.
figure_dir = cfg.path("outputs.figures")
figure_dir.mkdir(parents=True, exist_ok=True)
spatial_metrics = ("PGA", "FAS")
spatial_metric_labels = {"PGA": "Peak acceleration (PGA)", "FAS": "Fourier amplitude spectrum (FAS), 3 s"}
spatial_value_column = "ln_residual"

Run time: 13.0 ms

Prepare Metric-Specific Spatial Fields

Spatial statistics use one value per event-station observation, plus station and event coordinates. Here you will build those fields separately for PGA and FAS. Use Z-component PGA in the 1–2 s passband and FAS at 3 s from the independent 1 Hz lowpass calculation. Each field has one observation per event and station; periods, components and passbands are not pooled.

# Use the native ln metric results calculated in Step 3.
from spatial_vtk.io import load_output_table
metrics = load_output_table("metrics_enriched")
metrics = metrics.loc[
    metrics["component"].eq("Z") & (
        (metrics["metric"].eq("PGA") & metrics["passband"].eq("1-2 sec")) |
        (metrics["metric"].eq("FAS") & metrics["period_s"].eq(3.0) & metrics["passband"].eq("lowpass 1 Hz"))
    )
].copy()
assert not metrics.duplicated(["event_id", "station", "metric"]).any()

# Read the configured site/geology metadata table.
site_metadata = read_config_table("paths.site_metadata")

# Normalize metric columns so spatial-statistics helpers see consistent names.
spatial_ready = normalize_metrics_table(metrics)

spatial_products = {}
for metric_name in spatial_metrics:
    display(Markdown(f"### {spatial_metric_labels[metric_name]}"))

    # Build one event-station field for this metric using ln(observed/synthetic) residuals.
    field = build_metric_field(spatial_ready, metric=metric_name, value_column=spatial_value_column)

    # Remove each event mean before comparing persistent station patterns.
    centered = center_field_by_event(field)

    # Summarize each station's mean event-centered residual for this metric.
    station_bias = summarize_station_bias(centered)

    spatial_products[metric_name] = {
        "field": field,
        "centered": centered,
        "station_bias": station_bias,
    }

    display(
        pd.DataFrame(
            {
                "Output": ["Metric field", "Event-centered field", "Station bias"],
                "Rows": [len(field), len(centered), len(station_bias)],
                "Events": [field["event_id"].nunique(), centered["event_id"].nunique(), ""],
                "Stations": [field["station"].nunique(), centered["station"].nunique(), station_bias["station"].nunique()],
            }
        )
    )
    display(station_bias.head())

metric_field = pd.concat([item["field"] for item in spatial_products.values()], ignore_index=True)
event_centered_residuals = pd.concat([item["centered"] for item in spatial_products.values()], ignore_index=True)
station_bias = pd.concat([item["station_bias"].assign(metric=metric_name) for metric_name, item in spatial_products.items()], ignore_index=True)

# Save the core spatial-statistics handoff tables for later notebooks.
write_output_tables(
    metric_field=metric_field,
    event_centered_residuals=event_centered_residuals,
    station_bias=station_bias,
)

Peak acceleration (PGA)

	Output	Rows	Events	Stations
0	Metric field	115	5	30
1	Event-centered field	115	5	30
2	Station bias	30		30

	station	lat	lon	n_events	mean_centered	median_centered	std_centered	sem_centered	abs_mean_centered	bias_zscore
7	CRN	33.876	-117.561	3	-1.989584	-2.043541	0.741681	0.428210	1.989584	-4.646286
17	IVY	33.760	-117.477	4	1.130737	1.021022	0.442097	0.221049	1.130737	5.115333
6	CRF	33.874	-118.039	1	0.860199	0.860199	NaN	NaN	0.860199	NaN
16	GSA	34.137	-118.128	3	0.792946	0.489647	0.717965	0.414517	0.792946	1.912937
13	FON	34.100	-117.439	4	-0.791733	-0.720329	0.371544	0.185772	0.791733	-4.261847

Fourier amplitude spectrum (FAS), 3 s

	Output	Rows	Events	Stations
0	Metric field	140	5	30
1	Event-centered field	140	5	30
2	Station bias	30		30

	station	lat	lon	n_events	mean_centered	median_centered	std_centered	sem_centered	abs_mean_centered	bias_zscore
2	BLC	34.24400	-118.67300	4	1.069910	0.818259	0.834631	0.417315	1.069910	2.563792
7	CRN	33.87600	-117.56100	4	-0.882931	-0.581021	0.685306	0.342653	0.882931	-2.576748
27	MWC	34.22400	-118.05800	3	-0.868699	-0.499090	0.692043	0.399551	0.868699	-2.174189
20	LCG	34.00033	-118.37794	5	0.834246	0.490915	0.796891	0.356381	0.834246	2.340884
26	MLS	34.00500	-117.56200	5	0.816817	0.736142	1.236567	0.553009	0.816817	1.477040

{'metric_field': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/metric_field.parquet'),
 'event_centered_residuals': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/event_centered_residuals.parquet'),
 'station_bias': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/station_bias.parquet')}

Run time: 189.3 ms

Station Bias Maps

These maps show the mean event-centered residual at each station. Positive values mean the observed amplitudes are larger than the synthetic amplitudes on average for that metric.

for metric_name, products in spatial_products.items():
    display(Markdown(f"### {spatial_metric_labels[metric_name]}"))

    # Map the mean event-centered residual at each station for this metric.
    station_bias_fig = plot_station_bias_map(
        products["station_bias"],
        title=f"{spatial_metric_labels[metric_name]} Station Bias",
        value_col="mean_centered",
        value_label="Mean event-centered ln(obs/syn)",
        add_basemap=True,
        showfig=True,
        savefig=True,
        outpath=figure_dir / f"step_04_{metric_name.lower()}_station_bias.png",
    )

Peak acceleration (PGA)

Fourier amplitude spectrum (FAS), 3 s

Run time: 3.13 s

Residual Grid Maps

A residual grid gives you a quick spatial overview of where residuals are broadly positive or negative. These examples use the same event-centered residual field as the station-bias maps.

for metric_name, products in spatial_products.items():
    display(Markdown(f"### {spatial_metric_labels[metric_name]}"))

    # Grid event-centered residuals onto lon/lat cells for a broader spatial view.
    residual_grid_fig = plot_residual_grid(
        products["centered"],
        lon_col="lon",
        lat_col="lat",
        value_col="field_centered",
        cell_size_deg=0.05,
        title=f"{spatial_metric_labels[metric_name]} Residual Grid",
        add_basemap=True,
        showfig=True,
        savefig=True,
        outpath=figure_dir / f"step_04_{metric_name.lower()}_residual_grid.png",
    )

Peak acceleration (PGA)

Fourier amplitude spectrum (FAS), 3 s

Run time: 3.52 s

Spatial Correlation Tests

Moran’s I and distance-bin correlations help you see whether residuals cluster in space. Run them separately for each metric so the results are easier to interpret, then plot the distance-binned correlations together to compare the spatial pattern.

moran_tables = {}
distance_bin_tables = {}

for metric_name, products in spatial_products.items():
    display(Markdown(f"### {spatial_metric_labels[metric_name]}"))

    # Compute global Moran's I for this metric's station-bias field.
    moran_result = compute_global_morans_i(products["station_bias"])

    # Convert the Moran result object into a table for saving and display.
    morans_i = moran_result_to_frame(moran_result).assign(metric=metric_name)

    # Summarize residual similarity as a function of station separation distance for this metric.
    distance_bins = build_distance_bin_summary(products["centered"]).assign(metric=metric_name)

    moran_tables[metric_name] = morans_i
    distance_bin_tables[metric_name] = distance_bins
    display(morans_i)
    display(distance_bins.head())

morans_i = pd.concat(moran_tables.values(), ignore_index=True)
distance_bins = pd.concat(distance_bin_tables.values(), ignore_index=True)

# Plot correlation as a function of station separation distance for the selected metrics.
correlation_distance_fig = plot_distance_correlation_by_metric(
    distance_bins,
    significance_df=morans_i,
    title="Spatial Correlation by Distance",
    showfig=True,
    savefig=True,
    outpath=figure_dir / "step_04_spatial_correlation_distance.png",
)

# Save the spatial-correlation result tables.
write_output_tables(
    morans_i=morans_i,
    permutation_moran=morans_i,
    distance_bin_correlations=distance_bins,
)

Peak acceleration (PGA)

	n	k	moran_i	p_two_sided	permutations	metric
0	30	2	0.122716	0.46	99	PGA

	distance_start_km	distance_end_km	distance_center_km	pair_count	event_count	mean_pair_correlation	mean_semivariance	metric
0	0.0	20.0	10.0	166	5	0.051911	0.371014	PGA
1	20.0	40.0	30.0	314	5	0.076451	0.312521	PGA
2	40.0	60.0	50.0	340	5	-0.051039	0.482357	PGA
3	60.0	80.0	70.0	235	5	-0.141052	0.629634	PGA
4	80.0	100.0	90.0	189	5	-0.167490	0.621686	PGA

Fourier amplitude spectrum (FAS), 3 s

	n	k	moran_i	p_two_sided	permutations	metric
0	30	2	-0.113995	0.45	99	FAS

	distance_start_km	distance_end_km	distance_center_km	pair_count	event_count	mean_pair_correlation	mean_semivariance	metric
0	0.0	20.0	10.0	248	5	0.007315	0.872195	FAS
1	20.0	40.0	30.0	487	5	-0.062901	0.928965	FAS
2	40.0	60.0	50.0	494	5	-0.017764	0.862055	FAS
3	60.0	80.0	70.0	371	5	-0.062985	0.913875	FAS
4	80.0	100.0	90.0	241	5	-0.035239	0.851155	FAS

{'morans_i': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/morans_i.csv'),
 'permutation_moran': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/permutation_moran.csv'),
 'distance_bin_correlations': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/distance_bin_correlations.csv')}

Run time: 572.3 ms

Clustering and PCA Spatial Modes

These summaries group stations with similar residual fingerprints and identify dominant station-level patterns. You will create a PCA summary figure for PGA and another one for FAS.

cluster_tables = []
cluster_score_tables = []
cluster_summary_tables = []
pca_score_tables = []
pca_loading_tables = []
pca_variance_tables = []

for metric_name, products in spatial_products.items():
    display(Markdown(f"### {spatial_metric_labels[metric_name]}"))

    # Build a station-by-event residual matrix for this metric.
    feature_table = build_station_feature_table(products["centered"])

    # Cluster stations with similar event-centered residual patterns for this metric.
    cluster_assignments, cluster_scores, cluster_features, cluster_summary, best_cluster, cluster_feature_columns = run_residual_feature_clustering(feature_table)

    # Decompose this metric's station residual patterns into PCA spatial modes.
    pca_result = compute_pca_spatial_modes(feature_table)

    # Plot the PC1 map, explained variance, and feature loading summary for this metric.
    pca_summary_fig = plot_pca_summary(
        pca_result.station_scores,
        pca_result.explained_variance,
        pca_result.feature_loadings,
        mode="PC1",
        title=f"{spatial_metric_labels[metric_name]} PCA Spatial Mode Summary",
        add_basemap=True,
        showfig=True,
        savefig=True,
        outpath=figure_dir / f"step_04_{metric_name.lower()}_pca_summary.png",
    )

    cluster_tables.append(cluster_assignments.assign(metric=metric_name))
    cluster_score_tables.append(cluster_scores.assign(metric=metric_name))
    cluster_summary_tables.append(cluster_summary.assign(metric=metric_name))
    pca_score_tables.append(pca_result.station_scores.assign(metric=metric_name))
    pca_loading_tables.append(pca_result.feature_loadings.assign(metric=metric_name))
    pca_variance_tables.append(pca_result.explained_variance.assign(metric=metric_name))
    display(pca_result.explained_variance)

clusters = pd.concat(cluster_tables, ignore_index=True)
cluster_scores = pd.concat(cluster_score_tables, ignore_index=True)
cluster_summary = pd.concat(cluster_summary_tables, ignore_index=True)
pca_station_scores = pd.concat(pca_score_tables, ignore_index=True)
pca_feature_loadings = pd.concat(pca_loading_tables, ignore_index=True)
pca_explained_variance = pd.concat(pca_variance_tables, ignore_index=True)

# Save clustering and PCA tables for the mapping notebook.
write_output_tables(
    clusters=clusters,
    cluster_scores=cluster_scores,
    cluster_summary=cluster_summary,
    pca_station_scores=pca_station_scores,
    pca_feature_loadings=pca_feature_loadings,
    pca_explained_variance=pca_explained_variance,
)

Peak acceleration (PGA)

	mode	mode_index	explained_variance	explained_variance_ratio	cumulative_explained_variance_ratio	singular_value	n_stations	n_features	standardized
0	PC1	1	2.592492	0.501215	0.501215	8.670771	30	5	True
1	PC2	2	0.986044	0.190635	0.691850	5.347456	30	5	True

Fourier amplitude spectrum (FAS), 3 s

	mode	mode_index	explained_variance	explained_variance_ratio	cumulative_explained_variance_ratio	singular_value	n_stations	n_features	standardized
0	PC1	1	1.567109	0.302974	0.302974	6.741376	30	5	True
1	PC2	2	1.209353	0.233808	0.536783	5.922098	30	5	True

{'clusters': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/clusters.parquet'),
 'cluster_scores': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/cluster_scores.csv'),
 'cluster_summary': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/cluster_summary.csv'),
 'pca_station_scores': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/pca_station_scores.parquet'),
 'pca_feature_loadings': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/pca_feature_loadings.csv'),
 'pca_explained_variance': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/pca_explained_variance.csv')}

Run time: 3.56 s

Geology Contrasts

Compare PGA and FAS residuals across the configured geology classes. GeoJSON regions and corridors are handled in the next notebook.

geology_tables = []

for metric_name, products in spatial_products.items():
    display(Markdown(f"### {spatial_metric_labels[metric_name]}"))

    # Compare configured geology classes with an event-bootstrap uncertainty estimate for this metric.
    geology_contrast = bootstrap_contrast_table(products["centered"], station_metadata=site_metadata).assign(metric=metric_name)

    # Plot the residual distributions and annotate the configured geology-class contrast for this metric.
    geology_contrast_fig = plot_geology_contrast(
        products["centered"],
        station_metadata=site_metadata,
        contrast_df=geology_contrast,
        title=f"{spatial_metric_labels[metric_name]} Residuals by Geology Class",
        showfig=True,
        savefig=True,
        outpath=figure_dir / f"step_04_{metric_name.lower()}_geology_contrast.png",
    )

    geology_tables.append(geology_contrast)
    display(geology_contrast)

geology_contrasts = pd.concat(geology_tables, ignore_index=True)

# Save the geology contrast output for later review.
write_output_tables(geology_contrasts=geology_contrasts)

Peak acceleration (PGA)

	group_col	left_values	right_values	comparison_values	baseline_values	contrast_label	effect_direction	value_col	statistic	min_stations_per_group	...	bootstrap_p	significant_95	significant_p05	n_left_stations	n_right_stations	n_comparison_stations	n_baseline_stations	n_events	bootstrap_samples_used	metric
0	mapped_region_type	Basin	Mountains	Basin	Mountains	Basin minus Mountains	mean(Basin) - mean(Mountains)	field_centered	mean	1	...	0.0	True	True	1	2	1	2	5	100	PGA

1 rows × 26 columns

Fourier amplitude spectrum (FAS), 3 s

	group_col	left_values	right_values	comparison_values	baseline_values	contrast_label	effect_direction	value_col	statistic	min_stations_per_group	...	bootstrap_p	significant_95	significant_p05	n_left_stations	n_right_stations	n_comparison_stations	n_baseline_stations	n_events	bootstrap_samples_used	metric
0	mapped_region_type	Basin	Mountains	Basin	Mountains	Basin minus Mountains	mean(Basin) - mean(Mountains)	field_centered	mean	1	...	0.72	False	False	1	2	1	2	5	100	FAS

1 rows × 26 columns

{'geology_contrasts': PosixPath('/home/runner/work/spatial-vtk/spatial-vtk/outputs/tutorials/tables/geology_contrasts.csv')}

Run time: 965.8 ms