Implement native SQL Alchemy MVT

docs/source/data-publishing/ogcapi-tiles.rst

@@ -140,9 +140,6 @@ MVT-postgresql
 .. note::
    Must have PostGIS installed with protobuf-c support 
 
-.. note::
-   Geometry must be using EPSG:4326
-
 This provider gives support to serving tiles generated using `PostgreSQL <https://www.postgresql.org/>`_ with `PostGIS <https://postgis.net/>`_.
 The tiles are rendered on-the-fly using `ST_AsMVT <https://postgis.net/docs/ST_AsMVT.html>`_ and related methods.
 
@@ -152,24 +149,24 @@ This code block shows how to configure pygeoapi to render Mapbox vector tiles fr
 
    providers:
        - type: tile
-           name: MVT-postgresql
-           data:
-               host: 127.0.0.1
-               port: 3010 # Default 5432 if not provided
-               dbname: test
-               user: postgres
-               password: postgres
-               search_path: [osm, public]
-           id_field: osm_id
-           table: hotosm_bdi_waterways
-           geom_field: foo_geom
-           options:
-               zoom:
-                   min: 0
-                   max: 15
-           format:
-               name: pbf
-               mimetype: application/vnd.mapbox-vector-tile
+         name: MVT-postgresql
+         data:
+             host: 127.0.0.1
+             port: 3010 # Default 5432 if not provided
+             dbname: test
+             user: postgres
+             password: postgres
+             search_path: [osm, public]
+         id_field: osm_id
+         table: hotosm_bdi_waterways
+         geom_field: foo_geom
+         options:
+             zoom:
+                 min: 0
+                 max: 15
+         format:
+             name: pbf
+             mimetype: application/vnd.mapbox-vector-tile
 
 PostgreSQL-related connection options can also be added to `options`. Please refer to the :ref:`PostgreSQL OGC Features Provider<PostgreSQL>` documentation for more information.
 

pygeoapi/api/tiles.py

@@ -246,7 +246,7 @@ def get_collection_tiles_data(
         if content is None:
             msg = 'identifier not found'
             return api.get_exception(
-                HTTPStatus.NO_CONTENT, headers, format_, 'NocContent', msg)
+                HTTPStatus.NO_CONTENT, headers, format_, 'NoContent', msg)
         else:
             return headers, HTTPStatus.OK, content
 

pygeoapi/provider/mvt_postgresql.py

@@ -3,10 +3,12 @@
 # Authors: Prajwal Amaravati <prajwal.s@satsure.co>
 #          Tanvi Prasad <tanvi.prasad@cdpg.org.in>
 #          Bryan Robert <bryan.robert@cdpg.org.in>
+#          Benjamin Webb <bwebb@lincolninst.edu>
 #
 # Copyright (c) 2025 Prajwal Amaravati
 # Copyright (c) 2025 Tanvi Prasad
 # Copyright (c) 2025 Bryan Robert
+# Copyright (c) 2025 Benjamin Webb
 #
 # Permission is hereby granted, free of charge, to any person
 # obtaining a copy of this software and associated documentation
@@ -34,20 +36,20 @@
 from copy import deepcopy
 import logging
 
-from sqlalchemy.sql import text
-
+from sqlalchemy.sql import func, select
+from sqlalchemy.orm import Session
 from pygeoapi.models.provider.base import (
     TileSetMetadata, TileMatrixSetEnum, LinkType)
 from pygeoapi.provider.base import ProviderConnectionError
 from pygeoapi.provider.base_mvt import BaseMVTProvider
 from pygeoapi.provider.sql import PostgreSQLProvider
 from pygeoapi.provider.tile import ProviderTileNotFoundError
-from pygeoapi.util import url_join
+from pygeoapi.util import url_join, get_crs_from_uri
 
 LOGGER = logging.getLogger(__name__)
 
 
