[BUG] `contiguity_graph` returns undirected edge table

[stefdesabbata]

Contributing guidelines

• I understand the contributing guidelines

Documentation

• My proposal is not addressed by the documentation or examples

Existing issues

• Nothing similar appears in an existing issue

Describe the bug

The contiguity_graph(..., as_nx=False) method returns an edges_gdf in which each contiguity edge appears only once, with the smaller (in the order of the type) node id as the first MultiIndex level. The total edge count is correct for undirected edges. However, when transformed into a PyG graph, the edge table is assumed to be undirected, leading to missing edges.

PyG treats every row of edge_index as a directed edge. The edge_index is a 2 x |E| matrix where (under the default flow="source_to_target") the messages from the nodes in edge_index[0] are sent to nodes in edge_index[1]. With gdf_to_pyg, level 0 of the multi-index becomes the source and level 1 becomes the target node that receives the message. As each edge appears only once in the edge table, only one direction is included in the resulting edge_index.

To Reproduce

import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon
import city2graph as c2g

# Create simple geodataframe
polys, codes = [], []
for r in range(4):
    for c in range(4):
        polys.append(Polygon([(c, r), (c + 1, r), (c + 1, r + 1), (c, r + 1)]))
        codes.append(f"id{r}{c}")
gdf = (
    gpd.GeoDataFrame({"v": range(16)}, geometry=polys, crs="EPSG:27700")
    .set_index(pd.Index(codes, name="id"))
)

# Plot
ax = gdf.plot(color="lightgray", categorical=True, edgecolor="black")
for idx, row in gdf.iterrows():
    ax.annotate(idx, xy=(row.geometry.centroid.x, row.geometry.centroid.y),
                ha="center", va="center")
    ax.set_axis_off()

nodes, edges = c2g.contiguity_graph(gdf, contiguity="queen", as_nx=False)

pyg_graph = c2g.gdf_to_pyg(nodes, edges=edges)

# Check output, using id22 as an example
print("\n\nEdges with id22 as level 0 in the multi-index", flush=True)
print(edges.xs("id22", level=0), flush=True)
print("\n\nEdges with id22 as focal", flush=True)
print(pyg_graph.edge_index[:, pyg_graph.edge_index[0] == nodes.loc["id22"]["v"]], flush=True)

# Plot id22 and its neighbors where id22 is level 0 in the multi-index
id22_edges = edges.xs("id22", level=0)
ax = gdf.plot(color="lightgray", categorical=True, edgecolor="black")
nodes[nodes.index.isin(id22_edges.index)].plot(ax=ax, color="tab:orange")
nodes[nodes.index == "id22"].plot(ax=ax, color="tab:blue")
id22_edges.plot(ax=ax, color="tab:blue")
for idx, row in gdf.iterrows():
    ax.annotate(idx, xy=(row.geometry.centroid.x, row.geometry.centroid.y),
                ha="center", va="center")
    ax.set_axis_off()

# Plot id22 and its neighbors (both levels in the multi-index)
# this is what PyG expects
ax = gdf.plot(color="lightgray", categorical=True, edgecolor="black")
nodes[nodes.index.isin(id22_edges.index)].plot(ax=ax, color="tab:orange")
nodes[nodes.index == "id22"].plot(ax=ax, color="tab:blue")
id22_edges.plot(ax=ax, color="tab:blue")

# Add edges where id22 is on level 1
id22_lvl1 = edges.xs("id22", level=1)
nodes[nodes.index.isin(id22_lvl1.index)].plot(ax=ax, color="tab:red")
id22_lvl1.plot(ax=ax, color="tab:blue")

for idx, row in gdf.iterrows():
    ax.annotate(idx, xy=(row.geometry.centroid.x, row.geometry.centroid.y),
                ha="center", va="center")    
ax.set_axis_off()

Expected behavior

All directed edges should be added for queen adjacency.

Environment (please complete the following information):

• Python version: 3.12
• city2graph version: 0.3.1
• torch version: 2.8
• torch_geometric version: 2.7

How did you install city2graph?

via PyPI (e.g., pip)

Additional context

No response

[yu-ta-sato]

Thank you so much for raising this issue with reproduction! As discussed, gdf_to_pyg and pyg_to_gdf will need to allow duplication of edges when undirected graph is fed. Let me investigate the source codes and identify which part should be tackled, before raising an PR.

[stefdesabbata]

