# What is the algorithm that Shapely used to check if two polygons intersect?

What is the algorithm that Shapely used to check if two polygons intersect?

``````from shapely.geometry import Polygon
p1 = Polygon([(0,0), (3,0), (3,1), (1,1), (1,2), (3,2), (3,3), (0,3)])
p2 = Polygon([(4,0), (5,0), (5,1.5), (2,1.5), (2,1.2), (4,1.2)])
print(p1.intersects(p2))
``````

I had read the source code of Shapely but didn't find the implementation.

• Well, it's open-source, so you can track it down, but according to this, your search will take you to JTS. Commented Nov 14, 2020 at 15:28
• Thank you! I have read the source code of the Shapely, but I didn't find the implementation. Thank again for the JTS. Commented Nov 14, 2020 at 17:32
• Please be sure to include all relevant facts in the Question body when you post. The Tour has additional details about how GIS SE works. Commented Nov 14, 2020 at 17:53
• Thanks! I have modified the question. Commented Nov 14, 2020 at 17:58

A quick search of the Shapely code base leads to `impl.py` which:

"""Implementation of the intermediary layer between Shapely and GEOS This is layer number 2 from the list below.

1. geometric objects: the Python OO API.
2. implementation map: an abstraction that permits different backends.
3. backend: callable objects that take Shapely geometric objects as arguments and, with GEOS as a backend, translate them to C data structures.
4. GEOS library: algorithms implemented in C++. Shapely 1.2 includes a GEOS backend and it is the default. """

So the next place to look is GEOS and `Geometry.ccp` and it's `intersection` method. This is a wrapper for `HeuristicOverlay` which calls OverlayNGRobust::overlay which is where we start to see references to JTS, where the API docs state:

The obvious naive algorithm for intersection detection (comparing every segment with every other) has unacceptably slow performance. There is a large literature of faster algorithms for intersection detection. Unfortunately, many of them involve substantial code complexity. JTS tries to balance code simplicity with performance gains. It uses some simple techniques to produce substantial performance gains for common types of input data.

After some more poking around in the JTS code we eventually come to `RelateComputer` which calculates the `IntersectionMatrix` between two geometries (which includes intersection).

Any deeper than that and you will need to wait for Dr JTS to answer.

• And answer I will, but not with anything better than "read the source code" ATM. The spatial relationship algorithm is quite complex, and currently there's no good high-level description of it available. Maybe sometime... Commented Dec 10, 2020 at 17:13
• I will point out that the OverlayNG code is only used for overlay ops (like `intersection`). The spatial predicates (like `intersects`) are evaluated using an older and completely different codebase (which `RelateComputer` is an entry point for). Commented Dec 10, 2020 at 17:14
• For future readers who are interested in how JTS/Shapely implements polygon intersection, see my answer here. Commented Sep 16, 2021 at 17:29