-class MVTPostgreSQLProvider(BaseMVTProvider):
+class MVTPostgreSQLProvider(PostgreSQLProvider, BaseMVTProvider):
     """
     MVT PostgreSQL Provider
     Provider for serving tiles rendered on-the-fly from
@@ -62,48 +64,38 @@ def __init__(self, provider_def):
 
         :returns: pygeoapi.provider.MVT.MVTPostgreSQLProvider
         """
-
-        super().__init__(provider_def)
-
         pg_def = deepcopy(provider_def)
         # delete the zoom option before initializing the PostgreSQL provider
         # that provider breaks otherwise
-        del pg_def["options"]["zoom"]
-        self.postgres = PostgreSQLProvider(pg_def)
-
-        self.layer_name = provider_def["table"]
-        self.table = provider_def['table']
-        self.id_field = provider_def['id_field']
-        self.geom = provider_def.get('geom_field', 'geom')
-
-        LOGGER.debug(f'DB connection: {repr(self.postgres._engine.url)}')
-
-    def __repr__(self):
-        return f'<MVTPostgreSQLProvider> {self.data}'
+        del pg_def['options']['zoom']
+        PostgreSQLProvider.__init__(self, pg_def)
+        BaseMVTProvider.__init__(self, provider_def)
 
-    @property
-    def service_url(self):
-        return self._service_url
+    def get_fields(self):
+        """
+        Get Postgres columns
 
-    @property
-    def service_metadata_url(self):
-        return self._service_metadata_url
+        :returns: `list` of columns
+        """
+        return [
+            c.label('id') if c.name == self.id_field else c
+            for c in self.table_model.__table__.columns
+            if c.name != self.geom
+        ]
 
     def get_layer(self):
         """
-        Extracts layer name from url
+        Use table name as layer name
 
-        :returns: layer name
+        :returns: `str` of layer name
         """
-
-        return self.layer_name
+        return self.table
 
     def get_tiling_schemes(self):
-
         return [
-                TileMatrixSetEnum.WEBMERCATORQUAD.value,
-                TileMatrixSetEnum.WORLDCRS84QUAD.value
-            ]
+            TileMatrixSetEnum.WEBMERCATORQUAD.value,
+            TileMatrixSetEnum.WORLDCRS84QUAD.value
+        ]
 
     def get_tiles_service(self, baseurl=None, servicepath=None,
                           dirpath=None, tile_type=None):