It might be worth checking out the approach taken by libpysal with libpysal.graph.Graph.to_networkx(), which if "Graph is symmetric, returns nx.Graph, otherwise returns a nx.DiGraph." and then how torch_geometric.utils.from_networkx() converts them.

[yu-ta-sato]

Thanks! As of now, contiguity_graph returns nx.Graph (undirected).

I replicated the issue and investigate the root cause, and yes it was both from the de-duplication of _generate_contiguity_edges and directed deal of gdf_to_pyg. Thanks for this report again!

Next action as PR

First and foremost, gdf_to_pyg's treatment of undirected graphs should be fixed, which is the prioritised one.
gdf_to_pyg should ensure the duplication of edge_index if undirected graphs are set as input, holding the original edge geometries in the metadata to be reconstructed by pyg_to_gdf. For this fix, the use of torch_geometric.utils.to_undirected would be the key that can automatically add reverse edges, without duplicating the geometries.

What about the exclusive separation in gdfs of undirected graphs?

In the latest update of v0.3.0 -> v0.3.1, I omitted the duplicated edges to reduce the burden of duplicated edge geometries. As of now, gdf_to_nx handles the indirected graph to nx.Graph, it would be better to keep this de-duplication to halve the volume of geometries if undirected edges are enough.

This would be a dilemma of two needs from potial users:

1. those who would like to reduce data volume because of duplications
2. those who would like to handle graph-based operation from gdfs (like filtering neighbours from MultiIndex)

Some solution is to set a parameter whether to allow duplicaiton or not in undirected graph scenarios.
Let me work on it after the PR for gdf_to_pyg mentioend above.

By the way, in a separated PR in the future, directed option could also be enhanced to contiguity_graph

[stefdesabbata]

I think getting the right graph for PyG could be done using PySAL's Graph.adjacency, rather than going necessarily through nx, which can be quite computationally expensive for large graphs. Could you double-check if the example makes sense based on PyG's MessagePassing.propagate?

import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon
from libpysal import graph
import numpy as np
import torch

# Create simple geodataframe
polys, codes = [], []
for r in range(4):
    for c in range(4):
        polys.append(Polygon([(c, r), (c + 1, r), (c + 1, r + 1), (c, r + 1)]))
        codes.append(f"id{r}{c}")
gdf = (
    gpd.GeoDataFrame({"v": range(16)}, geometry=polys, crs="EPSG:27700")
    .set_index(pd.Index(codes, name="id"))
)
print(gdf)
print("\n\n")

gdf = gdf.reset_index().set_index("v")
print(gdf)
print("\n\n")

g_queen   = graph.Graph.build_contiguity(gdf, rook=False)
adj_queen = g_queen.adjacency

focal_idx    = adj_queen.index.get_level_values(0).to_numpy()
neighbor_idx = adj_queen.index.get_level_values(1).to_numpy()

edge_index_np = np.stack([
        # source_to_target: 
        # messages go neighbour -> focal
        neighbor_idx, focal_idx
    ], axis=0)
edge_index = torch.from_numpy(edge_index_np).long()

print(edge_index)
print("\n\n")

## Neighbour id22 receives messages from
id22_mess_from = edge_index[:, edge_index[1] == gdf[gdf["id"] == "id22"].index[0]]
print(id22_mess_from)
print("\n\n")
print(gdf.loc[id22_mess_from[0]])
print("\n\n")

[yu-ta-sato]

@stefdesabbata Sorry it's a coincidence that I issued the PR at the same timing now!
And thanks a lot for the idea! As gdf_to_pyg and pyg_to_gdf deal with graphs of non-polygon geometries (especially points) in general, contiguity-based approach from Graph.adjacency may not work for the other functions.

In #157, edge_index was ensured to be duplicated if a new parameter directed is False in gdf_to_pyg (I'll enhance nx_to_pyg too where nx.(Multi)Graph or nx.(Multi)Digraph would be the indicators). Once the objects are converted, the flag of directed will be stored in the class attrbutes as metadata for reconstruction in pyg_to_gdf.

In this edge_index duplication, only indices (tourch.Tensor) are manipulated, while the geometries in WKT are held as intact. When reconstructed by pyg_to_gdf, the original geometries will be recovered.

I prioritised the correctness of PyG objects when users intend to convert gdf or nx (nx_to_gdf will be fixed in the same PR). In this case, neighbouring queries in gdfs would still have exclusive separation. In each graph construction functions including contiguity_graph, a new parameter would be added (i.e. duplicate_edges=False) for those who want to query neighbours from gdfs.