Building Adjacency Graphs from GeoDataFrames for GNNs

Convert a GeoDataFrame to a torch-geometric graph: build Queen/Rook contiguity, KNN, or distance-band adjacency and pack a COO edge_index with node features.

TL;DR Answer

Build a spatial weights object from your GeoDataFrame with libpysal.weights — Queen, Rook, KNN, or DistanceBand — convert its symmetric adjacency matrix to a COO edge_index tensor, attach a standardised node-feature matrix x, and wrap both in a torch_geometric.data.Data object. The one subtlety that breaks training is islands: units with no neighbours must be repaired before you hand the graph to a GNN. This is a foundational step in Graph Neural Networks for Spatial Data, where the graph topology is as much a modelling choice as the network architecture.

import numpy as np
import torch
from libpysal.weights import Queen
from torch_geometric.data import Data

w = Queen.from_dataframe(gdf, use_index=True)
w.transform = "b"                      # binary (0/1) weights
adj = w.sparse.tocoo()                 # symmetric COO adjacency
edge_index = torch.tensor(
    np.vstack([adj.row, adj.col]), dtype=torch.long
)
x = torch.tensor(gdf[feature_cols].to_numpy(), dtype=torch.float)
data = Data(x=x, edge_index=edge_index)

Part of: Graph Neural Networks for Spatial Data

Why the Graph Is a Modelling Decision, Not a Formality

A graph neural network learns by passing messages along edges. If two spatial units are connected, the model treats their features as mutually informative; if they are not, it treats them as independent. That means the adjacency graph is your spatial prior. Choosing Queen contiguity over KNN, or a 5 km distance band over a 10 km one, changes which relationships the network can exploit — often more than swapping the GNN layer type.

The mistake that silently degrades results is treating graph construction as boilerplate. Engineers copy a KNN(k=8) call from a tutorial, never inspect the resulting topology, and then wonder why their node-classification accuracy plateaus. The graph encodes the same neighbourhood structure that drives spatial lag and neighborhood statistics in classical spatial econometrics, and the same spatial autocorrelation that biases ordinary models. If the graph misrepresents which units interact, message passing propagates the wrong signal and no amount of hyperparameter tuning recovers it.

Three failure modes recur in production:

  • Silent islands. A DistanceBand graph over sparsely sampled rural units leaves some nodes with zero edges. Those nodes receive no messages, and layers such as GCNConv can emit NaN for them, poisoning the loss.
  • Asymmetric KNN. K-nearest-neighbour relationships are not reciprocal — A can be among B’s neighbours without B being among A’s. A directed edge_index breaks the symmetry most spatial GNNs assume.
  • Index drift. After a dropna() or a spatial filter, the GeoDataFrame index no longer matches row positions. Weights built on the old index produce edge_index entries that point at the wrong nodes.

Core Principles for a Sound Spatial Graph

  • Match adjacency to geometry. Polygons that tile space cleanly favour Queen or Rook contiguity; points and irregularly sized polygons favour KNN or a distance band.
  • Reset the index before building weights. Call gdf = gdf.reset_index(drop=True) so node IDs 0…N-1 align exactly with row positions and with the x matrix.
  • Force symmetry. A spatial GNN expects an undirected graph. Symmetrise KNN and distance graphs so every edge (i, j) has its mirror (j, i) in edge_index.
  • Repair islands explicitly. Never ship a graph with isolated nodes. Add the nearest neighbour or a self-loop, and log how many nodes were repaired.
  • Standardise node features first. GNNs are sensitive to feature scale exactly as neural nets are; apply the same feature scaling you would use for any deep model before packing features into x.
  • Keep coordinates for validation, not as the only feature. Store easting/northing so you can visualise and verify the graph, but let the network learn from substantive attributes, not raw position.

Production-Ready Code

The function below converts a GeoDataFrame into a torch_geometric.data.Data object. It validates inputs, supports contiguity, KNN, and distance-band methods, symmetrises the adjacency, repairs islands, and optionally adds self-loops. It returns a graph guaranteed to have no isolated nodes and a node count that matches the input.

import logging
from typing import Literal

import numpy as np
import torch
from libpysal.weights import Queen, Rook, KNN, DistanceBand, W
from torch_geometric.data import Data
from torch_geometric.utils import to_undirected, add_self_loops, contains_isolated_nodes

logger = logging.getLogger(__name__)

AdjMethod = Literal["queen", "rook", "knn", "distanceband"]