@@ -118,13 +110,14 @@ def get_tiles_service(self, baseurl=None, servicepath=None,
         :returns: `dict` of item tile service
         """
 
-        super().get_tiles_service(baseurl, servicepath,
-                                  dirpath, tile_type)
+        BaseMVTProvider.get_tiles_service(self,
+                                          baseurl, servicepath,
+                                          dirpath, tile_type)
 
         self._service_url = servicepath
         return self.get_tms_links()
 
-    def get_tiles(self, layer=None, tileset=None,
+    def get_tiles(self, layer='default', tileset=None,
                   z=None, y=None, x=None, format_=None):
         """
         Gets tile
@@ -141,64 +134,55 @@ def get_tiles(self, layer=None, tileset=None,
         if format_ == 'mvt':
             format_ = self.format_type
 
-        fields_arr = self.postgres.get_fields().keys()
-        fields = ', '.join(['"' + f + '"' for f in fields_arr])
-        if len(fields) != 0:
-            fields = ',' + fields
-
-        query = ''
-        if tileset == TileMatrixSetEnum.WEBMERCATORQUAD.value.tileMatrixSet:
-            if not self.is_in_limits(TileMatrixSetEnum.WEBMERCATORQUAD.value, z, x, y): # noqa
-                raise ProviderTileNotFoundError
-
-            query = text("""
-                WITH
-                    bounds AS (
-                        SELECT ST_TileEnvelope(:z, :x, :y) AS boundgeom
-                    ),
-                    mvtgeom AS (
-                        SELECT ST_AsMVTGeom(ST_Transform(ST_CurveToLine({geom}), 3857), bounds.boundgeom) AS geom {fields}
-                        FROM "{table}", bounds
-                        WHERE ST_Intersects({geom}, ST_Transform(bounds.boundgeom, 4326))
-                    )
-                SELECT ST_AsMVT(mvtgeom, 'default') FROM mvtgeom;
-            """.format(geom=self.geom, table=self.table, fields=fields)) # noqa
-
-        if tileset == TileMatrixSetEnum.WORLDCRS84QUAD.value.tileMatrixSet:
-            if not self.is_in_limits(TileMatrixSetEnum.WORLDCRS84QUAD.value, z, x, y): # noqa
-                raise ProviderTileNotFoundError
-
-            query = text("""
-                WITH
-                    bounds AS (
-                        SELECT ST_TileEnvelope(:z, :x, :y,
-                                'SRID=4326;POLYGON((-180 -90,-180 90,180 90,180 -90,-180 -90))'::geometry) AS boundgeom
-                    ),
-                    mvtgeom AS (
-                        SELECT ST_AsMVTGeom(ST_CurveToLine({geom}), bounds.boundgeom) AS geom {fields}
-                        FROM "{table}", bounds
-                        WHERE ST_Intersects({geom}, bounds.boundgeom)
-                    )
-                SELECT ST_AsMVT(mvtgeom, 'default') FROM mvtgeom;
-            """.format(geom=self.geom, table=self.table, fields=fields)) # noqa
-
-        with self.postgres._engine.connect() as session:
-            result = session.execute(query, {
-                    'z': z,
-                    'y': y,
-                    'x': x
-            }).fetchone()
-
-            if len(bytes(result[0])) == 0:
-                return None
-            return bytes(result[0])
+        [tileset_schema] = [
+            schema for schema in self.get_tiling_schemes()
+            if tileset == schema.tileMatrixSet
+        ]
+        if not self.is_in_limits(tileset_schema, z, x, y):
+            LOGGER.warning(f'Tile {z}/{x}/{y} not found')
+            return ProviderTileNotFoundError
+
+        storage_srid = get_crs_from_uri(self.storage_crs).to_string()
+        out_srid = get_crs_from_uri(tileset_schema.crs).to_string()
+        envelope = func.ST_TileEnvelope(z, x, y).label('bounds')
+
+        geom_column = getattr(self.table_model, self.geom)
+        geom_filter = geom_column.intersects(
+            func.ST_Transform(envelope, storage_srid)
+        )
+
+        mvtgeom = (
+            func.ST_AsMVTGeom(
+                func.ST_Transform(func.ST_CurveToLine(geom_column), out_srid),
+                func.ST_Transform(envelope, out_srid))
+            .label('mvtgeom')
+        )
+
+        mvtrow = (
+            select(mvtgeom, *self.fields)
+            .filter(geom_filter)
+            .cte('mvtrow')
+            .table_valued()
+        )
+
+        mvtquery = select(
+            func.ST_AsMVT(mvtrow, layer)
+        )
+
+        with Session(self._engine) as session:
+            result = bytes(
+                session.execute(mvtquery).scalar()
+            ) or None
+
+        return result
 
     def get_html_metadata(self, dataset, server_url, layer, tileset,
                           title, description, keywords, **kwargs):
 
         service_url = url_join(
             server_url,
-            f'collections/{dataset}/tiles/{tileset}/{{tileMatrix}}/{{tileRow}}/{{tileCol}}?f=mvt')  # noqa
+            f'collections/{dataset}/tiles/{tileset}'
+            '{tileMatrix}/{tileRow}/{tileCol}?f=mvt')
         metadata_url = url_join(
             server_url,
             f'collections/{dataset}/tiles/{tileset}/metadata')
@@ -217,9 +201,10 @@ def get_default_metadata(self, dataset, server_url, layer, tileset,
 
         service_url = url_join(
             server_url,
-            f'collections/{dataset}/tiles/{tileset}/{{tileMatrix}}/{{tileRow}}/{{tileCol}}?f=mvt')  # noqa
+            f'collections/{dataset}/tiles/{tileset}',
+            '{tileMatrix}/{tileRow}/{tileCol}?f=mvt'
+            )
 
-        content = {}
         tiling_schemes = self.get_tiling_schemes()
         # Default values
         tileMatrixSetURI = tiling_schemes[0].tileMatrixSetURI
@@ -231,17 +216,18 @@ def get_default_metadata(self, dataset, server_url, layer, tileset,
                 tileMatrixSetURI = schema.tileMatrixSetURI
 
                 tiling_scheme_url = url_join(
-                    server_url, f'/TileMatrixSets/{schema.tileMatrixSet}')
-                tiling_scheme_url_type = "application/json"
+                    server_url, 'TileMatrixSets', schema.tileMatrixSet)
+                tiling_scheme_url_type = 'application/json'
                 tiling_scheme_url_title = f'{schema.tileMatrixSet} tile matrix set definition' # noqa
 
-                tiling_scheme = LinkType(href=tiling_scheme_url,
-                                         rel="http://www.opengis.net/def/rel/ogc/1.0/tiling-scheme", # noqa
-                                         type_=tiling_scheme_url_type,
-                                         title=tiling_scheme_url_title)
+                tiling_scheme = LinkType(
+                    href=tiling_scheme_url,
+                    rel='http://www.opengis.net/def/rel/ogc/1.0/tiling-scheme',
+                    type_=tiling_scheme_url_type,
+                    title=tiling_scheme_url_title)
 
         if tiling_scheme is None:
-            msg = f'Could not identify a valid tiling schema'  # noqa
+            msg = 'Could not identify a valid tiling schema'
             LOGGER.error(msg)
             raise ProviderConnectionError(msg)
 
@@ -250,9 +236,9 @@ def get_default_metadata(self, dataset, server_url, layer, tileset,
                                   tileMatrixSetURI=tileMatrixSetURI)
 
         links = []
-        service_url_link_type = "application/vnd.mapbox-vector-tile"
+        service_url_link_type = 'application/vnd.mapbox-vector-tile'
         service_url_link_title = f'{tileset} vector tiles for {layer}'
-        service_url_link = LinkType(href=service_url, rel="item",
+        service_url_link = LinkType(href=service_url, rel='item',
                                     type_=service_url_link_type,
                                     title=service_url_link_title)
 
@@ -261,4 +247,7 @@ def get_default_metadata(self, dataset, server_url, layer, tileset,
 
         content.links = links
 
-        return content.dict(exclude_none=True)
+        return content.model_dump(exclude_none=True)
+
+    def __repr__(self):
+        return f'<MVTPostgreSQLProvider> {self.data}'

pygeoapi/provider/tile.py

@@ -56,7 +56,6 @@ def __init__(self, provider_def):
         self.mimetype = provider_def['format']['mimetype']
         self.options = provider_def.get('options')
         self.tile_type = None
-        self.fields = {}
 
     def get_layer(self):
         """

pygeoapi/templates/collections/tiles/index.html

@@ -108,7 +108,7 @@ <h3>Tiles</h3>
             maxZoom: {{ data['maxzoom'] }},
             indexMaxZoom: 5,
             getFeatureId: function(feat) {
-                return feat.properties["id"]
+                return feat.properties.id || feat.properties.fid || feat.properties.uri;
             }
         };
 
@@ -129,7 +129,7 @@ <h3>Tiles</h3>
                     .openOn(map);
 
                 clearHighlight();
-                highlight = e.layer.properties.id;
+                highlight = e.layer.properties.id || e.layer.properties.fid || e.layer.properties.uri;
                 tilesPbfLayer.setFeatureStyle(highlight, {
                     weight: 2,
                     color: 'red',