Hybrid Properties for Geometry in SQLAlchemy

-- 1. Confirm PostGIS is installed
SELECT extname, extversion
FROM pg_extension
WHERE extname = 'postgis';

-- 2. Confirm the geometry column exists with the expected type and SRID
SELECT f_table_name, f_geometry_column, type, srid, coord_dimension
FROM geometry_columns
WHERE f_table_name = 'parcels';

-- 3. Confirm a GiST index exists on the geometry column
-- (required for ST_DWithin and && bounding-box operators used alongside hybrids)
SELECT indexname, indexdef
FROM pg_indexes
WHERE tablename = 'parcels'
  AND indexdef ILIKE '%gist%';

CREATE INDEX CONCURRENTLY idx_parcels_boundary_gist
    ON parcels USING GIST (boundary);
ANALYZE parcels;

from sqlalchemy.orm import DeclarativeBase
from sqlalchemy import Column, Integer, String
from geoalchemy2 import Geometry

class Base(DeclarativeBase):
    pass

class Parcel(Base):
    __tablename__ = "parcels"

    id = Column(Integer, primary_key=True)
    name = Column(String(100), nullable=False)
    # Explicit SRID=4326 and geometry_type prevents silent coercion errors
    boundary = Column(
        Geometry(geometry_type="POLYGON", srid=4326),
        nullable=False,
    )

from sqlalchemy.ext.hybrid import hybrid_property
from sqlalchemy import func
from geoalchemy2.shape import to_shape
from shapely.ops import transform
import pyproj

# Build the transformer once at module load time; reuse across instances
_wgs84_to_equal_area = pyproj.Transformer.from_crs(
    "EPSG:4326", "EPSG:6933", always_xy=True
).transform

class Parcel(Base):
    __tablename__ = "parcels"

    id = Column(Integer, primary_key=True)
    name = Column(String(100), nullable=False)
    boundary = Column(
        Geometry(geometry_type="POLYGON", srid=4326),
        nullable=False,
    )

    @hybrid_property
    def area_sqkm(self) -> float:
        """Area in km² computed on a loaded instance using an equal-area CRS.

        Projects from WGS84 to EPSG:6933 (World Equal Area) before calling
        Shapely .area so the result is in square metres, then converts to km².
        """
        geom = to_shape(self.boundary)           # WKBElement → Shapely geometry
        projected = transform(_wgs84_to_equal_area, geom)
        return projected.area / 1_000_000        # m² → km²

    @area_sqkm.expression
    @classmethod
    def area_sqkm(cls):
        """SQL expression: cast to geography triggers spheroidal ST_Area in m².

        The division by 1 000 000 converts to km² inline, so ORDER BY and
        WHERE clauses see the same scale as the Python instance method.
        """
        return func.ST_Area(
            func.cast(cls.boundary, "geography")
        ) / 1_000_000

    @hybrid_property
    def perimeter_km(self) -> float:
        """Perimeter in km computed on a loaded instance."""
        geom = to_shape(self.boundary)
        projected = transform(_wgs84_to_equal_area, geom)
        return projected.length / 1_000             # m → km

    @perimeter_km.expression
    @classmethod
    def perimeter_km(cls):
        """SQL expression: ST_Perimeter on geography returns metres."""
        return func.ST_Perimeter(
            func.cast(cls.boundary, "geography")
        ) / 1_000

from sqlalchemy import select
from sqlalchemy.orm import Session

with Session(engine) as session:
    # Filter: area_sqkm > 50 compiles to a WHERE clause with ST_Area
    large_parcels = session.scalars(
        select(Parcel).where(Parcel.area_sqkm > 50)
    ).all()

    # Project specific columns — PostGIS computes both metrics in one query
    ranked = session.execute(
        select(Parcel.name, Parcel.area_sqkm, Parcel.perimeter_km)
        .order_by(Parcel.area_sqkm.desc())
        .limit(20)
    ).all()

    # Access on a loaded instance calls the Python method
    parcel = session.get(Parcel, 1)
    print(parcel.area_sqkm)   # uses to_shape + pyproj, no extra SQL

SELECT parcels.name,
       (ST_Area(CAST(parcels.boundary AS geography)) / 1000000) AS area_sqkm,
       (ST_Perimeter(CAST(parcels.boundary AS geography)) / 1000) AS perimeter_km
