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?

First, SQLAlchemy is actually two separate pieces, the Core and the ORM. The former is just an abstraction over SQL. And the latter is not actually a traditional ORM that automatically statically maps classes to tables; it's a framework for building data mappers between classes and a relational algebra engine (like the SQLAlchemy Core). — abarnert, May 23 '13 at 01:17
You should probably read [Key Features of SQLAlchemy](http://www.sqlalchemy.org/features.html) and [Object Relational Tutorial](http://docs.sqlalchemy.org/en/rel_0_8/orm/tutorial.html), if you haven't yet, and then read through some of the examples. — abarnert, May 23 '13 at 01:19
I read those pages before but in this case it actually seams to be an GEOAlchemy issue. With Strings and numbers the ORM magic works as expected but the GIS stuff is crazy. session.query(Office) results in a real Office object with all its workers and what not but with geometries it doesn't work that way. — Giszmo, May 23 '13 at 05:42
I've never used GEOAlchemy. You can always build mappers yourself to connect the geometry objects to the relational engine. But hopefully you'll get a better answer. — abarnert, May 23 '13 at 17:24
I'm currently working on a parser to get in and out of the binary protocol. Guess it simply doesn't exist. I'm planning to share my result on GitHub. Thank you abarnert anyway for your reply. — Giszmo, May 25 '13 at 02:39
You can use ST_GeometryN to access the polygons within a MultiPolygon, see http://geoalchemy-2.readthedocs.org/en/0.2.4/spatial_functions.html — John Powell, May 31 '14 at 17:40

score 1 · Answer 1 · answered Feb 12 '22 at 09:38

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"

SQLAlchemy \ GeoAlchemy support of PostGIS geometry data types

1 Answers1