[Bug report] `ST_Intersects` not returning all expected entries
bug
### 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.
<img width="640" height="352" alt="Image" src="https://github.com/user-attachments/assets/c6f79e4c-416f-4c48-ae33-bc19db0165a9" />
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
```Python
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
```Python traceback
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
6 条评论