SQLAlchemy \ GeoAlchemy对PostGIS几何数据类型的支持

Question

在ORM中，我希望能够做到这一点：

(湖(“你好”，多边形((0，0)，(1,0)，(1，1)session.add=session.get(session.get).first()断言实例(lake.geometry，Polygon)断言session.add点)

相反，我看到的访问湖点的唯一方法是通过一个相当复杂的代码怪物：

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( node_n_y).all() lake_coasts = {}表示行中的行: lake_coasts: lake_coastslake = [] lake_coast = lake_coastslake lake_coast.append((row1，row2))用于lake_coasts: lake_coast = lake_coastslake print(“Lake#{0}：\{1}”在{2}“.format”(lake.id，lake.name )，( lake_coast))

虽然这个循环获得了我想要的坐标，但是我弄不清楚如何实现一些返回这个实例的Lake.get_coast()。

此外，我也放弃了用MULTIPOLYGON实现“Lake”的方法，因为对postgis来说，嵌套的时间太长了(至少我是这样读取错误消息的)。

我对postgis、gis、python和sqla很陌生，但在谷歌搜索的两天内，我找不到任何在sqla 2中看起来像ORM的东西，而只是一些用于解析WKB的SQL函数(postgis)，而在数据库中只找到了。我需要一些额外的框架吗？我看到了gdal，ogr，fiona，但是我觉得我看错了方向。是否有一些使用结构良好的SQL 2的开源项目？实体道斯什么的？在这个极简主义的例子之外，你如何使用这个怪物呢？

Answer 1

您可以使用GeoAlchemy2和修长

    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"