0

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']]

6
  • I'd use following workflow. Create combo column name_addindex. Group by it. Iterate every group, itertools combination by 2. If distance in pair less than 90, add pair to list. This list IS you graph, populate it by add_edges_from
    – FelixIP
    Jun 8 at 21:13
  • Note some nodes, might be islands, take it into account during component no transfer
    – FelixIP
    Jun 8 at 21:16
  • Hey @FelixIP not entirely sure what you mean. Maybe add some code to explain?
    – MapThug
    Jun 8 at 21:47
  • Vectorization of distance is only possible for points using scipy pdist method. Unfortunately count of points is limited by available memory
    – FelixIP
    Jun 9 at 4:36
  • Can you add a screenshot showing some parcels? Or even better share some test data
    – BERA
    Jun 9 at 17:02

1 Answer 1

0

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:

enter image description here

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

2
  • Hey Felix, thanks for taking the time to write this up. I am trying to stick to open source tools for this solution which is why I was using geopandas. If I get a license in the future I'll give this a go.
    – MapThug
    Jun 9 at 15:09
  • The only difference is distance calculation. I was hoping you understand idea: no need in 2 graphs, no need to add nodes to graph.
    – FelixIP
    Jun 9 at 20:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.