Efficient Spatial Joins on Large Geometry Sets

We consider the following standard spatial-join problem: Given two sets of geometric objects in 2D (points, lines, polygons, and collections of these), compute the spatial relations of all pairs of intersecting objects as a standard DE-9IM matrix. Most previous work focuses on one aspect of the problem, like the candidate generation, candidate reduction heuristics, efficient data structures, or parallelization. We provide a complete, fully functional, and carefully engineered implementation, as well as an extensive experimental evaluation of the relevance of various heuristics and of three variants for the exact geometry comparisons: our own implementation which preprocesses the geometries, one using the standard geometries from the GEOS library (which powers spatial joins in the widely used PostgreSQL+PostGIS), and one using the implementation based on preprocessed (so-called “prepared”) geometries offered by GEOS (the last variant is possible since GEOS 3.13). In particular, we find that compared to the standard GEOS implementation, our own implementation speeds up spatial joins by more than an order of magnitude when complex geometries are involved and no prefilter heuristics are used, and is still consistently faster than the “prepared” variant offered by GEOS . Our best approach can compute the full DE-9IM self-join of the 1.4 billion geometries from OpenStreetMap in 2 hours on a commodity PC. This was out of reach for any existing implementation we tried. We also extend our approach to work with within-distance( \(\hat{d} \) ) joins, which compute the distance of all object pairs within a user-defined distance \(\hat{d} \) , using both Euclidean and great-circle meter distances. Our code and all the materials needed to reproduce our results are freely available on GitHub.