def gdf_to_pyg_data(
    gdf,
    feature_cols: list[str],
    method: AdjMethod = "queen",
    k: int = 8,
    threshold: float | None = None,
    y_col: str | None = None,
    add_loops: bool = False,
) -> Data:
    """Convert a GeoDataFrame into a torch-geometric Data graph.

    Parameters
    ----------
    gdf : geopandas.GeoDataFrame
        Spatial units (polygons or points) in a projected metric CRS.
    feature_cols : list[str]
        Columns packed into the node-feature matrix ``x``.
    method : {"queen", "rook", "knn", "distanceband"}
        Adjacency construction rule.
    k : int, default=8
        Number of neighbours when ``method="knn"``.
    threshold : float or None
        Distance band radius (CRS units) when ``method="distanceband"``.
    y_col : str or None
        Optional target column packed into ``data.y``.
    add_loops : bool, default=False
        If True, add a self-loop to every node after island repair.

    Returns
    -------
    torch_geometric.data.Data
        Graph with symmetric ``edge_index`` and no isolated nodes.
    """
    if gdf.crs is None or gdf.crs.is_geographic:
        raise ValueError(
            "gdf must be in a projected metric CRS (e.g. UTM) so KNN and "
            "distance thresholds are measured in metres, not degrees."
        )
    missing = [c for c in feature_cols if c not in gdf.columns]
    if missing:
        raise KeyError(f"feature_cols not found in gdf: {missing}")

    # Node IDs 0..N-1 must align with row positions and the x matrix.
    gdf = gdf.reset_index(drop=True)
    n_nodes = len(gdf)

    if method == "queen":
        w = Queen.from_dataframe(gdf, use_index=False)
    elif method == "rook":
        w = Rook.from_dataframe(gdf, use_index=False)
    elif method == "knn":
        w = KNN.from_dataframe(gdf, k=min(k, n_nodes - 1))
    elif method == "distanceband":
        if threshold is None:
            raise ValueError("distanceband requires a threshold in CRS units.")
        w = DistanceBand.from_dataframe(gdf, threshold=threshold, silence_warnings=True)
    else:
        raise ValueError(f"Unknown method: {method}")

    w.transform = "b"  # binary weights; edge presence, not row-standardised

    # Repair islands before building the tensor.
    if w.islands:
        logger.warning("%d island(s) detected; repairing with nearest neighbour.",
                       len(w.islands))
        w = _repair_islands(gdf, w)

    adj = w.sparse.tocoo()
    edge_index = torch.tensor(np.vstack([adj.row, adj.col]), dtype=torch.long)

    # Force an undirected graph (KNN is not reciprocal by construction).
    edge_index = to_undirected(edge_index, num_nodes=n_nodes)

    if add_loops:
        edge_index, _ = add_self_loops(edge_index, num_nodes=n_nodes)

    x = torch.tensor(gdf[feature_cols].to_numpy(dtype="float32"), dtype=torch.float)
    data = Data(x=x, edge_index=edge_index, num_nodes=n_nodes)

    if y_col is not None:
        data.y = torch.tensor(gdf[y_col].to_numpy(), dtype=torch.long)

    assert not contains_isolated_nodes(edge_index, num_nodes=n_nodes), \
        "Graph still contains isolated nodes after repair."
    return data


def _repair_islands(gdf, w: W) -> W:
    """Add each island's single nearest neighbour to the weights object."""
    from libpysal.weights import KNN as _KNN

    knn1 = _KNN.from_dataframe(gdf, k=1)
    neighbors = {i: list(w.neighbors[i]) for i in w.neighbors}
    for island in w.islands:
        nn = knn1.neighbors[island][0]
        neighbors[island].append(nn)
        neighbors[nn].append(island)  # keep it symmetric
    return W(neighbors, silence_warnings=True)

Setting w.transform = "b" keeps weights binary so edge_index records edge presence; if you need edge weights (for example inverse-distance), read them from w.sparse.data and pass them as edge_attr instead of discarding them.

Step-by-Step Walkthrough

1. Load regions and project to a metric CRS

import geopandas as gpd

# Administrative regions with attached attributes and a target class
gdf = gpd.read_file("regions.gpkg").to_crs("EPSG:32633")  # UTM 33N
feature_cols = ["pop_density", "mean_ndvi", "elevation", "road_length"]

print(f"Regions: {len(gdf)}")
print(f"CRS: {gdf.crs}  (projected metric — safe for KNN/distance)")

2. Standardise the node features

from sklearn.preprocessing import StandardScaler

gdf[feature_cols] = StandardScaler().fit_transform(gdf[feature_cols])

3. Build the graph with Queen contiguity

data = gdf_to_pyg_data(
    gdf,
    feature_cols=feature_cols,
    method="queen",
    y_col="land_class",
)
print(data)
# Data(x=[184, 4], edge_index=[2, 1012], y=[184], num_nodes=184)

4. Compare adjacency methods on the same units

Different rules produce different edge counts and mean degrees. Inspecting them before training tells you how dense a neighbourhood the GNN will see.

for method, kwargs in [
    ("queen", {}),
    ("rook", {}),
    ("knn", {"k": 6}),
    ("distanceband", {"threshold": 15000.0}),
]:
    g = gdf_to_pyg_data(gdf, feature_cols, method=method, **kwargs)
    mean_degree = g.edge_index.shape[1] / g.num_nodes
    print(f"{method:13s} edges={g.edge_index.shape[1]:5d}  mean_degree={mean_degree:.2f}")
queen         edges= 1012  mean_degree=5.50
rook          edges=  876  mean_degree=4.76
knn           edges= 1104  mean_degree=6.00
distanceband  edges= 1338  mean_degree=7.27

