[How to use] `ST_Intersects` not returning all expected entries

[pjonsson]

Describe the bug

The attached code is a reproducer with the latest GeoAlchemy2/SQLAlchemy 2.0/psycopg3 running against a postgis:18-3.6 Docker container. The issue was present in June with postgis:16-3.5 as well, but I ran out of time to look into this back then.

The reproducer includes 2 polygons, the small test set that the reproducer was extracted from contained 198 polygons, and the code without the cast returns 192 overlapping polygons, with the cast it returns all 198 polygons.

Here is a screenshot of the polygon that (supposedly) doesn't overlap, with the large polygon extending beyond the screenshot in north/east/west directions.

Unfortunately I inherited the code, so I don't have any background information on the MyGeography class, but the code did contain a comment that links to issue 371.

Optional link from https://geoalchemy-2.readthedocs.io which documents the behavior that is expected

No response

To Reproduce

import os
from collections.abc import Sequence
from typing import Any, override

import shapely.wkt
from geoalchemy2 import Geography, WKBElement, WKTElement, Geometry
from geoalchemy2.functions import ST_Intersects
from geoalchemy2.shape import from_shape, to_shape
from shapely import Polygon
from shapely.geometry.base import BaseGeometry
from sqlalchemy import TypeDecorator, create_engine, func, select, text, cast
from sqlalchemy.orm import (
    DeclarativeBase,
    Mapped,
    Session,
    mapped_column,
    registry,
    sessionmaker,
)

db_user = os.getenv("DB_USER")
db_password = os.getenv("DB_PASSWORD")
db_host = os.getenv("DB_HOST")

db_connection_base = f"postgresql+psycopg://{db_user}:{db_password}@{db_host}:5432/"


def create_session() -> Session:
    """Connect to postgres database, create the test database, and return a session"""
    engine = create_engine(
        db_connection_base + "postgres", isolation_level="AUTOCOMMIT"
    )
    with engine.connect() as conn:
        conn.execute(text("DROP DATABASE IF EXISTS repro_intersect"))
        conn.execute(text("CREATE DATABASE repro_intersect"))
        conn.commit()
    engine.dispose(close=True)
    engine = create_engine(
        db_connection_base + "repro_intersect", plugins=["geoalchemy2"], echo=True
    )
    with engine.connect() as conn:
        conn.execute(text("CREATE EXTENSION IF NOT EXISTS postgis"))
        conn.commit()
    return sessionmaker(bind=engine)()


class Base(DeclarativeBase):
    pass


mapper_registry = registry(metadata=Base.metadata)


class MyGeography(TypeDecorator):
    # Keep a private reference with proper type for type checking.
    _impl = Geography
    impl = _impl
    cache_ok = True

    # https://github.com/geoalchemy/geoalchemy2/issues/371
    @override
    def column_expression(self, column: Any) -> Any:
        st_asbinary = self._impl.as_binary
        assert st_asbinary is not None
        return getattr(func, st_asbinary)(column, type_=self)

    @override
    def process_bind_param(
        self, value: WKBElement | Any, dialect: Any
    ) -> WKBElement | None:
        if value is None:
            return value  # pragma: no cover
        return from_shape(value, srid=4326)

    @override
    def process_result_value(
        self, value: WKBElement | WKTElement | None, dialect: Any
    ) -> BaseGeometry | None:
        if value is None:
            return value  # pragma: no cover
        return to_shape(value)


class Thing(Base):
    __tablename__ = "things"

    id: Mapped[int] = mapped_column(primary_key=True, autoincrement=True)
    area: Mapped[Polygon] = mapped_column(
        MyGeography(geometry_type="GEOMETRY", srid=4326)
    )


session = create_session()
conn = session.connection()
mapper_registry.metadata.drop_all(conn, checkfirst=True)
mapper_registry.metadata.create_all(conn)

