Grouping Polygons by Attributes and by Distance Using GeoPandas and networkx

Question

I am working on creating a unique group ID for a parcel dataset I have containing 115k records. I have been attempting to do this with PostGIS (Grouping Parcels with an ID based on conditions, Query Speed assistance), but have not been able to get it working correctly using the cluster functions available.

Since I don't have much experience with PostGIS and SQL, I decided to stop my efforts and attempt doing this with geopandas and networkx.

My requirements: I am trying to group parcels and create a unique ID to represent these groups. Parcels should be grouped if they have the same parcel_owner value OR if they have the same parcel_owner_address value AND the parcels are within 90 meters of each other.

Here is my code.

import pandas as pd
import networkx as nx
import geopandas as gpd
from shapely.ops import unary_union
from shapely.geometry import Polygon

gdf = gpd.read_file('parcels_testing.gpkg')

# create a graph and add nodes
G = nx.Graph()
G.add_nodes_from(gdf.index)

# add edges based on name and address
for name, data in gdf.groupby('parcel_owner'):
    G.add_edges_from([(data.index[i], data.index[j]) for i in range(len(data)) for j in range(i + 1, len(data))])
for address, data in gdf.groupby('parcel_owner_address'):
    G.add_edges_from([(data.index[i], data.index[j]) for i in range(len(data)) for j in range(i + 1, len(data))])

components = list(nx.connected_components(G))

# create a new graph to handle distance-based grouping
G_distance = nx.Graph()
G_distance.add_nodes_from(gdf.index)

# add edges based on distance within each component
for component in components:
    sub_gdf = gdf.loc[list(component)]
    for i, row1 in sub_gdf.iterrows():
        for j, row2 in sub_gdf.iterrows():
            if i != j and row1['geometry'].distance(row2['geometry']) <= 90:
                G_distance.add_edge(i, j)

# find connected components based on distance
group_mapping = {}
for group_id, component in enumerate(nx.connected_components(G_distance)):
    for node in component:
        group_mapping[node] = group_id

# map the group id back to the GeoDataFrame
gdf['group_id'] = gdf.index.map(group_mapping)

It actually works really well from what I can tell. However, while testing data it appears to take about 1 minute to run a file containing 20k records, but for a file containing 50k records it is taking more than 4 minutes to run. Obviously, this does not scale for my purposes. I was hoping to eventually do this analysis on hundreds of thousands of records. Since I am iterating over the dataframe several times here (to create the edges based on name and address, and then again when I am adding edges based on distance), I imagine there is a much faster way to do this while being more organized.

Is there any way to do the above in a vectorized manner? I know the speed lag is most likely because I am iterating over all of the parcels, but I just can't think of a way around iteration in this scenario, since I need to check the parcels to see if they are within the distance parameter. Any suggestions are greatly appreciated.

EDIT: Updated Code

I was able to speed it up significantly making the following changes, which resulted in a 50k record analysis going from 4 minutes to 1 minute and 30 seconds. However, I am still doing a lot of things inefficiently I fear, and 115k record analysis is taking around 6 minutes. Any further suggestions would be greatly apricated.

Current version:

import itertools
import numpy as np
import pandas as pd
import networkx as nx
import geopandas as gpd
from shapely.ops import unary_union
from shapely.geometry import Polygon

# Load data
gdf = gpd.read_file('parcels.gpkg')
gdf = gdf[['combined_owner_name', 'combined_owner_address', 'geometry']]

# create graph
G = nx.Graph()

# Add nodes
G.add_nodes_from(gdf.index)

# Group the data and create edges
for col in ['parcel_owner', 'parcel_owner_address']:
    grouped = gdf.groupby(col).indices
    G.add_edges_from([edge for group in grouped.values() for edge in itertools.combinations(group, 2)])

# find connected components
components = list(nx.connected_components(G))

# Create graph
G_distance = nx.Graph()

# Add nodes
G_distance.add_nodes_from(gdf.index)

# Pre-calculate the buffer geometries
gdf['buffered_geometry'] = gdf['geometry'].buffer(90)

# Add edges based on distance within each component
for component in components:
    sub_gdf = gdf.loc[list(component)]
    spatial_index = sub_gdf.sindex

    for i, row in sub_gdf.iterrows():
        buffered_geometry = row['buffered_geometry']

        possible_matches_index = list(spatial_index.intersection(buffered_geometry.bounds))
        possible_matches = sub_gdf.iloc[possible_matches_index]
        precise_matches = possible_matches[possible_matches.intersects(buffered_geometry)]
        G_distance.add_edges_from([(i, j) for j in precise_matches.index if i != j])

# Find connected components based on distance
group_mapping = {}
for group_id, component in enumerate(nx.connected_components(G_distance)):
    for node in component:
        group_mapping[node] = group_id

# Map the group id back to the GeoDataFrame
gdf['group_id'] = gdf.index.map(group_mapping)
gdf = gdf[['group_id', 'combined_owner_name', 'combined_owner_address', 'geometry']]

Answer 1

I don't have geopandas, so I'll try to get close to it using ArcGis. This script

import pandas as pd
import itertools as itt

## create dataframe with geometries and names
aList =[]
with arcpy.da.SearchCursor("PARCELS",("OID@","Name","Shape@X","Shape@Y")) as cursor:
    for row in cursor:aList.append(row)
DF = pd.DataFrame(aList)
aList=[]
DF.columns = "FID","Name","X","Y"
DF.set_index("FID",inplace=True)

import networkx as nx
G = nx.Graph()
## iterate potential groups (repeat for address thing from here to end)
for name in DF.Name.unique():
    df = DF[DF.Name==name]
    for i,j in itt.combinations(df.index, 2):
        distance = math.hypot(df.X.loc[i]-df.X.loc[j],df.Y.loc[i]-df.Y.loc[j])
        if distance>=90.0:continue
        aList.append((i,j))
## add edges
G.add_edges_from(aList)

prints number of connected component equal to 21 for points shown:

Way more efficient solution can be done in ArcGIS using near table with limited search radius.