FROM parcels
ORDER BY area_sqkm DESC
LIMIT 20;

from sqlalchemy.orm import Session, load_only

with Session(engine) as session:
    # Eager-load boundary so the hybrid property works after session closes
    parcel = session.scalars(
        select(Parcel)
        .where(Parcel.id == 42)
        .options(load_only(Parcel.id, Parcel.name, Parcel.boundary))
    ).one()

    session.expunge(parcel)

# Safe to call after session close: boundary is already in memory
print(parcel.area_sqkm)

EXPLAIN (ANALYZE, BUFFERS, FORMAT TEXT)
SELECT name,
       ST_Area(boundary::geography) / 1000000 AS area_sqkm
FROM parcels
WHERE ST_Area(boundary::geography) / 1000000 > 50;

-- Functional index on the same expression as the hybrid property
CREATE INDEX idx_parcels_area_sqkm
    ON parcels ((ST_Area(boundary::geography) / 1000000));

ANALYZE parcels;

ALTER TABLE parcels
    ALTER COLUMN boundary SET STATISTICS 500;
ANALYZE parcels;

SELECT f_table_name, f_geometry_column, srid
FROM geometry_columns
WHERE f_table_name IN ('parcels', 'zones');

@area_sqkm.expression
@classmethod
def area_sqkm(cls):
    # Transform to EPSG:4326 first if the stored SRID could vary
    return func.ST_Area(
        func.ST_Transform(cls.boundary, 4326).__cast__("geography")
    ) / 1_000_000

from geoalchemy2.shape import to_shape

@hybrid_property
def area_sqkm(self) -> float:
    geom = to_shape(self.boundary)   # converts WKBElement → shapely.geometry.Polygon
    projected = transform(_wgs84_to_equal_area, geom)
    return projected.area / 1_000_000

# Option A: project scalars — geometry never crosses the wire
results = session.execute(
    select(Parcel.id, Parcel.name, Parcel.area_sqkm)
).all()

# Option B: stream full instances in batches of 500
for parcel in session.scalars(select(Parcel).yield_per(500)):
    process(parcel.area_sqkm)

-- 1. Update statistics
ANALYZE parcels;

-- 2. Force planner reconsideration in your session
SET enable_seqscan = off;
EXPLAIN (ANALYZE, BUFFERS)
SELECT name FROM parcels
WHERE (ST_Area(boundary::geography) / 1000000) > 50;

-- 3. Re-enable after testing
SET enable_seqscan = on;

from sqlalchemy import select
from sqlalchemy.dialects import postgresql

stmt = select(Parcel).where(Parcel.area_sqkm > 50)
compiled = stmt.compile(dialect=postgresql.dialect())
sql_text = str(compiled)

assert "ST_Area" in sql_text,      "ST_Area missing from compiled SQL"
assert "geography" in sql_text,    "geography cast missing"
assert "1000000" in sql_text,      "unit conversion missing"
print("Compilation OK:", sql_text[:120])

import pytest
from sqlalchemy.orm import Session

def test_area_consistency(engine, sample_parcel_id):
    """Instance value and SQL value should agree within floating-point tolerance."""
    with Session(engine) as session:
        parcel = session.get(Parcel, sample_parcel_id)
        instance_area = parcel.area_sqkm

        sql_area = session.scalar(
            select(Parcel.area_sqkm).where(Parcel.id == sample_parcel_id)
        )

    # Allow 0.5% tolerance: pyproj equal-area vs PostGIS spheroidal differ slightly
    assert abs(instance_area - sql_area) / sql_area < 0.005, (
        f"Instance ({instance_area:.4f}) and SQL ({sql_area:.4f}) diverge by "
        f"more than 0.5%"
    )

-- After creating the functional index and running ANALYZE,
-- this should show Index Scan (not Seq Scan) for selective queries
EXPLAIN (ANALYZE, BUFFERS)
SELECT name, (ST_Area(boundary::geography) / 1000000) AS area_sqkm
FROM parcels
WHERE (ST_Area(boundary::geography) / 1000000) BETWEEN 1 AND 5
ORDER BY area_sqkm;

import asyncio
from sqlalchemy.ext.asyncio import AsyncSession

async def get_large_parcels(session: AsyncSession):
    result = await session.scalars(
        select(Parcel).where(Parcel.area_sqkm > 50)
    )
    return result.all()