Main question
What is an efficient (i.e., good balance of memory & processing time) way to compute adjacency using a GeoPandas GeoDataFrame?
For my the rest of my question, I'll focus on the case where the GeoDataFrame contains Polygon shapes, but my question is fundamentally agnostic towards the type of geometries contained in the GeoDataFrame.
Reproducible example
# Loading important libraries
import geopandas as gpd
from shapely.geometry import Polygon
# Creating a function that makes square polygons
def make_square_polygon(centroid_x, centroid_y, side_length):
return (Polygon(
[(centroid_x - side_length/2, centroid_y - side_length/2),
(centroid_x + side_length/2, centroid_y - side_length/2),
(centroid_x + side_length/2, centroid_y + side_length/2),
(centroid_x - side_length/2, centroid_y + side_length/2)]))
# Creating the GeoDataFrame
gdf = gpd.GeoDataFrame({'id':list(range(1,10)),
'geometry':[make_square_polygon(0.5, 0.5, 1),
make_square_polygon(1.5, 0.5, 1),
make_square_polygon(2.5, 0.5, 1),
make_square_polygon(0.5, 1.5, 1),
make_square_polygon(1.5, 1.5, 1),
make_square_polygon(2.5, 1.5, 1),
make_square_polygon(0.5, 2.5, 1),
make_square_polygon(1.5, 2.5, 1),
make_square_polygon(2.5, 2.5, 1),]},
geometry='geometry')
# Plotting the GeoDataFrame
gdf.plot(color='white', edgecolor='black')
The code above produces the following plot (I've manually added the IDs in red):
Ultimately, my goal is to find out what polygons are neighbors (i.e., which polygons are touching each other). So for the example above, the results would represent the following:
- for the square with ID 1, the IDs of the squares that touch it are: 2, 4, and 5
- for the square with ID 2, the IDs of the squares that touch it are: 1, 3, 4, 5, and 6
- for the square with ID 3, the IDs of the squares that touch it are: 2, 5, and 6
- for the square with ID 4, the IDs of the squares that touch it are: 1, 2, 5, 7, and 8
- for the square with ID 5, the IDs of the squares that touch it are: 1, 2, 3, 4, 6, 7, 8 and 9
- for the square with ID 6, the IDs of the squares that touch it are: 2, 3, 5, 8 and 9
- for the square with ID 7, the IDs of the squares that touch it are: 4, 5, and 8
- for the square with ID 8, the IDs of the squares that touch it are: 4, 5, 6, 7, and 9
- for the square with ID 9, the IDs of the squares that touch it are: 5, 6, and 8
What the results should look like
The result could come out in multiple different ways.
Optimal results
Optimally, the results would look like this:
- an array with 9 positions, each one representing the rows in the original GeoDataFrame. This array is populated with lists of the IDs or indices of the adjacent geometries (similar to what I have described in the example above).
- a long Mx2 matrix with each row containing a pair of IDs or indices from the original GeoDataFrame that touch each other. In this case, M represents the number of pairs of features that touch each other.
Note that the list in the example above has several repeats (e.g., 1 & 2 and 2 & 1). Ideally, this list might even be further simplified to remove those duplicates.
2nd best results
Less optimally, the results would look like this:
- a 9x9 matrix with
Trues andFalses
- a long (R²)x3 matrix representing all possible combinations of pairs of elements (eg: 1-1, 1-2, 1-3, 2-1, 2-2, 2-3, ...). Here, the first and second column would represent the IDs or indices of the elements of the original GeoDataFrame and the third column would be filled with boolean values representing whether or not the features intersect. In this case, R represents the number of elements in the original GeoDataFrame.
I classify these as "less optimal" because they require a lot of memory to store (especially when dealing with large datasets).
My attempt
Here's my attempt at writing a function that performs the operation described above:
def get_adj_matrix(input_gdf):
nrows = input_gdf.shape[0]
geometry_column_name = input_gdf._geometry_column_name
temp = input_gdf.merge(input_gdf, how='cross')
adj_matrix = (gpd.GeoSeries(temp[f'{geometry_column_name}_x'])
.touches(gpd.GeoSeries(temp[f'{geometry_column_name}_y']))
.values.reshape((nrows, nrows)))
return adj_matrix
print(get_adj_matrix(gdf))
# array([[False, True, False, True, True, False, False, False, False],
# [ True, False, True, True, True, True, False, False, False],
# [False, True, False, False, True, True, False, False, False],
# [ True, True, False, False, True, False, True, True, False],
# [ True, True, True, True, False, True, True, True, True],
# [False, True, True, False, True, False, False, True, True],
# [False, False, False, True, True, False, False, True, False],
# [False, False, False, True, True, True, True, False, True],
# [False, False, False, False, True, True, False, True, False]])
The code above works just fine. The problem, however, is two-fold:
- The code is messy/ugly. It's weird that I need to do a
crossmerge, create a vector that has (R²) elements and then reshape it back to be an RxR matrix. - The function is very fast because it uses NumPy's/Pandas'/GeoPandas' vectorized methods. However, for larger datasets, it is going to consume A LOT of memory because this approach requires us to make a dataframe of size (R²)x(2*C) and then a matrix of size RxR. (R represents the number of elements/rows in the original GeoDataFrame and C represents the number of columns in the original GeoDataFrame).
Simply put, my approach generates a 2nd best answer (as defined above) and I'd like to compute an Optimal answer which not only looks optimal in the end, but also avoids going through the large memory-consuming intermediate steps where huge RxR matrices or (R²)x(2*C) DataFrames are generated.
Back to my main question
How can I quickly compute adjacency among elements within a GeoDataFrame and have the results be computed in a way that doesn't take up much unnecessary memory (as exposed in the "Optimal results" section above) and doesn't depend on Python's rather slow for-loops? If possible, it'd be great to avoid the use of apply methods as well.