# Insert test data to database.
for p in [
    # Intersects, and is found.
    "POLYGON((34.63890642844616 56.83305528502311, "
    "34.59716193895972 55.84695640311729, "
    "36.34836128644486 55.81173603255935, "
    "36.43571634947202 56.79650515210518, "
    "34.63890642844616 56.83305528502311))",
    # Intersects, but not found without cast.
    "POLYGON((22.65868601431409 53.58550153864066,"
    "22.68033586832527 54.13637742873316, "
    "20.99969379936479 54.14810410396127, "
    "20.99970086834074 53.16117354480671, "
    "22.41446155119707 53.15142479978014, "
    "22.48062749772035 53.27022660537294, "
    "22.56125560535012 53.41365952061176, "
    "22.6424609927254 53.55701784114553, "
    "22.65868601431409 53.58550153864066))",
]:
    t = Thing(area=shapely.wkt.loads(p))
    session.add(t)
    session.flush()
    session.expire(t)


# Query parts.
def get_intersects(geometry: Polygon, use_cast: bool) -> Sequence[Thing]:
    column = cast(Thing.area, Geometry(srid=4326)) if use_cast else Thing.area
    result = select(Thing).filter(ST_Intersects(column, from_shape(geometry, srid=4326)))
    return session.execute(result).scalars().all()

things = get_intersects(
    shapely.wkt.loads("POLYGON((4.3 53.5, 4.3 71.5, 35 71.5, 35 53.5, 4.3 53.5))"), False
)

print(f"\n\nFound {len(things)}/2 things that overlap (expected 2/2)\n\n")

more_things = get_intersects(
    shapely.wkt.loads("POLYGON((4.3 53.5, 4.3 71.5, 35 71.5, 35 53.5, 4.3 53.5))"), True
)

print(f"\n\nFound {len(more_things)}/2 things that overlap (expected 2/2)\n\n")

Error

2025-10-22 12:04:11,098 INFO sqlalchemy.engine.Engine SELECT things.id, ST_AsBinary(things.area) AS area
FROM things
WHERE ST_Intersects(things.area, ST_GeomFromWKB(%(ST_GeomFromWKB_1)s, 4326))
2025-10-22 12:04:11,098 INFO sqlalchemy.engine.Engine [generated in 0.00007s] {'ST_GeomFromWKB_1': <memory at 0x795e0048fe80>}


Found 1/2 things that overlap (expected 2/2)


2025-10-22 12:04:11,102 INFO sqlalchemy.engine.Engine SELECT things.id, ST_AsBinary(things.area) AS area
FROM things
WHERE ST_Intersects(CAST(things.area AS geometry(GEOMETRY,4326)), ST_GeomFromWKB(%(ST_GeomFromWKB_1)s, 4326))
2025-10-22 12:04:11,102 INFO sqlalchemy.engine.Engine [generated in 0.00009s] {'ST_GeomFromWKB_1': <memory at 0x795e0048f940>}


Found 2/2 things that overlap (expected 2/2)

Additional context

No response

GeoAlchemy 2 Version in Use

0.18.0

Python Version

3.12

Operating system

Linux

[adrien-berchet]

Hi @pjonsson

Sorry for the very late answer, I'm currently on vacation and I won't be able to take much time on this for the next 2 weeks. I hope it's not a blocker for you...

Anyway, thank you very much for this very detailed bug report!

At first glance I think it's not related to GeoAlchemy2 since the queries reported by the engine look correct so you would probably have the same issue with raw queries. But I will check when I have some time.

[adrien-berchet]

Oh wait, I think it's just because you use a Geography type instead of Geometry. Using geography means that all the calculations use spherical coordinates. In your case you define a giant rectangle with only 4 points (the corners), so the edges are very distorted in spherical coordinates. You can add intermediate points on the edge if you want an almost-Cartesian rectangle. Or maybe it's actually normal that these rectangles don't intersect, it depends on your use case?

[pjonsson]

I do get the same behavior when running the generated SQL query in psql, so you're right that the same issue appears with raw queries.

Your suggestion that the square would be distorted and the results not intersect sounds feasible. I have a few other more urgent things on my plate right now, but I will look into that later this week.

And don't worry about this bug for now, you deserve your vacation.

[adrien-berchet]

Hi @pjonsson
Were you able to find some time to check this?

[pjonsson]

@adrien-berchet sorry, too many things piled up at once, I'll get back to this in late January.

Merry Christmas and happy new year!

[adrien-berchet]

Hi @pjonsson
No problem, take your time.
Merry Christmas to you and your family and see you next year! 😊

[pjonsson]

I'm trying to allocate some time from someone who knows more about this stuff than I do, but things are dragging out on time. Sorry for the long silence, the issue is still on my radar.