SQLAlchemy \ GeoAlchemy support of PostGIS geometry data types

Question

In an ORM I would expect to be able to do this:

session.add(Lake("hello",Polygon([(0,0),(1,0),(1,1)])))
lake = session.get(Lake).first()
assert isinstance(lake.geometry, Polygon)
assert isinstance(lake.geometry.get_exterior_ring().get_points()[0], Point)
print(lake.geometry.get_exterior_ring().get_points()[0].x)

Instead the only way I see to access the Points of my lake is through a rather complex monster of code:

ring = func.ST_ExteriorRing(func.ST_GeomFromWKB(Lake.geometry))
node_count = func.ST_NPoints(ring)
node_series = func.generate_series(1, node_count)
node_n = func.ST_PointN(ring, node_series)
node_n_x = func.ST_X(node_n)
node_n_y = func.ST_Y(node_n)
rows = session.query(Lake, node_n_x, node_n_y).all()
lake_coasts = {}
for row in rows:
    lake = row[0]
    if not lake in lake_coasts:
        lake_coasts[lake] = []
    lake_coast = lake_coasts[lake]
    lake_coast.append((row[1], row[2]))
for lake in lake_coasts:
    lake_coast = lake_coasts[lake]
    print("Lake #{0}: \"{1}\" has its coasts at {2}"
          .format(lake.id, lake.name, lake_coast))

and while this loop gets me the coordinates I want, I'm lost figuring out how to implement some Lake.get_coast() that returns this instance's coast.

Also I gave up on implementing the same for ""Lake""s with MULTIPOLYGONs as the nesting to get down to the points was too much for postgis (at least that's how I read the Error message)

I'm new to postgis, gis, python and sqla but in two days of googling I could not find anything that looks like an ORM within SQL-Alchemy 2 but only some SQL-helper functions (postgis) to parse the WKB but only within the database. Do I need some additional framework? I saw gdal, ogr, fiona but I feel like looking in the wrong direction. Is there some open source project using SQL-Alchemy2 with a nice structure? Entities, DAOs, whatever? How do you actually use this monster beyond this minimalist example?

Answer 1

You may use GeoAlchemy2 and Shapely:

    from geoalchemy2 import Geography
    from geoalchemy2.shape import from_shape
    from shapely.geometry import Polygon, Point, shape
    from sqlalchemy import create_engine, Column, Integer, func
    from sqlalchemy.ext.declarative import declarative_base
    from sqlalchemy.orm import scoped_session
    from sqlalchemy.orm import sessionmaker

    Base = declarative_base()
    db_session = scoped_session(sessionmaker())
    db_session.configure(bind=create_engine(connection_string))

    """
    create table polygon (
        id serial not null primary key,
        polygon geography(POLYGON) not null
    );
    """

    class TablePolygon(Base):
        __tablename__ = 'polygon'

        id = Column(Integer, primary_key=True)
        polygon = Column(Geography('POLYGON'), nullable=False)

    db_session.add(
        TablePolygon(
            id=1,
            polygon=from_shape(
                Polygon([
                    [-74.0133476, 40.7013069], [-73.9776421, 40.7121134], [-74.012661, 40.7230482],
                    [-74.0133476, 40.7013069]
                ])
            )
        )
    )

    geojson = """{"coordinates": [[[-73.9597893,40.8024021],[-73.9846802,40.7668997],[-73.9726639,40.7623468],[-73.9460564,40.7969415],[-73.9597893,40.8024021]]],"type": "Polygon"}"""

    db_session.add(
        TablePolygon(
            id=2,
            polygon=from_shape(
                shape(json.loads(geojson))
            )
        )
    )

    def find_polygon_id(point: Point):
        return db_session.query(
            TablePolygon.id
        ).filter(
            func.ST_Intersects(TablePolygon.polygon, from_shape(point))
        ).scalar()

    print(find_polygon_id(Point(-73.9822769, 40.7564925)))  # prints "None"
    print(find_polygon_id(Point(-73.9625359, 40.7858887)))  # prints "2"
    print(find_polygon_id(Point(-74.0059662, 40.7138058)))  # prints "1"