TopologyException: side location conflict while intersects on valid polygons

Question

I am cheking if two valid polygons intersects using GEOS. The error is the same when I use Djangos GEOS binding or PostGIS which make me think this is a GEOS related issue.

Here is the error when using Django binding:

GEOS_ERROR: TopologyException: side location conflict at 1.839012980156925 43.169860517728324. This can occur if the input geometry is invalid.

Traceback (most recent call last):
  File "/app/manage.py", line 21, in <module>
    main()
  File "/app/manage.py", line 17, in main
    execute_from_command_line(sys.argv)
  File "/venv/lib/python3.10/site-packages/django/core/management/__init__.py", line 442, in execute_from_command_line
    utility.execute()
  File "/venv/lib/python3.10/site-packages/django/core/management/__init__.py", line 436, in execute
    self.fetch_command(subcommand).run_from_argv(self.argv)
  File "/venv/lib/python3.10/site-packages/django/core/management/base.py", line 413, in run_from_argv
    self.execute(*args, **cmd_options)
  File "/venv/lib/python3.10/site-packages/django/core/management/base.py", line 459, in execute
    output = self.handle(*args, **options)
  File "/app/batid/management/commands/sandbox.py", line 41, in handle
    p1.intersects(p2)
  File "/venv/lib/python3.10/site-packages/django/contrib/gis/geos/geometry.py", line 334, in intersects
    return capi.geos_intersects(self.ptr, other.ptr)
  File "/venv/lib/python3.10/site-packages/django/contrib/gis/geos/libgeos.py", line 154, in __call__
    return self.func(*args)
  File "/venv/lib/python3.10/site-packages/django/contrib/gis/geos/prototypes/threadsafe.py", line 47, in __call__
    return self.cfunc(self.thread_context.handle.ptr, *args)
  File "/venv/lib/python3.10/site-packages/django/contrib/gis/geos/prototypes/errcheck.py", line 57, in check_predicate
    raise GEOSException(
django.contrib.gis.geos.error.GEOSException: Error encountered on GEOS C predicate function "GEOSIntersects_r".

Here is the Python code triggering the error and also some working examples below:

The problematic shapes

from django.contrib.gis.geos import GEOSGeometry

# Problematic case
# This triggers error
p1 = GEOSGeometry("POLYGON ((1.839012980156925 43.169860517728324, 1.838983490127865 43.169860200336274, 1.838898525601717 43.169868281549725, 1.838918565176068 43.1699719478626, 1.838920733577112 43.16998636433192, 1.838978629555589 43.16997979090823, 1.838982586839382 43.169966339940714, 1.838974943184281 43.169918580432174, 1.839020497362873 43.169914572864634, 1.839012980156925 43.169860517728324))")
p2 = GEOSGeometry("POLYGON ((1.8391355300979277 43.16987802887805, 1.83913336164737 43.16986361241434, 1.8390129801569248 43.169860517728324, 1.8390790978572837 43.16987292371998, 1.8390909520103162 43.16995581178317, 1.8391377530291442 43.16995091801345, 1.8391293863398452 43.16987796276235, 1.8391355300979277 43.16987802887805))")
p1.intersects(p2)

Here is what the shapes look like:

The reported point in the error message (1.839012980156925 43.169860517728324) is the tip of the shard

Both shapes are marked as valid.

Similar shape but no error

I thought GEOS had a specific problem on checking intersection of two polygon sharing one exact point at the end of the shard. I created those two shapes but the intersects works well

from django.contrib.gis.geos import GEOSGeometry

# Similar shapes, sharing one point at the end of a shard
# No errors
p1 = GEOSGeometry("POLYGON((1.8391709402203558 43.16998351248077,1.8390354886651037 43.169982045324815,1.8389851972460745 43.170070074617314,1.8391260132193563 43.17009892863558,1.8391803279519081 43.17008670235833,1.8391709402203558 43.16998351248077))")
p2 = GEOSGeometry("POLYGON((1.8391709402203558 43.16998351248077,1.8394726887345314 43.169986935844406,1.8394968286156654 43.17004806730611,1.839553825557232 43.17004122058543,1.8395008519291878 43.169959059877556,1.8391709402203558 43.16998351248077))")
p1.intersects(p2)

Trying with two square sharing one point

from django.contrib.gis.geos import GEOSGeometry

# Simple squares sharing one point
# no error on simple square with only one touching point
square1 = GEOSGeometry("POLYGON ((1 1, 1 2, 2 2, 2 1, 1 1))")
square2 = GEOSGeometry("POLYGON ((2 2, 2 3, 3 3, 3 2, 2 2))")
square1.intersects(square2)

Hacky workarounds

Some workarounds remove the error but I prefer to find a clean solution. Here is what I found

• Converting the geometry into GeoJSON reduces the decimals number. Converting this GeoJSON back to GEOS geometry remove the error but the intersects verifiation return False
• Same result using geom.buffer(0.000001).buffer(-0.000001)
• using buffer(0) does nothing

Answer 1

This is a (rare) instance of a robustness problem with the current topological predicate algorithm in JTS/GEOS. Luckily, there is a new implementation in development which solves all of these problems. See this blog post and this GEOS issue for some background.

Also, current versions of PostGIS and GEOS (3.12) have a fix in place for this problem as well, for the specific case of intersects. So upgrading may solve this as well.