Rook contiguity is sparsest because it requires a shared edge, not just a shared vertex; KNN sits between the two by construction; the distance band is densest here because 15 km captures several rings of small regions.

5. Hand the graph to a GNN

import torch.nn.functional as F
from torch_geometric.nn import GCNConv


class SpatialGCN(torch.nn.Module):
    def __init__(self, in_dim: int, hidden: int, n_classes: int) -> None:
        super().__init__()
        self.conv1 = GCNConv(in_dim, hidden)
        self.conv2 = GCNConv(hidden, n_classes)

    def forward(self, x, edge_index):
        x = F.relu(self.conv1(x, edge_index))
        x = F.dropout(x, p=0.5, training=self.training)
        return self.conv2(x, edge_index)


model = SpatialGCN(in_dim=data.num_node_features, hidden=32, n_classes=int(data.y.max()) + 1)
out = model(data.x, data.edge_index)
print(out.shape)  # (184, n_classes)

How the Conversion Flows

The diagram below traces a GeoDataFrame through weights construction, island repair, symmetrisation, and tensor packing into the final Data object.

From GeoDataFrame to torch_geometric.data.Data A left-to-right pipeline: a GeoDataFrame of regions becomes a libpysal weights object, which is checked for islands and repaired, symmetrised into an undirected adjacency, and split into an edge_index tensor plus a node-feature matrix x that combine into a Data object. GeoDataFrame polygons / points + attributes libpysal.weights Queen / Rook KNN / DistanceBand Repair detect islands add nearest edge Symmetrise (i,j) and (j,i) undirected edge_index COO tensor [2, E] x (node features) float [N, F] Data torch_geometric ready for GNN edge_index and x are packed together into a single Data object with no isolated nodes.

Verification

Never trust a graph you have not audited. After construction, assert three invariants: the node count matches the input, edge_index is symmetric (every edge has its mirror), and no isolated nodes remain.

import torch
from torch_geometric.utils import contains_isolated_nodes, is_undirected


def verify_graph(data: Data, expected_nodes: int) -> None:
    """Fail loudly if the graph violates a spatial-GNN invariant."""
    assert data.num_nodes == expected_nodes, (
        f"node count {data.num_nodes} != expected {expected_nodes}"
    )

    # Symmetry: an undirected edge_index equals its own reverse.
    assert is_undirected(data.edge_index, num_nodes=data.num_nodes), \
        "edge_index is directed — symmetrise before training."

    # No index points outside the valid node range.
    assert int(data.edge_index.max()) < data.num_nodes, "edge index out of range"

    # No node is cut off from message passing.
    assert not contains_isolated_nodes(data.edge_index, num_nodes=data.num_nodes), \
        "isolated node(s) present — repair islands first."

    mean_degree = data.edge_index.shape[1] / data.num_nodes
    print(f"OK: {data.num_nodes} nodes, {data.edge_index.shape[1]} directed "
          f"edge entries, mean degree {mean_degree:.2f}")


verify_graph(data, expected_nodes=len(gdf))

If is_undirected fails on a KNN graph, you skipped the to_undirected call; if contains_isolated_nodes fails, a DistanceBand threshold left units unconnected. Both are construction bugs, not model bugs, and catching them here saves hours of debugging NaN losses later. The same defensive mindset underpins sound handling of spatial autocorrelation elsewhere in the training pipeline.

FAQ

Should I use contiguity, KNN, or a distance band for my adjacency graph?

Use Queen or Rook contiguity when your units are polygons that tile space with no gaps, such as census tracts or administrative regions — a shared border defines a natural topology. Use KNN when your units are points, or when polygons vary wildly in size, because it guarantees every node has exactly k neighbours and produces no islands. Use a distance band when a fixed interaction radius is physically meaningful, for example a sensor range or a species dispersal distance, but be ready to repair isolated nodes in sparse areas. When in doubt, compare mean degree across all three as shown above before committing.

How do I convert a libpysal weights object into a torch-geometric edge_index?

Set w.transform = "b", take w.sparse.tocoo() to get the symmetric adjacency in coordinate format, and stack its row and col arrays into a 2 × E LongTensor. Because libpysal weights are symmetric, both directions (i, j) and (j, i) are already present, so message passing flows both ways. Assign the result to Data.edge_index and confirm no entry exceeds num_nodes - 1. For KNN graphs, run the tensor through torch_geometric.utils.to_undirected because nearest-neighbour relationships are not reciprocal.

What do I do about islands (nodes with no neighbours)?

Islands break message passing: an isolated node receives no neighbourhood information and some layers emit NaN for it. Detect them with w.islands, then repair by adding each island’s nearest neighbour (keeping the edge symmetric) or, as a fallback, by adding self-loops so the node retains its own features. Always assert not contains_isolated_nodes(...) before building the Data object, and log how many nodes you repaired so the change is visible in your training records.

Part of: Graph Neural Networks for Spatial Data — Training Geospatial Predictive Models in Python